A Hartmann layer

# Heat and momentum transfer for magnetoconvection in a vertical external magnetic field

## Abstract

The scaling theory of Grossmann and Lohse (J. Fluid Mech. 407, 27 (2000)) for the turbulent heat and momentum transfer is extended to the magnetoconvection case in the presence of a (strong) vertical magnetic field. The comparison with existing laboratory experiments and direct numerical simulations in the quasistatic limit allows to restrict the parameter space to very low Prandtl and magnetic Prandtl numbers and thus to reduce the number of unknown parameters in the model. Also included is the Chandrasekhar limit for which the outer magnetic induction field is large enough such that convective motion is suppressed and heat is transported by diffusion. Our theory identifies four distinct regimes of magnetoconvection which are distinguished by the strength of the outer magnetic field and the level of turbulence in the flow, respectively.

Rayleigh-Bénard convection; turbulent convection; magnetoconvection
###### pacs:
47.20.Bp, 47.27.te, 47.65.Cb, 44.25.+f

## I Introduction

One of the central questions in turbulent convection is that of the global transport of heat and momentum as a function of the thermal driving and the properties of the working fluid (1); (2); (3). In the simplest setting of turbulent convection – the Rayleigh-Bénard case – one considers an infinitely extended horizontal layer of fluid which is uniformly heated from below and cooled from above. The thermal driving of the turbulent convection in the layer is then established by the temperature difference between the top and the bottom, , and directly proportional to the dimensionless Rayleigh number . The properties of the working fluid are determined by the Prandtl number , defined as the ratio of the kinematic viscosity to the thermal diffusivity . Turbulent heat and momentum transfer are quantified by the dimensionless Nusselt, , and Reynolds, , numbers, respectively. In a nutshell, one seeks for and being functions of and .

One of the oldest scaling theories that aimed at predicting at fixed dates back to Malkus (4); (5) and is based on a marginal stability argument for the turbulent mean profiles. More recently, scaling theories by Shraiman and Siggia (6); (7) and Grossmann and Lohse (GL) (8); (9) have been developed. The central idea of the GL theory is a decomposition of the thermal and kinetic energy dissipation into contributions from the bulk and the boundary layers in the vicinity of the plates. These contributions have to be weighted with the volume fractions that the boundary layers of the temperature and velocity fields occupy. The theory is adapted to doubly diffusive convection (10) and horizontal convection (11).

In astrophysical systems, thermal convection is often tightly coupled to magnetic fields (and rotation) which is known as magnetoconvection (12). Examples are sunspots in the solar chromosphere (13) or the X-ray flaring activity of some young neutron stars which are termed magnetars (14). Less spectacular, but not less important are numerous industrial applications reaching from materials processing, such as crystal growth by the Czochralski method (15) or dendritic solidification in alloys (16), to fusion technology (17). In case of a strong prescribed magnetic induction the secondary magnetic induction , which is generated by flow-induced eddy currents, remains very small. While a strong vertical external field can then damp and even suppress the convective fluid motion (18), convection rolls can be stabilized when the magnetic field is applied in horizontal direction (19). From a standard magnetohydrodynamic (MHD) perspective, the turbulence of coupled velocity and magnetic fields is then constrained. This regime is known as the quasistatic regime of MHD: the Lorentz force enters the momentum equation, the induction equation which describes the temporal evolution of the magnetic induction field is however neglected (20).

The aim of the present work is to extend the GL theory of turbulent transport to the case of magnetoconvection. First attempts in this direction have been reported by Chakraborty (21). He showed that an Ohmic dissipation rate, , has to be incorporated beside the thermal and kinetic energy dissipation rates, and . One is thus left with eight different regimes of boundary-layer- and/or bulk-dominated dissipation rates. Together with free parameters for the viscous boundary layer thickness and a critical Reynolds number for the crossover from low to high Prandtl numbers (9), one ends up with at least ten parameters to fit. Furthermore, dimensionless parameters have to be added that relate the electrical conductivity either to the kinematic viscosity or the thermal diffusivity and quantify the strength of the outer magnetic field. In view to this significant extension of the parameter space, one has to seek for regimes of magnetoconvection that can be studied with a reduced set of fit parameters.

