Scaling relations in quasi-static magnetoconvection with a strong vertical magnetic field

Abstract

The scaling law for the horizontal length scale $\ell$ relative to the domain height $L$, originating from the linear theory of quasi-static magnetoconvection, $\ell /L \sim Q^{-1/6}$, has been verified through two-dimensional (2-D) direct numerical simulation (DNS), particularly at high values of the Chandrasekhar number ($Q$). This relationship remains valid within a specific flow regime characterized by columnar structures aligned with the magnetic field. Expanding upon the $Q$-dependence of the horizontal length scale, we have derived scaling laws for the Nusselt number ($Nu$) and the Reynolds number ($Re$) as functions of the driving forces (Rayleigh number ($Ra$) and $Q$) in quasi-static magnetoconvection influenced by a strong magnetic field. These scaling relations, $Nu \sim Ra/Q$ and $Re \sim Ra Q^{-5/6}$, have been successfully validated using 2-D DNS data spanning a wide range of five decades in $Q$, ranging from $10^5$ to $10^9$. The successful validation of the relations at large $Q$ values, combined with our theoretical analysis of dissipation rates and the incorporation of the horizontal length scale's influence on scaling behaviour, presents a valid approach for deriving scaling laws under various conditions.

1. Introduction

The convective motion of electrically conducting fluids in the presence of an externally imposed magnetic field plays a crucial role in both natural phenomena and practical flows. In astrophysical environments, fluid motions are intricately coupled with magnetic fields, such as in the outer layers of the Sun and other late-type stars. The interaction between thermal convection and magnetic fields in these scenarios has been extensively studied (Proctor & Weiss 1982; Cattaneo, Emonet & Weiss 2003; Schüssler & Vögler 2006). In industrial systems, magnetic damping is commonly employed in metallurgical applications, while magnetic levitation, pumping and heating of liquid metal are vital aspects of nuclear engineering (Davidson 2001). Despite its prevalence, understanding the dynamics of magnetoconvection in extreme parameter regimes, particularly in geophysical and astrophysical systems, remains a formidable challenge. Therefore, it is crucial to identify relevant flow regimes and establish corresponding scaling relations to characterize these complex flows. By discerning the points at which flow transitions occur, we can carefully extrapolate scaling laws to a reasonable range of input parameters, thereby providing practical utility in understanding and predicting these phenomena.

The planar Rayleigh–Bénard convection (RBC) is a useful model for investigating the fundamental dynamics of magnetic fields on convection. The canonical RBC consists of a fluid layer confined between a pair of plane, parallel boundaries separated by a distance $L$. The boundaries are maintained at a constant temperature difference $\varDelta$. In the absence of a magnetic field, the Rayleigh number $Ra={g\beta \Delta L^3}/{\nu \kappa }$ and the Prandtl number $Pr={\nu }/{\kappa }$ characterize the thermal forcing and the ratio of the diffusivities present in the system, where $\beta$ is the coefficient of thermal expansion, $g$ the acceleration due to gravity, with $\nu$ and $\kappa$ being the momentum and thermal diffusivity, respectively. In the presence of an external magnetic field of strength $B_0$, the damping effect by virtue of the Lorentz force must be considered. The strength of the Lorentz force to the viscosity is quantified by the Chandrasekhar number $Q= {B^2_0L^2}/({\rho \nu \mu \eta })$, where $\rho$ is the density, $\mu$ the magnetic permeability and $\eta$ magnetic diffusivity.

Early studies on magnetoconvection, both experimental as well as theoretical (Thompson 1951; Jirlow 1956; Nakagawa 1957; Chandrasekhar 1961; Busse & Clever 1982), are primarily devoted to the linear stability analysis. In recent decades, much focus has been on the scaling relations between the dimensionless heat transport as a function of $Ra$ and $Q$. The dimensionless convective heat transport is characterized by the Nusselt number $Nu$, the ratio of the total heat flux to the conduction heat transfer. In recent decades, experiments by Cioni, Chaumat & Sommeria (2000), Aurnou & Olson (2001), Burr & Müller (2001) and King & Aurnou (2013) have suggested different scalings. Aurnou & Olson (2001) found $Nu \sim (Ra/Q)^{1/2}$ in contrast with the findings of Burr & Müller (2001) who report a scaling law $Nu \sim (Ra/Q)^{2/3}$. Both of these studies consider $0 < Q \lesssim 10^4$. Yu, Zhang & Ni (2018) conducted numerical simulations over a similar range of values of $Q$, and found the scaling relations between $Nu$ and $Ra/Q$ are in good agreement with the empirical formulae obtained in the experimental studies of Aurnou & Olson (2001). Experimental investigations by Cioni et al. (2000) pushed the range of Chandrasekhar number by two decades up to $Q \approx 10^6$, and comprehensively characterized two major regimes influenced by the magnetic effects, with an additional transitional regime between them. The first regime, interpreted by them as a condition of marginal stability for the thermal boundary layer, the heat transport scaling $Nu \sim Ra/Q$ holds. For the second major regime, which is characterized by higher thermal forcing, they report $Nu \sim Ra^{0.43}Q^{-0.25}$. Recently, direct numerical simulations (DNSs) of quasi-static magnetoconvection by Yan et al. (2019) with magnetic strengths up to $Q=10^8$ have shown the scaling laws of Aurnou & Olson (2001), Burr & Müller (2001) and Yu et al. (2018) to be limited to low values of $Q$. The proposed convective regimes are distinguished by flow characteristics, with the first regime being reminiscent of the linear convection, characterized by laminar, cellular structures. The heat transport scaling law for the second regime, characterized by the existence of quasi-laminar columnar structures, is suggested as $Nu \sim (Ra/Q)^{\gamma }$, where $\gamma \rightarrow 1$ as $Q \rightarrow \infty$. Based on an exponential fit for $\gamma$, a value of $0.95$ is suggested at $Q=10^{16}$. However, as the authors mention, it must be noted that the exponential fit has no physical basis and is shown only to provide a guide for the behaviour at large values of $Q$.

