Unifying heat transport model for the transition between buoyancy-dominated and Lorentz-force-dominated regimes in quasistatic magnetoconvection

Abstract In magnetoconvection, the flow of an electromagnetically conductive fluid is driven by a combination of buoyancy forces, which create the fluid motion due to thermal expansion and contraction, and Lorentz forces, which distort the convective flow structure in the presence of a magnetic field. The differences in the global flow structures in the buoyancy-dominated and Lorentz-force-dominated regimes lead to different heat transport properties in these regimes, reflected in distinct dimensionless scaling relations of the global heat flux (Nusselt number $Nu$) versus the strength of buoyancy (Rayleigh number $Ra$) and electromagnetic forces (Hartmann number $Ha$). Here, we propose a theoretical model for the transition between these two regimes for the case of a static vertical magnetic field applied across a convective fluid layer confined between two isothermal, a lower warmer and an upper colder, horizontal surfaces. The model suggests that the scaling exponents $\gamma$ in the buoyancy-dominated regime, $Nu\sim Ra ^\gamma$, and $\xi$ in the Lorentz-force-dominated regime, $Nu\sim (Ha^{-2}Ra)^\xi$, are related as $\xi =\gamma /(1-2\gamma )$, and the onset of the transition scales with $Ha^{-1/\gamma }Ra$. These theoretical results are supported by our direct numerical simulations for $10\leq Ha\leq 2000$, Prandtl number $Pr=0.025$ and $Ra$ up to $10^9$ and data from the literature.


Introduction
Magnetoconvection (MC) governs most astro-and geophysical systems and is relevant to various engineering applications [25,7].The former include, for instance, outer layers of stars and liquid metal planetary cores [12], examples of the latter comprise liquid-metal batteries, induction heating, casting, liquid-metal cooling for nuclear fusion reactors and semiconductor crystal growth [6].MC occurs in an electrically conducting fluid that is subjected both to a magnetic field and an imposed temperature gradient.The buoyancy forces induce convective fluid motion due to thermal expansion and contraction, while the magnetic field affects this motion and distorts the global flow structure through the Lorentz force, which eventually influences the heat transport in the system.The resulting main two control parameters, the strength of the imposed thermal driving and that of the external magnetic field, are encoded in independent dimensionless groups, the Rayleigh number Ra and Hartmann number Ha, respectively.
One of the key objectives in MC research is to provide scaling relations for the heat transport through the system, represented in dimensionless form by the Nusselt number Nu, as a function of Ra and Ha.However, the heat transport scaling relations also depend on the flow configuration, including the angle between the magnetic field and gravity, the geometry of the container and the boundary conditions (BCs), and on whether the buoyancy forces dominate over the Lorentz forces in the system or vice versa.This inherent complexity results in the need, at least in principle, to derive separate heat transport scaling relations to describe each specific flow regime itself and transitions between distinct regimes.The considerable difficulty of doing so in a coherent manner is exacerbated by non-universal scaling relations even within specific regimes -the scaling relations in the buoyancy-dominated and Lorentz-force-dominated regimes themselves change with the control parameters, and transitions between the different regimes are also non-universal.
The objective of this paper is to offer a unifying heat transport model for the transition between the buoyancy-dominated and Lorentz-force-dominated regimes in quasistatic MC.We focus on Rayleigh-Bénard convection (RBC) [2] with an applied vertical magnetic field and assume that the magnetic field is constant in the entire domain, without being affected by a fluid motion or finite magnetic diffusion.The model uses the theoretical predictions by Grossmann and Lohse [9,10,24] for RBC without magnetic field and transfers the approach by Ecke and Shishkina [8, §3.3] for transitions in rotating RBC to the case of RBC with a vertical magnetic field.To verify the proposed model, we compare the theoretical predictions with results for liquid metal MC obtained by direct numerical simulation (DNS) carried out by us and others [16,27,3,26], as well as experiments [5,13,28,26].In addition, we carried out simulations for a different working fluid at higher Prandtl number Pr to compare our DNS data with that from Ref. [15].The predictions of the proposed model agree well the experimental and DNS data.