We will therefore restrict the turbulent magnetoconvection to a specific parameter range. In view to a comparison with laboratory experiments of magnetoconvection, which are typically conducted in liquid metals, one can restrict the Prandtl number range to

 Pr=νκ≲10−2. (1)

Also the range of the magnetic Prandtl number can be limited to

 Pm=νη=RmRe≲10−5. (2)

with the diffusivity of the magnetic induction and being the permeability. In many laboratory flows the magnetic Reynolds number will thus remain small, . This regime is termed the quasistatic case of magnetohydrodynamics. The magnetic field lines cannot be bended significantly by the fluid motion since the magnetic diffusion time scale is very short. This excludes some astrophysical applications such as interstellar turbulent gases in which (22).

Similar to standard GL theory, our predictions have to be fitted to one reference data set. Our adjustment of the free coefficients will be based on an experiment by Cioni et al. (23) which is to the best of our knowledge the only experiment that was operated at a sufficiently high Rayleigh number. Further data records by Burr and Müller (24) and Aurnou and Olson (25) have been conducted at smaller Rayleigh numbers and will be discussed only briefly. In addition, our own direct numerical simulations of magnetoconvection in the quasistatic regime will be included to obtain (at least one) data point with known Reynolds and Nusselt numbers at given Rayleigh, Hartmann (the dimensionless measure for magnetic field strength which will be defined in section II) and Prandtl numbers.

The outline of the work is as follows. In the next section, the set of magnetoconvective equations of motion is discussed, the characteristic scales, dimensionless parameters and dissipation rates are defined. Also the numerical method and a short description of the data sets will be presented. This section is followed by a derivation of the nonlinear equations for and . Finally the free parameters of the scaling theory are fitted to data records. The results are summarized and discussed in brief at the end of the work.

## Ii Equations and parameters

### ii.1 Quasistatic equations of magnetoconvection in Boussinesq approximation

We solve the three-dimensional Boussinesq equations for turbulent magnetoconvection in a rectangular cell of height and side lengths in the quasistatic limit. The equations for the velocity field and the temperature field are given by

 ∇⋅u =0, (3) ∂u∂t+(u⋅∇)u =−1ρ0∇p+ν∇2u =+gα(T−T0)ez+1ρ0(j×B0), (4) ∂T∂t+(u⋅∇)T =κ∇2T. (5)

The pressure field is denoted , is a reference temperature, the constant mass density and the magnetic field. The Ohm law for the current density is given by

 j =σ(−∇ϕ+u×B0), (6)

where the electric potential follows from . The Rayleigh number is given by

 Ra=gαΔTH3νκ, (7)

and the Hartmann number by

 Ha=B0H√σρ0ν=√Q. (8)

The square of is also known as the Chandrasekhar number . The variables , and denote the acceleration due to gravity, the electrical conductivity and the thermal expansion coefficient, respectively. In a dimensionless form length scales are expressed in units of , velocities in units of the free-fall velocity , temperature in units of the outer difference and magnetic induction in units of . The configuration is sketched in figure 1.

### ii.2 Direct numerical simulations

The equations (3) – (6) are solved for a closed Cartesian cell with a second-order finite difference scheme. The projection-type scheme is nearly fully conservative. The advection-diffusion equation is solved by semi-implicit scheme in which nonlinear terms are treated explicitly and diffusion terms implicitly. The program applies MPI and Open MP. More details are found in (26). For the fit, we will use two series of direct numerical simulations (DNS):

• Series 1: , , . The aspect ratios are and . The grid is non-uniform and contains points.

• Series 2: , , . The aspect ratios are and . The grid is non-uniform and contains points.

The boundary conditions are as follows: all walls are electrically insulated walls, i.e. the field lines of the current density are closed inside the fluid volume. No-slip boundary conditions hold for the velocity at all walls, the top and bottom walls are additionally isothermal with prescribed temperatures and , respectively. The side walls are thermally insulated. The grid is clustered at the top and bottom walls to resolve the Hartmann layers at the top and bottom and first order quantities. We have also performed grid-sensitivity studies to make sure that the Nusselt number remains constant plane-by-plane (plane at constant height ).