Based on the scaling approach introduced by Grossmann & Lohse (2000), Zürner et al. (2016) decomposed the dissipation rates in magnetoconvection into their bulk and boundary layer counterparts and identified four distinct regimes distinguished by the strength of the external magnetic field and the level of turbulence in the flow. Studies by Zürner et al. (2020) and Akhmedagaev et al. (2020) report a detailed analysis of the spatial structure of magnetoconvective flows in cylindrical geometry with no-slip walls. Wall modes similar to those in RBC with rotation are observed near the linear stability limit. Furthermore, in these studies, the scaling relations for normalized Nusselt and Reynolds numbers reveal that the global transport properties approach a universal power law at larger degrees of supercriticality.

Following Grossmann & Lohse (2000), Bhattacharyya (2006) obtained $Nu \sim Nu(Ra,Q)$ scaling relations for different regimes characterized by weak and strong magnetic fields. It must be noted here that Bhattacharyya (2006) and Zürner et al. (2016) assumed the length scale associated with the magnetic induction to be the domain height $L$. The presence of a magnetic field opposes fluid motions perpendicular to the field lines while leaving the parallel component unopposed, thus the classical RBC under the effect of an externally applied magnetic field becomes highly anisotropic. Due to anisotropy, the length scale of the flow is significantly altered; the scaling of which remains vital in the scaling laws associated with the heat and momentum transport and Ohmic dissipation. The primary objective of our study is to explore the impact of the horizontal length scale on the scaling relations previously discussed, focusing specifically on the framework of two-dimensional (2-D) quasi-static assumptions. Although the real-world applications of magnetoconvection are three-dimensional (3-D), 2-D simulations are substantially cheaper in terms of computational time and thus can be utilized to understand certain aspects of the flow, including scalings. Goluskin et al. (2014) and Wang et al. (2020) are some noteworthy 2-D studies that have explored the ability of the convection to drive vertically sheared, large-scale horizontal flow. Building on these works, recently, quasi-2-D simulations of magnetoconvection have been used to study the magnetic damping of jet flows which finds applications in probing the mechanisms relevant to damping of large-scale azimuthally directed jets on Jupiter (Aggarwal, Aurnou & Horn 2022). Here, our 2-D investigation extends beyond previous studies by considering even larger values of $Q$, reaching up to $10^9$. This range of values is particularly relevant for geophysical and astrophysical systems, such as the Earth's outer core, where $Q$ is estimated to be of the order of $10^{16}$ (Gillet et al. 2010). By identifying scaling relations that converge at values of $Q$ approaching realistic estimates, we aim to enhance our understanding of the fundamental aspects of convection in these systems.

2. Physical model and numerical settings

2.1. Governing equations

The dimensionless equations governing the flow of a conducting fluid in Boussinesq approximation, under the action of a quasi-static, vertically imposed external magnetic field $\boldsymbol {B}= B_0 \boldsymbol {\hat {e}_x}$ are

(2.1a)\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{u} ={-}\boldsymbol{\nabla} p + Q\sqrt{\frac{Pr}{Ra}} (\boldsymbol{J} \times \boldsymbol{\hat{e}_x}) + \theta \boldsymbol{\hat{e}_x} + \sqrt{\frac{Pr}{Ra}} \nabla^2 \boldsymbol{u}, \end{gather}