Model for the transition between the buoyancy dominated and Lorentz force dominated regimes
We consider a layer of electrically conducting fluid confined between two infinitely wide and long plates, driven by a buoyancy force generated by an imposed vertical temperature difference between the top and bottom plates, and subjected to a uniform vertically orientated magnetic field.When the magnetic field is weak compared to the buoyancy force, we recover classical RBC scaling for the dimensionless convective heat flux Nu − 1, that is, the total dimensionless heat flux Nu less its conductive contribution, Nu − 1 ∼ (Ra/Ra c,b ) γ ∼ Ra γ , for an exponent γ, where where Ra c,b is the critical Ra for RBC bulk onset for a given container geometry.When the Lorentz force is strong compared to the buoyancy force, we expect a similar scaling law Nu − 1 ∼ (Ra/Ra c,L ) ξ , for an exponent ξ.Here, the dependence on the critical Rayleigh number Ra c,L in the Lorentz-force-dominated regime is kept due to its dependence on Ha, that can be obtained from the linear stability theory [4] Ra c,L ∼ Ha 2 , with no dependence on the Prandtl number.
Although these buoyancy-dominated and Lorentz-force-dominated scaling laws appear disconnected, they are intrinsically linked under the assumption that they must overlap at some intermediate region between the two extreme regimes, where neither the influence of the Lorentz force or buoyancy force on the convective heat transport can be ignored.At the transition between these regimes, the two corresponding scaling laws will scale in the same way, To construct a relationship between the two scaling exponents γ and ξ, we assume that this cross-over occurs when the thicknesses of the thermal and viscous boundary layers (BLs) scale in the same way, motivated by the observation that the viscous layer is nested within the thermal BL in the Lorentz-force-dominated regime and vice versa in the buoyancy-dominated regime.In a domain of height H over a semi-infinite plate, the average thermal BL thickness from laminar Prandtl-Blasius boundary layer theory [18,21] is δ θ = H/(2Nu) ∼ Nu −1 , while the laminar viscous BL influenced by a vertical magnetic field is the Hartmann layer [11,7], δ Ha = cH/Ha ∼ Ha −1 , where c is a constant.Assuming that δ θ ∼ δ Ha at the transition implies Nu ∼ Ha.Since this transition is seen to typically occur at high Nusselt numbers, meaning Nu ≈ Nu − 1, we obtain resulting in the following relationship between the exponents, One can see that a larger (smaller) exponent in one regime requires a larger (smaller) exponent the other regime, and that γ is always smaller than 1/2.In Fig. 1 we present a sketch of the proposed scaling relations for (Nu − 1) versus Ha (Fig. 1 a) and Ra (Fig. 1 b), according to the relations (3).Once γ is known for any specific Pr, the exponent ξ can be calculated from Eq. ( 3).These scalings can then be used define coordinates (Nu − 1)Ra −γ and Ha −1/γ Ra with respect to which the heat transport dependence for different values of Ha and Ra collapse onto a master curve, as sketched in Fig. 2. The transition then should take place in a Rayleigh-number range that scales as Ha 1/γ .Figure 2: Schematic representation of the normalized convective heat transport Nu − 1 displaying the transition from the Lorentz-force-dominated regime, Nu − 1 ∼ (Ha −2 Ra) ξ , to the buoyancy-dominated regime, Nu − 1 ∼ Ra γ , according to our model.The scaling exponents ξ and γ follow Eq. ( 3), while the transition scales with Ha −1/γ Ra.
To close the model, a theoretical prediction for either γ or ξ must be made.The former is readily available through Grossmann-Lohse (GL) theory [9,10,24] for RBC without magnetic field, applicable here in the buoyancy-dominated regime.For any given Pr and Ra-range, the theory provides accurate predictions of the value of γ, for containers of aspect ratio Γ ≳ 1.For Γ ≪ 1, the data can be rescaled according to the method suggested in [22,1], which we do not discuss here, as in the present study Γ = 1.
3 Status Quo -experimental and numerical data.
To verify the theoretical model we compare its predictions against data obtained from experiments of liquidmetal MC [5,13,28,26] and DNS conducted by us and others [16,27,3,26].However, before doing so we provide a brief overview of the data collated from the literature and produced by us to demonstrate the considerable challenges that arise when trying to draw firm conclusions on the scaling of the heat transport with Ha and Ra.
We simulate an incompressible, viscous buoyancy-driven flow of an electrically conducting fluid in the presence of an imposed magnetic field for very small magnetic Reynolds number Re m ≪ 1 and magnetic Prandtl number Pr m ≪ 1 by numerically solving the MC equations within the Oberbeck-Boussinessq and quasistatic approximation, where u is the velocity, T the temperature, p the kinematic pressure, j the electric current density, ϕ the electric potential, and e z and e B are unit vectors that point, respectively, upward (opposite to gravity) and in the direction of the magnetic field B = B 0 e B .The magnetic field is aligned with the buoyancy force, e B = e z .Equations ( 4)-( 8) have been non-dimensionalized using the container height H, the free-fall velocity u ff ≡ (αgH∆) 1/2 , the free-fall time t ff ≡ H/u ff , the temperature difference between the bottom and top plates, ∆ ≡ T + − T − , and the external magnetic field strength, B 0 , as units of length, velocity, time, temperature and magnetic field strength, respectively.The dimensionless control parameters are the Rayleigh number Ra, the Prandtl number Pr, and the Hartmann number Ha, where σ is the electrical conductivity, ρ the mass density, α the thermal expansion coefficient, g the acceleration due to gravity, ν the kinematic viscosity, and κ the thermal diffusivity.We apply no-slip BCs for the velocity at all boundaries, u = 0, constant temperatures at the end-faces, i.e., T = T + at the bottom plate at z = 0 and T = T − at the top plate at z = H, and adiabatic BC at the side walls, ∂T /∂n = 0, where n is the vector orthogonal to the surface.All solid boundaries are considered electrically insulated; Neumann BCs for the electric potential are ∂ϕ/∂n = 0.The simulation domain is cubic of height H, width W , length L, H = W = L, i.e. has aspect ratio Γ ≡ L/H = 1.Most simulations are carried out for liquid metals such as GaInSn at Pr = 0.025, for Ra up to 10 9 and Ha up to 2000.Some DNS are conducted also for Pr = 8, to extend the parameter range studied in Ref. [15].In the bulk, spatial flow fluctuations are resolved down to 2-5, occasionally 10 Kolmogorov miscroscales, near the rigid walls we resolve the thermal and Hartmann BLs [23].Our DNSs have been carried out with an MC extension of goldfish [14,20,19].Further details of the simulations are provided in appendix A.
In Fig. 3 we present Nu − 1 as a function of Ra (Fig. 3 a-b) and Ha (Fig. 3 c-d), respectively, from our DNS, experimental [5,13,28,26] and DNS data [16,3,26] for liquid metals, 0.025 ≤ Pr ≤ 0.029 (Fig. 3 b, d) and a fluid with Pr = 8 (Fig. 3 a, c).In addition, in Fig. 3 (b, d), we plot for comparison the DNS data [27] for free-slip BCs.In Fig. 3 (a-b) one can see that Nu generally increases with growing Ra, but with different slopes for different Ha, which are steeper for larger Ha.In the double logarithmic plots of Fig. 3 (a-b), the curves of the (Nu − 1)-vs.-Radependences for different Ha approach each other when Ra increases.In Fig. 3 (c-d) one can see that Nu remains to be almost unaffected by the magnetic field for relatively small Ha, but for a strong Lorentz force (large Ha) they gradually decrease with growing Ha.Here, again, the decreasing slopes are different for different Ra, and the transition to the regime, where the heat transport is affected by the magnetic field, depends on Ra.In summary, the data in Fig. 3 look rather different in