### ii.3 Dissipation rate balances

In correspondence with classical Rayleigh-Bénard convection, we can derive exact relations for the mean kinetic energy dissipation rate, , the mean magnetic dissipation rate, , and the mean thermal dissipation rate, . The fields are defined as

 ε(x,t) =ν2(∂iuj+∂jui)2, (9) εB(x,t) =η2μρ0(∂ibj−∂jbi)2, (10) εT(x,t) =κ(∂iT)2, (11)

with . Since is constant equation (10) contains derivatives of the induced magnetic induction only which arise from the eddy currents . In the statistically stationary regime we obtain

 ε+εB =ν3H4(Nu−1)RaPr2, (12) εT =κ(ΔT)2H2Nu. (13)

The Nusselt number, which quantifies the turbulent heat transfer, is given by

 Nu=1+H⟨uzT⟩κΔT. (14)

The global momentum transfer in the magnetoconvective system is quantified by the Reynolds number which is defined as

 Re=⟨u2i⟩1/2Hν. (15)

In both definitions stands for volume-time average or ensemble average. While the thermal balance remains unchanged in comparison to the classical Rayleigh-Bénard case, the kinetic energy balance differs by the addition of on the left hand side of Eq. (12). It results from the Joule dissipation in the presence of a magnetic field. For completeness, we also list the definition of the magnetic Reynolds number

 Rm=UHη=μσUH, (16)

where is again given by the root mean square velocity, .

## Iii Extension of the scaling theory of Grossmann and Lohse

The central idea of the scaling theory is to combine Eqns. (12) and (13) with a decomposition of dissipation rates into contributions coming from the bulk and the boundary layers (BL) (8); (9). The following modifications are made to predict and for our case at hand:

1. The relevant boundary layer for the velocity field is the Hartmann layer (20) (see also appendix),

 δv=HHa, (17)

while the thermal boundary layer thickness remains . Contrary to the original GL theory we do not have the free parameter that appears in the Prandtl-Blasius-type expression .

2. We limit the study to low Prandtl numbers as already mentioned in the introduction. Thus the modification for the limit of large Prandtl numbers which has been developed in (9) and the related parameter are not necessary here. This saves a second fit parameter.

3. It is well-known from the linear stability analysis (18) that the critical Rayleigh number scales as

 Rac=π2Ha2. (18)

If is too big at a given , convection is suppressed completely.

The mean energy dissipation rates will be composed of a boundary layer contribution and a bulk contribution. This results to

 ε =εbulk+εBL, (19) εB =εB,bulk+εB,BL, (20) εT =εT,bulk+εT,BL. (21)

The dimensional estimates of the different contributions are given by

 εbulk ∼U3H =ν3H4Re3, (22) εBL ∼νU2δ2vδvH =ν3H4Re2Ha, (23) εB,bulk ∼ημρ0Rm2B20H2 =ν3H4Re2Ha2, (24) εB,BL ∼ημρ0Rm2B20δ2vδvH =ν3H4Re2Ha3, (25) εT,bulk ∼(ΔT)2UH =κ(ΔT)2H2RePr, (26) εT,BL ∼κ(ΔT)2H2√RePr. (27)

The bulk scalings of the kinetic and thermal dissipation rates in (22) and (26) are the same as in the original GL theory (9). The argumentation in (9) that leads to (27) remains valid for the present case. However, the scaling relation in (23) differs to the original case. Instead of the original BL expression , we insert the Hartmann layer thickness (17). For the new estimates in (24) and (25) we use the definition of which is given in (10) and measure the induced magnetic field in units of .

Following Grossmann and Lohse (9), we introduce interpolation functions to account for changes of the scaling laws in different parameter regimes. Once becomes smaller than the dominant velocity in the thermal BL changes from to . This is accounted for by replacing with in (26) and (27), where

 f(xT)=1(1+xnT)1/n (28)

with the argument and . For this interpolation function follows that and .

Close to the critical Rayleigh number the bulk of the fluid becomes laminar and scales with rather than as in (22) for the turbulent regime. This change is modelled by multiplying (22) with

 g(x∗) =1f(1/x∗) (29)