(2.1b)\begin{gather}\frac{\partial {\theta}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \theta = \frac{1}{\sqrt{RaPr}} \nabla^2 {\theta}, \end{gather}

(2.1c)\begin{gather}\boldsymbol{J} ={-}\boldsymbol{\nabla} \varPhi + \boldsymbol{u} \times \boldsymbol{\hat{e}_x}, \end{gather}

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

(2.1e)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{J} = 0, \end{gather}

where $\boldsymbol {u}$, $p$ and $\theta$ represent the velocity, pressure and temperature fields. The vertical unit vector antiparallel to the direction of gravity is $\boldsymbol {\hat {e}_x}$, and $\varPhi$ and $\boldsymbol {J}$ represent the electric potential and the current density, respectively. At all places in this paper, $x$ and $y$ represent the vertical and horizontal directions, respectively. The input parameters $Ra$, $Pr$ and $Q$ are defined in § 1. The rest of the symbols have their usual meanings. In (2.1), domain height $L$, free-fall velocity $u_{ff} = {(\beta g \Delta L)}^{1/2}$, free-fall time scale $t_{ff} = L/u_{ff}$, imposed temperature difference $\varDelta$, free-fall velocity-based dynamic pressure $\rho _0 u^2_{ff}$ and the magnitude of the imposed magnetic field $B_0$, are adopted for non-dimensionalization. All lateral boundaries are set to periodic. At the top and bottom walls, stress-free boundary conditions for $\boldsymbol {u}$, isothermal for $\theta$ and an insulating boundary condition on $\varPhi$ are imposed, $u_x = {\partial u_y}/{\partial x} = 0 \ \text {at} \ x = 0,1; \quad \theta = 1(0) \text { at } x = 0(1); \quad {\partial \varPhi }/{\partial x} = 0 \ \text {at} \ x=0, 1$. To compute $\boldsymbol {J}$, the constraint of charge conservation, given by (2.1e), is invoked to obtain a Poisson equation for $\varPhi$, $\nabla ^2 \varPhi = \boldsymbol {\nabla } \boldsymbol {\cdot } (\boldsymbol {u} \times \boldsymbol {\hat {e}_x}).$

The Reynolds number, $Re$, referred to in the subsequent sections, is defined as $Re = {\mathcal {U}L}/{\nu }$, where $\mathcal {U}$ is the characteristic velocity. Depending on if the characteristic velocity is $u_{rms} = ( \overline {\langle u^2_x \rangle _A + \langle u^2_y \rangle _A} )^{1/2}$, $u_{x,rms}= ( \overline {\langle u^2_x \rangle _A } )^{1/2}$ or $u_{y,rms} = (\overline {\langle u^2_y \rangle _A } )^{1/2}$, we further define $Re$, $Re_x$ and $Re_y$ to characterize the total, vertical and horizontal Reynolds numbers, respectively. Here $\langle {\cdot } \rangle _A$ and $\overline {({\cdot })}$ denote the area-average over the whole 2-D domain, and time-averaging, respectively.

2.2. Numerics

A conservative, second-order centred spatial discretization is employed on a staggered grid (van der Poel et al. 2015). Time marching is performed with a fractional-step third-order Runge–Kutta scheme, in combination with the Crank–Nicholson scheme for the implicit terms. In the following section, the implementation of the subroutines for $\varPhi$ and the validation of the code are described in detail.

2.2.1. Code validation

In this work, we have used AFiD (van der Poel et al. 2015; Zhu et al. 2018) as the base solver, and implemented the subroutines for solving a Poisson equation for the electric potential $\varPhi$, and an explicit Lorentz body force term. Temperature and velocities are computed on cell faces, with $\theta$ and $u_x$ located at the same face (to avoid interpolation errors in calculating the term $\theta \boldsymbol {\hat {e}_x}$). The scalar variables $p$ and $\varPhi$ are computed at the cell centres. Using the cell-centred values of $\varPhi$, the components of the current density are then constructed on the cell faces, which are used to finally compute the components of the Lorentz force.

The solver was validated by comparing the Nusselt number with the values reported in Yan et al. (2019) at different Rayleigh numbers away from the onset. Table 1 shows the Nusselt number from 3-D test runs at $Q=10^6$ to reproduce results from the literature. The resolutions employed for these cases in order, from smallest to largest $Ra$, are $(96 \times 180 \times 180)$, $(256 \times 384 \times 384)$ and $(512 \times 1024 \times 1024)$, where the first entry corresponds to the grid number in the vertical $(x)$ direction. Here $Nu_p = \overline { \langle {\rm d} \theta /{{\rm d}\kern0.7pt x} \rangle }_{x=0}$ is the plate Nusselt number, $Nu_{\nu + \eta } = 1+ \sqrt {Ra\,Pr} ( \overline {\langle \epsilon _{\nu } \rangle } +\overline {\langle \epsilon _{\eta } \rangle } )$ and $Nu_{\kappa } = \sqrt {Ra\,Pr} \ \overline {\langle \epsilon _{\kappa } \rangle }$ represent the Nusselt numbers calculated from the global, time- and area-averaged viscous $(\overline {\langle \epsilon _{\nu } \rangle } = \sqrt {Pr/Ra} \ \overline {\langle |\boldsymbol {\nabla } \boldsymbol {u}|^2 \rangle } )$ and Ohmic dissipation $(\overline {\langle \epsilon _{\eta } \rangle } = Q\sqrt {Pr/Ra} \ \overline {\langle |\boldsymbol {\nabla } \boldsymbol {J}|^2 \rangle } )$ and thermal dissipation $(\overline {\langle \epsilon _{\kappa } \rangle } = \sqrt {Pr/Ra} \ \overline {\langle |\boldsymbol {\nabla } {\theta }|^2 \rangle } )$, respectively.

Table 1. Nusselt number from 3-D test runs at $Q=10^6$ to reproduce results from the literature.

In addition to that, the dominant length scale at $Ra \gtrapprox Ra_c$ (very close to the onset) was also calculated from the energy spectra of the vertical velocity, and validated against the predictions from the linear theory. In figure 1, the normalized kinetic energy spectra of the vertical velocity $E_{uu}(k) = \sum [\hat {u}_x(k) \hat {u}^*_x(k)$] computed at mid-depth are shown. Here $\hat {u}_x$ and $\hat {u}^*_x$ are the Fourier transform of $u_x$ and its complex conjugate, respectively. The dominant integer wavenumber represented by the highest peak of the spectrum, indicates the number of (near-onset) critical horizontal wavelengths ($\lambda _c$) present in the domain. The aspect ratio for $Q=10^3$, $10^5$ and $10^7$ cases was set to $\varGamma = 11.0 , 5.0$ and $2.25$, respectively, which admit $\approx 10\lambda _c$. For the $Q=10^6$ case, $\varGamma = 3.5$, which allows for $\approx 11\lambda _c$ to be resolved. The number of critical horizontal wavelengths observed in our simulations match reasonably well with the linear theory predictions for a given aspect ratio.

Figure 1. Normalized kinetic energy spectra of the vertical velocity at (a) $Q=10^3, Ra = 1.530 \times 10^4 \gtrapprox Ra_c = 1.521 \times 10^4$, (b) $Q=10^5, Ra = 1.088 \times 10^6 \gtrapprox Ra_c = 1.078 \times 10^6$, (c) $Q=10^6, Ra = 1.038 \times 10^7 \gtrapprox Ra_c = 1.028 \times 10^7$ and (d) $Q=10^7, Ra = 1.016 \times 10^8 \gtrapprox Ra_c = 1.006 \times 10^8$. Vertical dashed red lines show the results from linear theory.

3. Results and discussion

3.1. Parameter range and flow structure

We have considered five sets of simulations with $Q$ ranging from $10^5$ to $10^9$. For each set, at a fixed $Q$, the Rayleigh number is varied up to $Ra \gtrsim 100Ra_c$, $Ra_c$ being the critical Rayleigh number at a given $Q$, which is calculated from the expression of the Rayleigh number corresponding to the marginal stability of the horizontal wavenumber. For details, the interested reader is referred to Chandrasekhar (1961). In all of our simulations, we have used $Pr=1$. More details on grid resolution and the parameter ranges explored are listed in the Appendix.

Figure 2 shows the visualization of the effect of increasing $Ra$ at a fixed $Q$ ($Q=10^6$). The flow structures at all $Ra$ are stabilized by the Lorentz force, resulting in the convection directed along the externally imposed magnetic field. For the cases in figure 2(a–c), flow is dominated by cellular structures which eventually morph into columns as $Ra$ is increased. The transport of heat in the cellular regime is primarily via conduction. In the columnar regime, the length scale of the structures stays more or less the same, until a third regime of magnetoplumes emerges. Upon further increasing the Rayleigh number, the plumes begin to merge, characterized by stronger buoyancy relative to the Lorentz force. We will not be dealing with the scaling laws in this merging-plume regime. With increasing $Ra$, the thermal boundary layer thickness also decreases. This has important consequences on the heat transfer scaling within the columnar regime, which we will discuss later. The flow in cellular, columnar and plume regimes is strictly steady. Unsteadiness kicks in when the stabilizing effect of the Lorentz force is weakened at around $Ra/Ra_c \sim 80$, which also marks the onset of the merging plume regime. Figure 2(r) shows such a case. Upon comparisons with the flow structures in 3-D simulations, we observe that the intercolumn distance in 2-D simulations is identical in the horizontal direction. However, no such uniformity is seen in 3-D simulations of Yan et al. (2019). To illustrate the differences, in figure 3, we have presented instantaneous temperature contours computed in a 3-D simulation in the columnar regime. Given insulating boundary conditions $(\partial _x \varPhi |_{x=0,1}=0)$ (non-conducting walls), electric currents induced by the horizontal motion $(\boldsymbol {u} \times \hat {\boldsymbol {e}}_x)$ will close within the convective roll itself. In 3-D simulations, this gives rise to currents in the horizontal plane $(y$–$z)$ near the boundaries. In the case of 2-D simulations, since one horizontal component of velocity is zero ($u_z$ in the present case), the current associated with that component also ceases to exist, and only one of the components, $u_y$, remains. This gives rise to the Lorentz force which is confined to the only horizontal degree of freedom $(y)$ existing in the 2-D flow. Since everything is constrained to a single horizontal dimension, there is a tendency of the flow structures to adhere to a specific equidistant pattern, until the stabilizing effect of the Lorentz force is weakened at relatively higher $Ra$. At that point, the interplume distance starts to vary.

Figure 2. Instantaneous temperature contours in 2-D quasi-static magnetoconvection at $Q=10^6$; $Ra$ increases from $(\textit {a})$ to $(r)$, as listed in the Appendix for $Q =10^6$, from the smallest $(1.1\times 10^7)$ to the largest $(8\times 10^8)$. From $(\textit {a})$ to $(r)$, $Ra/Ra_c =$ $1.07$, $1.26$, $1.46$, $1.65$, $1.94$, $2.43$, $2.92$, $3.89$, $4.86$, $5.84$, $6.81$, $8.75$, $9.73$, $14.59$, $19.45$, $24.32$, $38.91$, $77.82$, respectively. The colourmap values have been adjusted to $[0.2,0.8]$ for better visibility.

Figure 3. Columnar structures in a 3-D simulation at $Ra = 1\times 10^8$, $Q=10^6$, $\varGamma =1.0$. Colour represents instantaneous temperature. In the 3-D simulation, these values lie in the columnar regime.

The anisotropy in the convective flow can be characterized by examining the Reynolds number based on the vertical and horizontal components of the root mean square (r.m.s.) velocity. Figure 4(e) clearly shows that the vertical component of the velocity dominates over the horizontal component. As $Ra$ is increased, the vertical component increases significantly relative to the horizontal component.

Figure 4. Vertical profiles of the (a) time- and horizontally averaged mean temperature $\overline {\langle \theta \rangle }_y$, (b) horizontally averaged r.m.s. fluctuations of temperature $\langle \theta _{rms} \rangle _y$, (c) viscous dissipation $\overline {\langle \epsilon _{\nu } \rangle }_y$, (d) Ohmic dissipation $\overline {\langle \epsilon _{\eta } \rangle }_y$ and (e) Reynolds number. Blue curves represent $Re_y$ and red curves represent $Re_x$; darker curves represent higher $Ra$. The data are plotted for $Q=10^6$. The Rayleigh numbers from smallest to largest are $1.1\times 10^7, 1.5\times 10^7, 2\times 10^7, 3\times 10^7, 5\times 10^7, 7\times 10^7$. The dashed black curve represents the vertical Reynolds number ($Re_x$) computed in a 3-D simulation at $Ra=7\times 10^7$. The aspect ratio for the 3-D simulation is $\varGamma _{3D}=1.0$. However, that is not expected to change the conclusion.

In figure 4(a,b), vertical profiles of time- and horizontally averaged mean temperature $\overline {\langle \theta \rangle }_y$, and horizontally averaged r.m.s. fluctuations of temperature $\langle \theta _{rms} \rangle _y$ at $Q=10^6$ are presented. For lower values of $Ra$, the r.m.s. fluctuations of temperature are small, and the mean temperature follows the conduction profile. At larger values of $Ra$, the mean temperature tends to be an isothermal profile in the bulk, and the near-wall peaks in the r.m.s. profiles indicate the presence of well-developed thermal boundary layers. The darkest shade in figure 4 represents $Ra=7\times 10^7$, which lies in the columnar regime (see figure 2k).

Figure 4(c,d) show the depth-dependence of the time- and horizontally averaged viscous $(\epsilon _{\nu } = (Pr/Ra)^{1/2} \langle \overline {|\boldsymbol {\nabla } \boldsymbol {u}|^2} \rangle _y)$ and Ohmic dissipation $(\epsilon _{\eta } = Q(Pr/Ra)^{1/2} \langle \overline {| \boldsymbol {J}|^2} \rangle _y)$. Since we have adopted stress-free boundary conditions, the viscous dissipation is larger in the bulk of the domain, away from the boundaries. Ohmic dissipation, on the other hand, is dominant near the boundaries and decreases in the interior of the domain. At larger $Ra$, an increase in the Ohmic dissipation is expected due to larger currents generated by enhanced convection.

3.2. Scaling analysis

Figure 5(a) gives a brief overview of the parameter space characterized by dimensionless heat transport. Three distinct scaling regimes with different slopes are observed at large values of $Q$: (i) a cellular regime close to the onset; (ii) $Nu \sim Ra$ regime, characterized by the existence of elongated, vertically aligned columns; and (iii) $Nu \sim Ra^{1/3}$ regime, which is characterized by the existence of plume-like structures with thin stems and broadened heads. Some select cases representing these regimes are shown in figure 6. The corresponding scaling regimes in $(Nu,Ra)$ space are represented by black dashed (columnar) and dotted (plume) lines in figure 5(a). It is important to note here that the states where transitions occur from cellular to columnar to plume regimes are different in 2-D simulations in comparison with 3-D simulations. This is because, for 2-D simulations, consistent with the mass conservation, we observed a higher $Re_x$ compared with that in 3-D simulations at the same $Ra$. To see a comparison of $Re_x$ in 2-D and 3-D simulations, refer to figure 4(e). Higher $Re_x$ is responsible for enhancing the buoyant fluxes in 2-D simulations, which ultimately translates into an increased $Nu$, at least for the regime that we are most interested in (the 2-D columnar regime), for example, see figure 5(a). The differing transition locations are evident. In addition to that, we can also clearly observe the columnar regime in both 2-D and 3-D simulations, albeit with slightly different $Nu$. The transition locations and the magnitude of the heat flux and flow velocity may differ; however, structural similarities exist, and for 2-D and 3-D simulations alike, the flow is characterized by cellular, columnar and plume-like structures for the considered $Ra$ range at different values of $Q$. This is important for scaling since a specific regime in 3-D simulations is expected and shown to have a counterpart in 2-D simulations.

Figure 5. (a) Nusselt number; columnar and plume regimes are marked by $Nu \propto Ra$ (dashed) and $Nu \propto Ra^{1/3}$ (dotted) lines for $Q=10^7$. (b) Thermal boundary layer thickness as a function of the supercriticality parameter. The 3-D data adapted from Yan et al. (2019) are shown for $Q = 10^5, 10^6, 10^7$ and $10^8$.

Figure 6. Visualization of the cellular, columnar and plume regimes at $Q=10^7$ with $Ra =1.1\times 10^8, 8\times 10^8~\text {and}~5\times 10^9$ from (a) to (c), respectively. The colourmap values have been adjusted to $[0.2,0.8]$ for better visibility.

The existence of the columnar and plume regimes can be explained very well by invoking the predictions from marginally stable thermal boundary layer analysis, which predicts $Nu \sim {L}/{\lambda _{\theta }} \sim Ra/Q$ at large values of $Q_{\lambda _{\theta }} = B^2_0 \lambda ^2_{\theta }/(\mu \rho \nu \eta ) = Q\lambda ^2_{\theta }/L^2$, and ${Nu \sim {L}/{\lambda _{\theta }} \sim Ra^{1/3}}$ at smaller values of $Q_{\lambda _{\theta }}$ due to thinner thermal boundary layer thickness in the plume regime (Bhattacharjee, Das & Banerjee 1991), see figure 5(b).

In the following sections, for the columnar regime, we will explain the scaling laws concerning not only the heat transport, but also the horizontal length scale, flow velocity, and Ohmic dissipation by using energetic arguments. To begin, we start with the exact relations for time- and volume-averaged (area-averaged in 2-D simulations) viscous dissipation rate $\epsilon _{\nu } = \nu \langle \overline {|\boldsymbol {\nabla } \boldsymbol {u}|^2} \rangle$, Ohmic dissipation rate $\epsilon _{\eta } = (\rho \sigma )^{-1} \langle \overline {|\boldsymbol {J}|^2} \rangle$ and the thermal dissipation rate $\epsilon _{\theta } = \kappa \langle \overline {|\boldsymbol {\nabla } \theta |^2} \rangle$ (Grossmann & Lohse 2000; Bhattacharyya 2006; Ahlers, Grossmann & Lohse 2009; Song, Shishkina & Zhu 2023):

(3.1)\begin{gather} \epsilon_{\nu} + \epsilon_{\eta} = \frac{\nu^3}{L^4} (Nu-1)Ra\,Pr^{{-}2}, \end{gather}

(3.2)\begin{gather}\epsilon_{\theta} = \kappa \frac{\varDelta^2}{L^2}Nu. \end{gather}

The convective part of $\epsilon _{\theta }$ can be expressed as

(3.3)\begin{equation} \tilde{\epsilon}_{\theta} = \epsilon_{\theta} - \epsilon_{cond} \implies Nu-1 \sim \frac{\tilde{\epsilon}_{\theta}}{\kappa \varDelta^2/L^2}. \end{equation}

Introducing $u, \theta \ \text {and} \ \ell$ as the representative scales for convective velocity, temperature and the length scale, the convective part of thermal dissipation rate can be expressed as $\tilde {\epsilon }_{\theta } \sim u\theta ^2/\ell$. Using this, after rearranging, we can rewrite (3.3) as (Song et al. 2023)

(3.4)\begin{equation} Nu-1 \sim \frac{\theta^2}{\varDelta^2} \frac{L}{\ell} Re\,Pr. \end{equation}

Nusselt number, $Nu$, can also be expressed in terms of convective heat transfer $q\sim u\theta$, as

(3.5)\begin{equation} Nu-1 \sim \frac{u \theta }{\kappa \varDelta/L}. \end{equation}

From (3.4) and (3.5), we get $\ell /L \sim \theta /\varDelta$, leading to

(3.6)\begin{equation} Nu-1 \sim \left ( \frac{\ell}{L} \right) Re\,Pr. \end{equation}

Assuming the induction $\boldsymbol {b}$, the total field can be written as $\boldsymbol {B} = \boldsymbol {B_0} + \boldsymbol {b}$. From Ohmic dissipation $\epsilon _{\eta } = (\rho \sigma )^{-1} \langle \overline {|\boldsymbol {J}|^2} \rangle$, and $\boldsymbol {J} = (\boldsymbol {\nabla } \times \boldsymbol {B})/\mu$, and assuming the length scale is $\ell$ we can write

(3.7)\begin{equation} \epsilon_{\eta} = \frac{1}{\rho \sigma \mu^2} \langle \overline{|\boldsymbol{\nabla} \times \boldsymbol{B}|^2} \rangle \sim \frac{1}{\rho \sigma \mu^2} \frac{b^2}{\ell^2}, \end{equation}

where $b$ represents the typical magnitude of the induced field. Under the quasi-static assumption, the induction equation can be simplified to

(3.8)\begin{equation} \eta \nabla^2 \boldsymbol{b} + B_0 \frac{\partial {u_x}}{\partial x} = 0, \end{equation}

which upon using the scale estimates for $\boldsymbol {b} \sim b$, ${u_x } \sim u$, $\partial ^2/\partial x^2 \sim 1/L^2$, $\partial ^2/\partial y^2 \sim 1/\ell ^2$, $1/\ell ^2\gg 1/L^2$ and $\partial /\partial x \sim 1/L$ becomes

(3.9)\begin{equation} \eta \frac{b}{\ell^2} \sim B_0 \frac{u}{L} \implies b \sim \frac{B_0}{\eta} \frac{u \ell^2}{L}. \end{equation}

Using this in (3.7), after rearranging and grouping relevant terms, we finally obtain a scaling law for the Ohmic dissipation with the horizontal length scale dependence,

(3.10)\begin{equation} \epsilon_{\eta} \sim \frac{\nu^3}{L^4} Re^2 Q \left ( \frac{\ell}{L} \right )^2. \end{equation}

With $Re$ and $\ell$ dependence of $Nu$ and $\epsilon _{\eta }$, we will proceed to derive the scaling relation of $Re \sim Re(Ra,Q,Pr,\ell /L)$ and $\ell /L \sim Q^{\gamma }$ in the columnar regime. It is to be noted the length scale of the structures in the columnar regime varies negligibly with $Ra$, which is why its $Ra$ dependency is not considered. In figure 7(a), in the columnar regime, we observe $\epsilon _{\eta }\gg \epsilon _{\nu }$. Hence from (3.1), we can consider $\epsilon _{\eta } \approx (\nu ^3/L^4) (Nu-1)Ra\,Pr^{-2}$. Combining this together with (3.10) and (3.6), we derive the scaling law for $Re$ in terms of the input parameters and the horizontal length scale,

(3.11)\begin{equation} Re \sim Ra Q^{{-}1} Pr^{{-}1}\left ( \frac{\ell}{L} \right )^{{-}1}. \end{equation}

Figure 7. (a) Fraction of viscous dissipation $f_{\nu } = \epsilon _{\nu }/(\epsilon _{\nu } + \epsilon _{\eta })$ (solid symbols) and Ohmic dissipation $f_{\eta } = \epsilon _{\eta }/(\epsilon _{\nu } + \epsilon _{\eta })$ (open symbols). (b) Dependence of the horizontal length scale on $Q$ in the columnar regime. The colour for both the open and solid symbols is the same as in figure 5.

Using this in (3.6), length scale and Prandtl number dependence drops out, yielding the Nusselt number scaling,

(3.12)\begin{equation} Nu -1 \sim \frac{Ra}{Q}. \end{equation}

This result is consistent with the scaling law obtained from the marginally stable thermal boundary layer analysis discussed above. For non-magnetic convection, a regime in which the entire fluid layer is turbulent (at asymptotically large values of $Ra$) is referred to as the ultimate regime. This regime, predicted by Kraichnan (1962), is characterized by the heat transport becoming independent of the diffusion coefficients. With additional constraining forces, like the Lorentz force in the present case, a heat transport regime independent of the diffusion coefficients $(\nu \text { and } \kappa )$ can be observed at relatively low Rayleigh numbers. It is noteworthy to point out that in the columnar regime ($Nu-1 \sim Ra/Q$), an inherently different flow state in comparison with the ultimate state of conventional RBC, the heat transport $Nu \times \kappa \varDelta /L$ in the limit of an asymptotically large $Q$, is independent of $\nu$ and $\kappa$. It does, however, depend on the magnetic diffusion coefficient $\eta$. It is also evident that the scaling of the length scale has no bearing upon the scaling law for the dimensionless convective heat transport in the columnar regime. Equation (3.12) is validated against our data in figure 8(a).

Figure 8. Scaling of the (a) Nusselt number, (b) Reynolds number and (c) Ohmic dissipation in the columnar regime.

One way to obtain the dependence of the horizontal length scale on the Chandrasekhar number is to obtain the scaling for the critical horizontal wavenumber $k_c$. From linear theory, it is well known that at large $Q$, $k_c \sim Q^{1/6}$, which implies $l_c \sim Q^{-1/6}$ (Chandrasekhar 1961). In figure 7(b), the $Q$-dependence of the length scale found in the data is compared with this result from the linear theory. We find a good collapse of data for all the simulation sets according to $\ell /L \sim Q^{-1/6}$ in the columnar regime. Using this in (3.10) and (3.11), we obtain the following relations for the Reynolds number and Ohmic dissipation scaling:

(3.13)\begin{gather} Re \sim RaQ^{{-}5/6} Pr^{{-}1}, \end{gather}

(3.14)\begin{gather}\epsilon_{\eta}\frac{L^4}{\nu^3} \sim Re^2 Q^{2/3} \sim Ra^2 Q^{{-}1} Pr^{{-}2}. \end{gather}

Data from our simulations show a good collapse for the above relations, see figure 8. In figure 8(a), at large $Q$, the $Nu$ data collapse in the columnar regime. Similarly, in figure 8(b), the $Re$ data at large $Q$ collapse, validating (3.13). It must, however, be noted for Reynolds number scaling, the exponent differs from one suggested by Yan et al. (2019) for 3-D quasi-static magnetoconvection. They derive the scaling relation for $Re$ by balancing the Lorentz force with the buoyancy in the columnar regime, which is evident from their 3-D data. Our derivations are based on energetic arguments involving the dissipation rates. Different assumptions invoked in these arguments is expected to lead to the differences observed in the exponent. The Ohmic dissipation data, shown in figure 8(c), at all $Q$ show a solid trend towards the scaling given by (3.14), collapsing for $Q$ in the asymptotic limit, $Q\rightarrow \infty$, in the columnar regime.

4. Concluding remarks

The effect of the horizontal length scale on the scaling behaviour of the response parameters in quasi-static magnetoconvection has not been explored in the literature. Previous studies concerning the energetic scaling arguments have considered the domain height $L$ as the sole length scale present in the system. In the present study, under the quasi-static assumption, we consider the horizontal width of the columns as a relevant length scale and explore the effect of its dependence on the Chandrasekhar number on the scaling of the dimensionless heat transport, flow velocity and Ohmic dissipation. Based on our DNS data, the derived scaling laws are validated successfully for Nusselt number, Reynolds number and Ohmic dissipation. The Nusselt number scaling obtained from our analysis is shown to be consistent with the marginal stability predictions, which further validates our theoretical analysis. The scaling law for the Reynolds number and Ohmic dissipation also shows a reasonable collapse in our data at large $Q$, implying the validity of our scaling arguments.

In the present work, we have considered periodicity in the horizontal direction and a fixed Prandtl number of unity. Compared with periodic boundaries, the presence of sidewalls has been shown to lead to convection at a much lower $Ra$ (Busse 2008). Numerical simulations of planar RBC by Liu, Krasnov & Schumacher (2018) confirmed the existence of wall-modes at $Ra \lesssim Ra_c$. Akhmedagaev et al. (2020) also confirmed the existence of these wall-modes at $Ra \gtrsim Ra_c$. The effect of sidewalls at even higher $Q$ and $Ra$ away from the onset will be something that we will explore in the future editions of this work. The presence of sidewalls may have implications in the design of liquid-metal batteries (Kelley & Weier 2018) and nuclear fusion technology (Ihli et al. 2008), wherein the geometries are often confined by walls.

Magnetoconvection is usually studied in liquid metals, which have very small Prandtl numbers and very high electrical conductivities. In the numerical and experimental low $Pr$ studies of Akhmedagaev et al. (2020) and Zürner et al. (2020), unlike the slender columns and magnetoplumes observed in our simulations, the dominant structures were observed to be a system of ascending and descending planar jets originating at the sidewalls. It must be noted, however, that the aforementioned studies consider convection in a cylindrical container in which sidewalls play a dominant role in organizing the flow. Hence, it remains to be seen if the slender columnar or plume structures can emerge in 2-D simulations by virtue of a low Prandtl number or exclusively in the presence of sidewalls. Whether or not the flow structures are different, we believe that it does not affect the scaling argument for the conditions that we have considered. That can be seen from the low $Pr$ ($Pr=0.025$) simulations of Yan et al. (2019), who performed the comparative calculations at $Pr = 1$ and $Pr = 0.025$. In their work, it was observed that the scaling behaviour of the dimensionless heat transport seems to follow the same trend as in the case of $Pr$ being unity, i.e. $Nu \sim Ra/Q$ at asymptotically large values of $Q$; since in the columnar regime, there is no $Pr$ dependence of $Nu$. Although, $Re$ and $\epsilon _{\eta }$ show $Pr$ dependence in this regime, in this work, our main goal was to test the scaling argument that we proposed for asymptotically large $Q$. The $Re$ and $\epsilon _{\eta }$ scaling relations show a good agreement with what we observed in our DNS data (for $Q$ and $Ra$ dependence). Exploration of the $Pr$ dependence of the response parameters and the flow structures is a full-fledged campaign in itself that we propose to study in future.