Model validation
We now validate the model using the data presented in Fig. 3. To calculate the scaling exponent γ in the buoyancy-dominated regime, for the considered Pr and Ra-ranges, we use GL theory, which gives γ about 0.30 for Pr = 8 and about 0.31 for Pr = 0.025.These values agree very well with fits to data for the Ra-ranges in the buoyancy-dominated regime for both values of Pr, see appendix A for further details.Using Eq. ( 3) with γ = 0.30 for Pr = 8 and γ = 0.31 for Pr = 0.025, we calculate the exponent ξ in the Lorentz-force-dominated regime, which equals ξ = 0.75 for Pr = 8 and ξ ≈ 0.82 for Pr = 0.025.In Fig. 4 we plot all data presented previously Fig. 3 using the coordinates suggested by our model and visualised in Fig. 2.This results in a clear collapse of the data onto master curves for Pr = 8 (Fig. 4 a) and 0.025 ≤ Pr ≤ 0.029 (Fig. 4 b).Some deviation of the data from the master curve in Fig. 4 (b) is observed when Ra is relatively small and the flow is in the wall-mode regime, that is, before the onset of bulk convection where heat transport is confined to the near-wall region.Since the influence of no-slip side-walls is outside of the scope of our model, the observed deviations are expected in this regime.3).Pink and blue lines show the predictions of the slopes in the buoyancy and Lorentz-force-dominated regimes, respectively.The symbols have the same meaning as in Fig. 3.