with the argument . From the definition of follows that and . The Reynolds number marks the range in which the transition from fully developed turbulence to weakly nonlinear time-dependent regime of velocity dynamics takes place. Combining all pure scaling laws with the interpolations as just described gives

 (Nu−1)RaPr2Re2 =c1Reg(ReRe∗) =+c2Ha+c3Ha2+c4Ha3, (30)
 Nu−1 =c5RePrf(2NuHa)+c6√RePrf(2NuHa). (31)

with the seven a priori unknown model parameters and to which have to be determined from a data record. The set of implicit equations can then be solved to obtain expressions and . While it is not possible to find a full solution analytically, can be calculated from (31) as a function of , , and :

 Re =(√c26+4c5(Nu−1)−c6)24c25Prf(2NuHa) (32)

Inserting (32) into (III) gives an equation independent of . However this new equation cannot be solved analytically for .

The stabilizing effect (iii) of large is included here in the following way: assuming we have found an analytical expression we can enforce the transition to the non-convective regime at by multiplying with

 h(xc) =1−f(xc), (33)

where . The function obeys the properties and , which ensures that in the purely diffusive equilibrium. The crossover function transitions smoothly between these two states, so that at we have instead of an abrupt jump to zero. Since we cannot determine directly we transform into and in the -independent equation we replace by . This gives the same result of in the non-convective regime once the equation is solved for by numerical methods. Thus the final model equations are (32) and

 (Nu−1)Raζ2Pr2h(Ra/Rac) =c1ζg(ζRe∗) =+c2Ha+c3Ha2+c4Ha3, (34)

with

 ζ=(√c26+4c5(Nu−1)h(Ra/Rac)−c6)24c25Prf(2NuHa). (35)

Now (III) can be used to determine the seven model parameters and to by fitting the equation to a data set . However examining (III) shows, that it is invariant for the following transformations:

 c1 →α6c1, c2 →α4c2, c3 →α4c3, c4 →α4c4, c5 →α2c5, c6 →αc6, Re∗ →Re∗/α2

for any . This means that the optimal values for the model parameters are ambiguous. To fix this ambiguity we need at least one data point which includes the Reynolds number. Then (32) can be used to calculate as a function of :

 c6 =Nu−1√RePrf(2NuHa)−c5√RePrf(2NuHa). (36)

With this step the optimal values of all remaining six model parameters and to are unique. It is absolutely clear that six parameters, which have to be adjusted, is still a large number. Nevertheless, one has to keep in mind that the number of free parameters has already been reduced significantly. We are not aware of any publications that report magnetoconvection data sets including . Therefore, we are using our own numerical simulations to determine data points for evaluating (36).

## Iv Results

Our numerical simulations are used to evaluate (36). After substituting (36) into (III), the resulting equation is fitted to the experimental data of Cioni et al. (23) in terms of and to , utilizing the Levenberg-Marquardt method (27). The experimental data have been obtained for convection in liquid mercury at a Prandtl number of . Our DNS are conducted at the same Prandtl number. With the known optimal model parameters we can calculate by solving (III) numerically for given , and and subsequently obtain from (32). The optimal model parameters are , , , , , and From (36) we get . The --phase diagrams for and calculated with these parameter values for are shown in figure 2. The top panel of the figure shows the magnitude of the Nusselt number as a function of and . The bottom figure displays the Reynolds number depending on both parameters. Also added are the experimental and DNS data. In figure 2, we also display the Chandrasekhar limit above which and .

Furthermore, the line is displayed for which . Above this line the Hartmann layer thickness will be smaller as the thermal boundary layer thickness. This characteristic line is crossed by a second line that shows . As mentioned already in section III (see equation (18)), on the left side of this line the convection flow is not fully developed turbulent, but in weakly nonlinear and time-dependent convection state. All data which are to the right of this line can be considered as fully turbulent convection data. It can be seen that only a few data points of (23) cross this threshold. The parameter space, thus, splits into four subregions by both lines:

• Region I: weakly nonlinear flow and strong magnetic field

• Region II: fully developed turbulent flow and strong magnetic field

• Region III: fully developed turbulent flow and weaker magnetic field

• Region IV: weakly nonlinear flow and weaker magnetic field

A few words about the quality of the fit should be addressed now. First, we mention that the size of the error bars of all fit coefficients (except ) is of the order of 100 %. In case of the coefficient this error level is even exceeded (see also next paragraphs). This is caused by the sparse record of data points. As can be seen in the figure, the data of Cioni et al. (23) are collected for three different Hartmann numbers that cover a small range. Also, these data reach only to the beginning of regime II. Regimes III and IV do not contain any data points. Stevens et al. (28) demonstrated in their recent update of GL theory that the uncertainties in the coefficients can be significantly reduced when the data cover a wide range of parameters. Furthermore, these three Hartmann numbers are much larger than those from our DNS. The additional data by Burr and Müller (24) or by Aurnou and Olson (25) have been conducted close to the onset regime of convection. Their experimental data are thus rather in the weakly nonlinear than in the fully turbulent range and will not be used for our study.

Secondly, it is observed that two fit coefficients, and , are negative. Although and thus practically zero, the corresponding term in (III) can give a non-negligible contribution to the scaling due to . Coefficient with the biggest error bar needs further consideration. Figure 3 displays the five coefficients in dependence of a fixed . To get this figure, we repeated the fits at each fixed value of the crossover Reynolds number. It is seen that the results for to are nearly insensitive for . Beyond this value, coefficient changes sign which is indicated by a dashed line in the plot. The eventual value of falls into a range, where small variations of cause large changes of (including sign changes).

The magnitude of in our fit corresponds to a Rayleigh number of . This estimate follows from recent numerical studies in liquid metal convection without magnetic field (30). It falls thus consistently into the range, for which convection develops into the fully developed turbulent regime which is also known as the hard convective turbulence regime (29). At the moment, we can only speculate that the inclusion of more data could lower the value of as it is expected in low- convection (see e.g. (31); (32)).

Thirdly, in order to quantify the impact of the error bars of the fit coefficients on and , we proceeded as follows. The six coefficients and were chosen randomly and statistically independently within their error bars. With these 6-tuples the parameter dependence and is reconstructed for 118 different cases. The superposition of these individual realizations results in an relative error around the original value in figure 2. The magnitudes of the relative error of both, Nusselt and Reynolds number, are plotted in logarithmic units in figure 4. The relative error of is highest along the border between regime I and II, but does not exceed 40 %. On the other hand the relative uncertainty of rises for smaller and reaches more than 100 % for below .

## V Summary

We have presented an extension of the scaling theory of Grossmann and Lohse (8); (9) to a convection layer in the presence of a vertical magnetic field. The discussion is restricted to magnetoconvection at low Prandtl and magnetic Prandtl numbers. In this regime the quasistatic approximation is applied that allows a significant reduction of the number of free parameters in the flow at hand and thus an application of the ideas of GL theory. Below the Chandrasekhar limit four different convection regimes are identified. On the one hand, they follow from the ratio of the Hartmann and thermal boundary layer thicknesses. On the other hand, the regions result from the critical Reynolds number , beyond which the convection flow is assumed to be fully turbulent.

In contrast to standard Rayleigh-Bénard convection, the data base is very small. In fact, there is only one data set from Cioni and co-workers, that can be used to fit the free parameters. The remaining data (24); (25) fall into a completely different section of the parameter plane. In particular, they remain close to the Chandrasekhar limit and cannot be used for turbulent magnetoconvection. This limits the predictive capabilities of our scaling results and calls for additional experimental data which are planed in the near future.

###### Acknowledgements.
TZ and WL are supported by the Research Training Group GK 1567 on Lorentz Force Velocimetry which is funded by the Deutsche Forschungsgemeinschaft. WL is additionally supported by a Fellowship of the China Scholarship Council. The work of DK is supported by the LIMTECH Research Alliance which is funded by the Helmholtz Association. We thank Jonathan Aurnou, Detlef Lohse and in particular Bruno Eckhardt for helpful discussions.

*

## Appendix A Hartmann layer

The Hartmann problem (33) describes an isothermal pressure-driven plane Poiseuille channel flow subject to a vertical homogeneous magnetic field (see also (20)). Starting point is equation (II.1) for . One seeks a steady solution in the quasistatic regime. This results in the inhomogeneous differential equation

 ρ0νd2ux(z)dz2−σB20ux(z)=−G, (37)