Conclusion
In conclusion, we have proposed a heat transport model for the transition between the buoyancy-and Lorentzforce-dominated regimes of vertical MC.We validated the model using our DNS and data available in the literature.We wish to emphasize that the proposed model is parameter-free.For a given Pr and Ra-range, one can calculate the scaling exponent in the buoyancy-dominated regime, using GL theory.Then, using Eq. ( 3), one can calculate the scaling exponent ξ in the Lorentz-force-dominated regime and collapse the data on a master curve by rescaling the coordinate axes as in Fig. 2. The model can in principle be extended to include the effect of a fluctuating magnetic field, this merely results in an adjustment of the prefactors.

A Appendix
A.1 Numerical simulations All direct numerical simulations have been carried using the direct numerical solver goldfish [14], which has been widely used in previous studies of different convective flows.The new version of the code that applies a fourth-order finite-volume discretisation on staggered grids and a third-order Runge-Kutta time marching scheme [20,19] has been extended to simulate magnetoconvective flows, where a consistent and conservative scheme [17] is utilised to calculate the current density and the Lorentz force.The three-dimensional DNS are performed using staggered grids that provide fine resolution in the core part of the domain and near the rigid walls [23], to resolve the thermal and Hartmann boundary layers.The main parameters and key observables of the simulations are summarised in table 1.
The DNS dataset comprises 37 simulations to cover the necessary ranges in Ha and Ra.To obtain a dataset this large within the available resources, a compromise had to be made in terms of the bulk resolution.Hence a number of simulations do not resolve scales close to the Kolmogorov microscale.To ensure accurate predictions of mean Nu, grid refinement studies have been done for key simulations (marked by * in the h DNS /h K column).For these cases, changing the grid resolution by a factor of at least 2.5 resulted in changes of Nu by less than 1% over a long-time average.Figure 5 shows a comparison between the predictions of Grossmann-Lohse theory [9,10,24] for the dependence of the dimensionless heat transport Nu − 1 on the Rayleigh number Ra, Nu − 1 = A Ra γ for Prdependent prefactor A. Fits to data have been carried out for 3 × 10 6 ⩽ 10 9 for Pr = 0.025 and 3 × 10 6 ⩽ 10 9 for Pr = 8.As can been seen from the data presented in fig.5, the theoretical predictions agree very well with the data.

Figure 1 :
Figure 1: Scalings of Nu − 1 versus (a) Ha and (b) Ra, according to the theory.

Figure 4 :
Figure 4: All data from Fig. 3 follow master scaling curves if plotted as Fig. 2 suggests, for (a) Pr = 8 and (b) 0.025 ≤ Pr ≤ 0.029.The values of γ are calculated from GL theory, and the values of ξ are calculated from Eq. (3).Pink and blue lines show the predictions of the slopes in the buoyancy and Lorentz-force-dominated regimes, respectively.The symbols have the same meaning as in Fig.3.

Figure 5 :
Figure 5: Predictions of GL theory [9, 10, 24] (thin lines) for Nu − 1 versus Ra scaling relations in RBC without magnetic field, for Pr = 8 (blue) and Pr = 0.025 (pink).The power fits for the considered Ra-ranges are highlighted with thick lines of the corresponding colors.

Figures 6
Figures 6 and 7 present visualisations the velocity magnitude and temperature, respectively, at an instant in time during statistically steady evolution in the magnetically dominated and the buoyancy-dominated

Table 1 :
Details on the conducted DNS, including the standard deviation σ Nu of the Nusselt number Nu, the number of nodes N x , N y , N z in the directions x, y and z, respectively; the number of simulation freefall times used for averaging T run ; the number of the nodes within the thermal boundary layer, and within the Hartmann boundary layer N Ha ; the Kolmogorov microscale, h K , and the relative mean grid stepping, h DNS /h K .Simulations marked by * are those for which grid refinement studies have been completed.