with The Hartmann layer thickness (17) arises as the characteristic length scale in the problem and is given by

 δv=√ρ0νσB20=HHa. (38)

### References

1. L. P. Kadanoff, Phys. Today 54 (8), 34 (2001).
2. G. Ahlers, S. Grossmann, and D. Lohse, Rev. Mod. Phys. 81, 503 (2009).
3. F. Chillà and J. Schumacher, Eur. Phys. J. E 35, article no. 58 (2012).
4. W. V. R. Malkus, Proc. R. Soc. London, Ser. A 225, 185 (1954).
5. W. V. R. Malkus, Proc. R. Soc. London, Ser. A 225, 196 (1954).
6. B. I. Shraiman and E. D. Siggia, Phys. Rev. A 28, 3650 (1990).
7. E. D. Siggia, Annu. Rev. Fluid Mech. 26, 137 (1994).
8. S. Grossmann and D. Lohse, J. Fluid Mech. 407, 27 (2000).
9. S. Grossmann and D. Lohse, Phys. Rev. Lett. 86, 3316 (2001).
10. Y. Yang, E. P. van der Poel, R. Ostilla-Monico, C. Sun, R. Verzicco, S. Grossmann and D. Lohse, J. Fluid Mech. 768, 476 (2015).
11. O. Shishkina, S. Grossmann and D. Lohse, Geophys. Res. Lett. 43, 1219 (2016).
12. N. O. Weiss and M. R. E. Procter, Magnetoconvection, Cambridge University Press, Cambridge, 2014.
13. M. Rempel and R. Schlichenmaier, Living Rev. Sol. Phys. 8, article no. 3, 60 pages (2011).
14. A. J. Castro-Tirado et al., Nature 455, 506 (2008).
15. R. W. Series and D. T. J. Hurle, J. Cryst. Growth 113, 305 (1991).
16. N. Shevchenko, O. Roshchupkina, O. Sokolova, and S. Eckert, J. Cryst. Growth 417, 1 (2015).
17. T. Ihli et al., Fusion Eng. Design 83, 912 (2008).
18. S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover, New York, 1961
19. Y. Tasaka, K. Igaki, T. Yanagisawa, T. Vogt, T. Zürner, and S. Eckert, Phys. Rev. E 93, 043109 (2016).
20. P. A. Davidson, An Introduction to Magnetohydrodynamics, Cambridge University Press, Cambridge, 2008.
21. S. Chakraborty, Physica D 237, 3233 (2008).
22. R. M. Kulsrud, Annu. Rev. Astron. Astrophys. 37, 37 (1999).
23. S. Cioni, S. Chaumat and J. Sommeria, Phys. Rev. E 62, R4520 (2000).
24. U. Burr and U. Müller, Phys. Fluids 13, 3247 (2001).
25. J. M. Aurnou, P. L. Olson, J. Fluid Mech. 430, 283 (2001).
26. D. Krasnov, O. Zikanov and T. Boeck, Comput. Fluids 50, 46 (2011).
27. J. J. Moré, The Levenberg-Marquardt algorithm: Implementation and theory, Springer, Berlin, Heidelberg 1978
28. R. J. A. M. Stevens, E. P. van der Poel, S. Grossmann and D. Lohse, J. Fluid Mech. 730, 295 (2013).
29. B. Castaing, G. Gunaratne, F. Heslot, L. Kadanoff, A. Libchaber, S. Thomae, X. Z. Wu, S. Zaleski and G. Zanetti, J. Fluid Mech. 204, 1 (1989).
30. J. D. Scheel and J. Schumacher, J. Fluid Mech. 802, 147 (2016).
31. T. Mashiko, Y. Tsuji, T. Mizuno, and M. Sano, Phys. Rev. E 69, 036306 (2004).
32. J. Schumacher, P. Götzfried, and J. D. Scheel, Proc. Natl. Acad. Sci. USA 112, 9530 (2015).
33. J. Hartmann, K. Dan. Vidensk. Selks. Mat. Phys. Medd. 15(7), 1 (1937).
121789