Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Inﬂuence of Azimuthal Magnetic Field

: Natural convection of liquid metal in an annular enclosure under the inﬂuence of azimuthal static magnetic ﬁeld was numerically studied. The liquid metal in the enclosure whose cross-sectional area is square was heated from an inner vertical wall and cooled from an outer vertical wall both isothermally whereas the other two horizontal walls were assumed to be adiabatic. The static azimuthal magnetic ﬁeld was imposed by a long straight electric coil that was located at the central axis of the annular enclosure. The computations were carried out for the Prandtl number 0.025, the Rayleigh number 10 4 , 5 × 10 5 and 10 7 , and the Hartmann number 0–100,000 by using an in-house code. It was found that the contour map of the electric potential was similar to that of the Stokes stream function of the velocity regardless of the Hartmann number. Likewise, the contour map of the pressure was similar to the Stokes stream function of the electric current density in the case of the high Hartmann number. The average Nusselt number was decreased in proportion to the square of the Hartmann number in the high Hartmann number regime.


Introduction
Natural convection of liquid metal under an external magnetic field is a phenomenon encountered in liquid blankets of a thermonuclear fusion reactor, metallurgy process, manufacturing single crystals of semiconductors, and so on [1,2].This is an important theme not only for improving the quality of such industrial products but also for the development of magnetohydrodynamics (MHD).One of the characteristics of liquid metal is that it has very high thermal conductivity and therefore it is classified as low-Prandtl-number fluid.Second, since it has high electric conductivity, it is possible to control the flow by applying a magnetic field from the outside.Third, it is usually opaque, which makes it difficult to observe the inside of flow field.The natural convection of liquid metal having such properties is affected by the inertial force because of its low Prandtl number characteristics, so it easily tends to be oscillatory or to exhibit turbulent in the high Rayleigh number regime.On the other hand, many previous studies have revealed that the flow field of natural convection is decelerated and stabilized by applying an external magnetic field.Hence, it is considered that the heat transfer rate is also reduced together with damping of the flow field.
However, contrary to this common recognition, in Fumizawa's experiment [3] using NaK in a rectangular enclosure, it was reported that the heat transfer was greatly increased by applying a static magnetic field.Okada and Ozoe [4,5] carried out experiments using liquid gallium in a cubic enclosure as well as three-dimensional numerical computations with Pr = 0.054.When the direction of the uniform magnetic field and the vortex axis of natural convection were parallel to each other (this case Energies 2020, 13, 2896 2 of 23 is called as Y-directional magnetic field in their papers), it was reported that the phenomena are significantly different from those of the other two cases (X-directional magnetic field or Z-directional magnetic field).Of particular note, the experiments suggested that heat transfer was enhanced slightly at a lower Hartmann number regime.However, in their experiment, since the size of the container was 30 mm, experiments with a higher Rayleigh number were not able to be performed and the confirmation was not obtained.Later, Tagawa and Ozoe [6] changed the container size to 64 mm and made experiments with a higher Rayleigh number with applying a uniform magnetic field in the Y-direction, and an obvious heat transfer enhancement phenomenon was observed.
In the experiment of Ozoe's group [4,6], the heat transfer coefficient (Nusselt number) was obtained by measuring the temperature at several points on the wall heated with uniform heat flux.However, the flow field was not observed because liquid gallium is an opaque fluid.Numerical calculation is an effective method to know the state of the liquid metal flow field.Tagawa and Ozoe [7] carried out a three-dimensional numerical computation of the natural convection of a liquid metal in a cubic enclosure and showed that the Nusselt number took its maximum value under a weak magnetic field.It was reported that in the absence of a magnetic field, the main vortex of natural convection was disturbed due to flow instability, but this disturbance tended to be suppressed as increase in the strength of the Y-direction magnetic field.Authié et al. [8] performed mercury experiments and the corresponding numerical computations of natural convection in a vertically long rectangular enclosure (aspect ratio 7.5) having a horizontal square cross-section.It was shown that the maximum Nusselt number was obtained in a certain weak magnetic field regime, and their experiments and computational results agreed well with each other qualitatively.In their experiments, the flow velocity of the core flow was estimated by measuring the electric potential on the Hartmann wall (the wall perpendicular to the applied uniform magnetic field) with using their own-developed probes.This is based on the fact that the flow field becomes quasi-two-dimensional when the magnetic field strength exceeds a certain level, and the stream line and the electric potential have similar tendency between them.In another experimental study, Wang et al. [9] used an Ultrasound Doppler velocimetry to measure the velocity of natural convection of liquid metal.Since their experiments were almost the same as the system of Tagawa and Ozoe [6,7], the Nusselt number took the maximum value under a weak magnetic field, and the Nusselt number decreased in a stronger range of magnetic field.
It has been reported that the natural convection of liquid metal under an external magnetic field is strongly influenced not only by the strength of the imposed magnetic field but also the electric conductivity of the container.Tagawa and Ozoe [10,11] performed numerical analyses for a cubic container wall whose conductivity ranged from zero to infinity, and showed that as the conductivity of the container wall increased, the decrease of Nusselt number became larger.Di Piazza and Ciofalo [12] performed numerical computations on a system similar to Tagawa and Ozoe [11].In their calculations, instead of having grids inside the solid container, the effect of the conducting wall was incorporated using the thin-wall approximation and a wall function.In recent years, Gajbhiye and Eswaran [13] have discussed the influence of the magnetic field direction for the similar system.In summary, in the case of natural convection of liquid metal in a rectangular container, the effect of enhancement of heat transfer by applying a magnetic field would be obtained when the conductivity of the container wall is small (insulating wall) and the Rayleigh number should be rather high, and likely be limited to the case where the axis of the main vortex of the convection and the applied magnetic field are parallel to each other.In such a case, all the current generated in the core region comes into the Hartmann layer, so the influence of the electric potential is large, which is also related to the quasi-two-dimensionality of the flow field.Generally, the Hartmann layer is a boundary layer formed in the vicinity of the wall perpendicular to the applied magnetic field direction, and its thickness is known to be inversely proportional to the Hartmann number [14].Therefore, the Hartmann layer becomes very thin, especially under the strong magnetic field, which makes numerical calculation difficult.One of the methods to avoid this difficulty is to employ the wall function, which was actually attempted in [6,10].
Energies 2020, 13, 2896 3 of 23 Natural convection of liquid metal under a magnetic field is of great engineering importance not only in a rectangular container but also in a cylindrical container in relation to the single crystal growth process such as the Czochralski method or to the liquid blanket of a fusion reactor.Kakarantzas et al. [15] carried out a numerical computation in which a horizontal magnetic field was applied to natural convection in a vertical annulus.The aspect ratio of the vertical cross-section of the container was fixed, and the effects of temperature difference between the inner and outer walls and internal volumetric heating were discussed.Later, the same group discussed the effects of radius ratio and aspect ratio only when there was no internal volumetric heating [16].Afrand et al. [17] discussed the importance of the induced electric field for a system similar to that of Kakarantzas et al.Mahfoud and Bessaïh [18] studied the effect of an axial magnetic field when the upper and lower disks rotated in a vertical cylindrical cavity.In that case, the bottom disk was heated and the top disk was cooled, both isothermally, resulting in a combined convection under the magnetic field.Wang et al. [19] numerically computed the natural convection in a horizontal cylindrical annulus whose inner and outer cylinders were heated and cooled, respectively, under the presence of an axial magnetic field.
As mentioned above, there have been several studies related to MHD natural convection in a cylindrical container, but it has been limited to the case that a uniform magnetic field parallel or perpendicular to the central axis of the container has been applied.A series of recent studies [20][21][22][23][24], motivated by liquid metal batteries, have indeed considered the combined effects of natural convection and azimuthal magnetic fields due to an axial current.Therefore, in this study, we attempted to consider natural convection of liquid metal in an annular vessel having a rectangular cross-sectional area under the influence of an azimuthal magnetic field.

Model Configuration
The configuration of the present analysis is shown in Figure 1.An annular enclosure with a square cross-section is filled with an electric conducting fluid under the influence of uniform gravitational field.The fluid is assumed to be an incompressible, Newtonian, and low-Prandtl-number fluid like usual liquid metals such as molten iron, mercury, or gallium.Unless otherwise specified, the Prandtl number, Pr, is limited to 0.025 in the present study.The inner wall of the enclosure is heated, and the opposing outer wall is cooled both isothermally whereas the upper and lower walls are thermally insulated.As a result of the interaction between the gravitational field and the fluid density difference, natural convection takes place due to buoyancy.We employed the Boussinesq approximation.The static azimuthal magnetic field is imposed by an infinitely long straight electric coil placed at the central axis of the annular enclosure.Both buoyancy and Lorentz force act on the working fluid as the external forces.The axisymmetric convection could occur with respect to the electric current coil placed at the central axis when the strength of the magnetic field is sufficient to stabilize the instability of liquid metal natural convection.
The radius ratio is defined as κ = r in /r out , where r in is the distance measured from the annulus center to the inner wall of the enclosure, and r out is that to the outer wall.In the present analysis, the radius ratio was set to κ = 0.5.The radius from the annulus center to the cross-section center of the enclosure is expressed as where the characteristic length, h, is the width or the height of the cross-section of enclosure.The magnetic flux density has only the circumferential component and it is derived from the Ampère's law as The reference value of magnetic field b 0 was defined at the location of averaged radius.
where, I is the constant electric current through the electric coil placed at the annulus center, and µ m is the magnetic permeability of fluid.
Energies 2020, 13, x FOR PEER REVIEW 4 of 25 Figure 1.A schematic model.The inner wall is heated and the outer wall is cooled, both isothermally, while the top and bottom horizontal walls are adiabatic.Static azimuthal magnetic field is imposed by an electric coil placed at the central axis of the annular enclosure.The cross-sectional area is supposed to be square and the characteristic length is h.
The radius ratio is defined as κ = rin/rout, where rin is the distance measured from the annulus center to the inner wall of the enclosure, and rout is that to the outer wall.In the present analysis, the radius ratio was set to κ = 0.5.The radius from the annulus center to the cross-section center of the enclosure is expressed as ( ) where the characteristic length, h, is the width or the height of the cross-section of enclosure.The magnetic flux density has only the circumferential component and it is derived from the Ampère's law as The reference value of magnetic field b0 was defined at the location of averaged radius.
where, I is the constant electric current through the electric coil placed at the annulus center, and μm is the magnetic permeability of fluid.

Governing Equations and Dimensionless Numbers
In considering the governing equations for magnetohydrodynamic (MHD) natural convection in the enclosure, the following assumptions were made.
A schematic model.The inner wall is heated and the outer wall is cooled, both isothermally, while the top and bottom horizontal walls are adiabatic.Static azimuthal magnetic field is imposed by an electric coil placed at the central axis of the annular enclosure.The cross-sectional area is supposed to be square and the characteristic length is h.

Governing Equations and Dimensionless Numbers
In considering the governing equations for magnetohydrodynamic (MHD) natural convection in the enclosure, the following assumptions were made.

1.
MHD approximation: since the displacement current or the convection current is much smaller than the conduction current, only the conduction current is taken into account.By this assumption, the electric charge conservation law is satisfied.2.
Ignoring induced magnetic field: electric current density is generated by the movement of liquid metal (the magnetic Prandtl number is very small, and its value is the order of 10 −6 ) under an external magnetic field, but the induced magnetic field created by this electric current density can be ignored compared to the strength of the external magnetic field.In this situation, the Ampère's law is no longer needed, and the current density must be obtained from Ohm's law.

3.
The Boussinesq approximation is assumed to be hold in the present study and the other physical properties such as viscosity, thermal conductivity, specific heat, thermal expansion coefficient, magnetic permeability, and electric conductivity are the constants irrespective of the fluid temperature.The values of such physical properties are taken at the reference temperature.

4.
The effect of viscous dissipation and Joule heat are neglected compared to the externally given temperature from the walls.
Energies 2020, 13, 2896 5 of 23 The governing equations employed in this study are summarized as follows.Where, T 0 is the reference temperature, which is taken as the average temperature between the hot and cold walls, and ρ 0 is the fluid density at the reference temperature.
In Equation ( 8), the static electric potential ψ e was employed to satisfy Faraday's law.The non-dimensional variables are as follows: where, the characteristic velocity, u 0 , was defined by the characteristic length, h, and the thermal diffusivity, α, as follows: The dimensionless governing equations expressed in the three-dimensional cylindrical coordinate system are summarized as follows.
Non-dimensional numbers appearing in this study are the Prandtl number, Grashof number, Rayleigh number, and Hartmann number.
Energies 2020, 13, 2896 6 of 23 The boundary conditions for the velocity, pressure, electric current density, electric potential, and temperature are as follows.
In this study, the natural convection was assumed to be uniform in the azimuthal direction, which was similar to the 2D calculation of the natural convection in a square duct at Pr = 0.025 [20] or at Pr = 0.71 [21].The result at Pr = 0.71 [21] has been regarded as the benchmark solution, and has been used for validation of CFD software.

Discretization in Axisymmetric Cases
The staggered grid system utilized in this study is shown in Figure 2. The grid number N is the number of cells in a coordinate direction.The grid number in the R direction and that in the Z direction were divided to be equal to each other.Either in a uniform or a non-uniform spacing grid system was used depending on the value of the Rayleigh number.The velocity was discretized at each interface of meshes.When the non-uniform spacing grid system was used, the grid points at the velocity definition in the R direction were set as where R in was the dimensionless radius from the annulus center to the inside wall.A similar equation was used for those in the Z direction.Equation (20) was only applied inside the fluid region including the wall surface.In order to set the boundary conditions at the walls, virtual grids were defined so as to be symmetrical with respect to the internal grids defined in the fluid domain.A double grid system was used in the outside the wall so that R i+1/2 − R in = 0 at i = 1.The definition points of the pressure were the midpoint between the definition points of the velocity.The grid space can be changed by the coefficient a, which can be selected from the range of 0 < a < 1.The closer a is to 0, the closer the grid is to the uniform spacing grid, and the closer a is to 1, the denser the grid is near the wall.In this study, a = 0.9 was used unless otherwise specified, when the minimum grid space of the non-uniform grid was about 1/3 of that of the uniform grid in the condition that both grid numbers were same.This rate is nearly the same as that used by Tagawa and Ozoe [25].
Energies 2020, 13, x FOR PEER REVIEW 7 of 25 of the pressure were the midpoint between the definition points of the velocity.The grid space can be changed by the coefficient a, which can be selected from the range of 0 < a < 1.The closer a is to 0, the closer the grid is to the uniform spacing grid, and the closer a is to 1, the denser the grid is near the wall.In this study, a = 0.9 was used unless otherwise specified, when the minimum grid space of the non-uniform grid was about 1/3 of that of the uniform grid in the condition that both grid numbers were same.This rate is nearly the same as that used by Tagawa and Ozoe [25].The third-order upwind difference method called UTOPIA scheme was used for the convection terms.Spacing discretization except the convection terms was approximated with the second-order central difference method.The conservation of the electric charge (Equation ( 16)) and Ohm's law (Equation ( 17)) were discretized as follows.The third-order upwind difference method called UTOPIA scheme was used for the convection terms.Spacing discretization except the convection terms was approximated with the second-order central difference method.The conservation of the electric charge (Equation ( 16)) and Ohm's law (Equation ( 17)) were discretized as follows.
The superscript indicates a time step, and the subscript does the definition point on the staggered grid, respectively.With respect to time derivative schemes, the second-order Runge-Kutta method for the momentum balance and the second-order Adams-Bashforth and Crank-Nicholson method for the energy equation were used, respectively.The combinations of the grid number, N, and a time step, ∆τ, for the Rayleigh number shown in Table 1 were used in most cases.The validity of the grid number was inspected in Sections 3.1-3.3.The largest time step at which the computation did not diverge was chosen from 1 × 10 n , 2 × 10 n , and 5 × 10 n , where n is an integer.Non-uniform 80 1 × 10 −5

Simultaneous Solution
The numerical simultaneous solution is as follows.At first, the velocity and the pressure at a certain time step were obtained simultaneously by the Newton-Raphson method, which is called the HSMAC algorithm.Furthermore, the temperature was obtained by the successive over-relaxation (SOR) method.After that, the electric current density, J R and J Z , and the electric potential, Ψ e , were calculated by the Newton-Raphson method, which was similar to the HSMAC algorithm for the velocity and the pressure.
The algorithm to obtain J R , J Z , and Ψ e is explained, where the superscript on the left indicates the iteration number.At the first step, J R and J Z were obtained temporarily.
Then, J R and J Z are regarded as a function of Ψ e and are corrected by Equations ( 29)-(32).
The iteration used from Equations ( 28)-( 32) was repeated until a convergence condition was reached.After that, the calculation was proceeded to a next time step.The boundary condition was considered at the wall.The convergence conditions were defined as follows.For the pressure and the velocity, the maximum value of divergence of the velocity was set as 10 −3 or less.For the electric potential and the electric current density, that of the electric current density was set as 10 −3 or less, and on the temperature, the L 2 norm of the energy equation was set as 10 −8 or less of the L 2 norm just before a time step discretized by the Crank-Nicholson method.The over-relaxation coefficient, ω, can be set for the Newton-Raphson method or the SOR method.In this study, ω D = 1.7 for the continuity of the mass, ω T = 1.5 for the energy equation, and ω E = 1.9 for the conservation of the electric charge.

Definition of Nusselt Number
In cylindrical coordinate, the ratio of local heat fluxes, Q = Q conv /Q cond , was defined, where Q conv was the local heat flux by the convective heat transfer as Equation (33) and Q cond was that by the heat conduction as Equation (34).The thermal field was transformed from −0.5 ≤ Θ ≤ 0.5 into 0 ≤ Φ ≤ 1, where Φ = Θ + 0.5.Q cond is in inversely proportion to the local radius R.
The definition of the average Nusselt number conformed to that by Davis [21].
Energies 2020, 13, 2896 9 of 23 To calculate Q conv , the temperature, Φ, and the temperature gradient, ∂Φ/∂R, were approximated as follows.Φ was obtained by the linear interpolation at the definition point of U, that is, It is understood by the Taylor expansion that this linear interpolation is the second order.∂Φ/∂R was obtained by the second-order central difference used two points on both sides: The truncation error of Equation ( 40) on a non-uniform grid is the first order, but the truncation error decreases with the second order accuracy if there are enough grid points.On the wall, the fourth-order central difference with four points on both sides was used as where The fourth-order central difference is considered to be nearly the second-order on the wall.In the benchmark solution by Davis [26], the second-order difference approximation was also used for the temperature gradient including on the wall.
The midpoint formula for a numerical integration was used to get the average Nusselt number, Nu R .In this formula, to obtain an integration value, the functional value at a midpoint on a grid side, Q conv i+1/2,j , was multiplied by the length of a grid side, ∆Z j .
Equation ( 43) is the second order accuracy if the functional value is known on a midpoint.To keep this accuracy, the approximation of Q conv i+1/2,j on a midpoint must be obtained with the second order.

Comparison with the Benchmark Solution in Cartesian Coordinates
The natural convection in Cartesian coordinates was computed by eliminating the inherent terms in cylindrical coordinates on the programing code.The results with our code were compared with the benchmark solution by Davis [21], so that the reliability of the code was verified.When the grid number increased gradually, it was confirmed that the error of the present results to the benchmark value estimated by Davis decreased monotonously.The grid number was defined as the number of the region divided.In regard to the time schemes, only in Section 3.1, the second-order Adams-Bashforth method for the momentum balance and the second-order Runge-Kutta method for the energy equation were used, respectively.
To calculate the Nusselt number, the Cartesian coordinate X was introduced instead of the cylindrical coordinate R. In a square duct, Q cond = 1, and Q conv can be obtained to change R into X in Equation (33).Nu v was defined at X = 0.5, and Nu h was done at X = 0.The computation was implemented at Ra = 10 3 , 10 4 , 10 5 , and 10 6 , and then the error, ε, was obtained.The transition of the error with respect to the grid number is shown in Figure 3.The error decreased exponentially with increasing the grid number except the result with N = 10 at Ra = 10 6 .The error was generally less than 1% with N = 20, 40, 80, and 160 at Ra = 10 3 , 10 4 , 10 5 , and 10 6 , respectively.cylindrical coordinate R. In a square duct, Qcond = 1, and Qconv can be obtained to change R into X in Equation (33).Nuv was defined at X = 0.5, and Nuh was done at X = 0.The computation was implemented at Ra = 10 3 , 10 4 , 10 5 , and 10 6 , and then the error, ε, was obtained.

− =
The transition of the error with respect to the grid number is shown in Figure 3.The error decreased exponentially with increasing the grid number except the result with N = 10 at Ra = 10 6 .The error was generally less than 1% with N = 20, 40, 80, and 160 at Ra = 10 3 , 10 4 , 10 5 , and 10 6 , respectively.

Unsteady Characteristics with Uniform Grids
Time development characteristics were verified at Pr = 0.025 and Ra = 5 × 10 5 .In this condition, Tagawa and Ozoe [25] investigated grid numbers and schemes of a convection term, and then reported that periodic oscillation flows occurred.
The temporal evolution of Nuh to the non-dimensional time, τ, are shown in Figure 4, where Pr = 0.025, κ = 0.99, Ra = 5 × 10 5 , and Ha = 0.It seems that the characteristics at κ = 0.99 are almost the same as those in a square duct.In the uniform grids, one-period simple oscillation was observed with N = 80, but two-period oscillation was recognized with N ≥ 120, which included double period oscillation in addition to fundamental oscillation.In the non-uniform grids, two-period oscillation was obtained with N ≥ 80.
The introduction of a non-uniform grid system has an advantage for the reduction of grid numbers, but the grid number cannot be reduced much because a certain grid number is required even in a core region of the convection for an unsteady flow.In addition, the time step must be also smaller for the limit of the Courant number or the diffusion number if a minimum grid interval becomes narrower.A total computational time in the uniform grid was less than that in the non-uniform grid.The uniform grid system was used in Ra ≤ 5 × 10

Unsteady Characteristics with Uniform Grids
Time development characteristics were verified at Pr = 0.025 and Ra = 5 × 10 5 .In this condition, Tagawa and Ozoe [25] investigated grid numbers and schemes of a convection term, and then reported that periodic oscillation flows occurred.
The temporal evolution of Nu h to the non-dimensional time, τ, are shown in Figure 4, where Pr = 0.025, κ = 0.99, Ra = 5 × 10 5 , and Ha = 0.It seems that the characteristics at κ = 0.99 are almost the same as those in a square duct.In the uniform grids, one-period simple oscillation was observed with N = 80, but two-period oscillation was recognized with N ≥ 120, which included double period oscillation in addition to fundamental oscillation.In the non-uniform grids, two-period oscillation was obtained with N ≥ 80.
The introduction of a non-uniform grid system has an advantage for the reduction of grid numbers, but the grid number cannot be reduced much because a certain grid number is required even in a core region of the convection for an unsteady flow.In addition, the time step must be also smaller for the limit of the Courant number or the diffusion number if a minimum grid interval becomes narrower.A total computational time in the uniform grid was less than that in the non-uniform grid.The uniform grid system was used in Ra ≤ 5 × 10 5 .

Unsteady Characteristics with Non-Uniform Grids
A gradient of a physical value is larger near the walls than that in the core region.If the uniform grid is used at Ra = 10 7 , based on the verification so far, the grid number is required at least N = 320.It is not realistic from the view point of a computational time, so the non-uniform grid system was employed at Ra = 10 7 .In order to verify that the characteristics of the flow did not depend on the number of the grids, the time development of Nu h was investigated and shown in Figure 5 for Pr = 0.025, κ = 0.5, and Ha = 0.The natural convection changed into a turbulent state oscillating irregularly at Ra = 10 7 .As long as the grid number was 60 ≤ N ≤ 120, the unsteady characteristics did not seem to be affected by the grid number.Grid dependence was also investigated by comparing the time-mean average Nusselt number, which was averaged for 0.5 non-dimensional time.The time-mean variables were expressed as < >.As shown in Figure 6, the relative error of <Nu m >, <Nu v >, and <Nu h > was within 1% although the data should have a certain error itself.In this simulation, the non-uniform grid of N = 80 was used at Ra = 10 7 for convenience of the computational time.
Time development characteristics were verified at Pr = 0.025 and Ra = 5 × 10 5 .In this condition, Tagawa and Ozoe [25] investigated grid numbers and schemes of a convection term, and then reported that periodic oscillation flows occurred.
The temporal evolution of Nuh to the non-dimensional time, τ, are shown in Figure 4, where Pr = 0.025, κ = 0.99, Ra = 5 × 10 5 , and Ha = 0.It seems that the characteristics at κ = 0.99 are almost the same as those in a square duct.In the uniform grids, one-period simple oscillation was observed with N = 80, but two-period oscillation was recognized with N ≥ 120, which included double period oscillation in addition to fundamental oscillation.In the non-uniform grids, two-period oscillation was obtained with N ≥ 80.
The introduction of a non-uniform grid system has an advantage for the reduction of grid numbers, but the grid number cannot be reduced much because a certain grid number is required even in a core region of the convection for an unsteady flow.In addition, the time step must be also smaller for the limit of the Courant number or the diffusion number if a minimum grid interval becomes narrower.A total computational time in the uniform grid was less than that in the non-uniform grid.The uniform grid system was used in Ra ≤ 5 × 10

Unsteady Characteristics with Non-Uniform Grids
A gradient of a physical value is larger near the walls than that in the core region.If the uniform grid is used at Ra = 10 7 , based on the verification so far, the grid number is required at least N = 320.It is not realistic from the view point of a computational time, so the non-uniform grid system was employed at Ra = 10 7 .In order to verify that the characteristics of the flow did not depend on the number of the grids, the time development of Nuh was investigated and shown in Figure 5 for Pr = 0.025, κ = 0.5, and Ha = 0.The natural convection changed into a turbulent state oscillating irregularly at Ra = 10 7 .As long as the grid number was 60 ≤ N ≤ 120, the unsteady characteristics did not seem to be affected by the grid number.Grid dependence was also investigated by comparing the time-mean average Nusselt number, which was averaged for 0.5 non-dimensional time.The time-mean variables were expressed as < >.As shown in Figure 6, the relative error of <Num>, <Nuv>, and <Nuh> was within 1% although the data should have a certain error itself.In this simulation, the non-uniform grid of N = 80 was used at Ra = 10 7 for convenience of the computational time.

Stream Function of Electric Current Density
Stokes stream function of the electric current density in a cylindrical coordinate, Ψc, was defined by changing the velocity components, U and W, into the electric current density components, JR and JZ, for stream function of the velocity, Ψ, as follows.

Unsteady Characteristics with Non-Uniform Grids
A gradient of a physical value is larger near the walls than that in the core region.If the uniform grid is used at Ra = 10 7 , based on the verification so far, the grid number is required at least N = 320.It is not realistic from the view point of a computational time, so the non-uniform grid system was employed at Ra = 10 7 .In order to verify that the characteristics of the flow did not depend on the number of the grids, the time development of Nuh was investigated and shown in Figure 5 for Pr = 0.025, κ = 0.5, and Ha = 0.The natural convection changed into a turbulent state oscillating irregularly at Ra = 10 7 .As long as the grid number was 60 ≤ N ≤ 120, the unsteady characteristics did not seem to be affected by the grid number.Grid dependence was also investigated by comparing the time-mean average Nusselt number, which was averaged for 0.5 non-dimensional time.The time-mean variables were expressed as < >.As shown in Figure 6, the relative error of <Num>, <Nuv>, and <Nuh> was within 1% although the data should have a certain error itself.In this simulation, the non-uniform grid of N = 80 was used at Ra = 10 7 for convenience of the computational time.

Stream Function of Electric Current Density
Stokes stream function of the electric current density in a cylindrical coordinate, Ψc, was defined by changing the velocity components, U and W, into the electric current density components, JR and JZ, for stream function of the velocity, Ψ, as follows.

Stream Function of Electric Current Density
Stokes stream function of the electric current density in a cylindrical coordinate, Ψ c , was defined by changing the velocity components, U and W, into the electric current density components, J R and J Z , for stream function of the velocity, Ψ, as follows.
The value of Ψ c was obtained by numerical integration based on the electrically insulating wall (Ψ c = 0) in order to draw a contour map.Although the direction of the numerical integration could be selected from four directions as R+, R−, Z+, and Z−, it was confirmed that the results were not different regardless of the choice.The numerical integration of J R was performed in the Z+ direction from Z = 0 there.
Figure 7 shows the contour map of streamline of electric current density Ψ c , electric current

Steady Flows (Ra = 10 4 )
The natural convection was steady at Pr = 0.025, κ = 0.5, and Ra = 10 4 .Figure 8 shows the contour figures of Ψ, P, Θ, Ψc, and Ψe at Ha = 10, 100, 200, and 1000.The contour patterns at Ha = 0 were almost the same as those at Ha = 10.Therefore, those at Ha = 0 were omitted.The natural convection was damped by the effect of the Lorentz force at high Hartmann numbers.As a result, the pattern of Θ was closed to the heat conduction state with increasing the value of Ha.In regard to the vortex of Ψ and Ψe, the center of those vortices moved from the center of the duct to near the cold wall, and the shape of those vortices changed from a circular to a flat shape.The patterns of Ψc were generally anti-symmetric with respect to the horizontal center line.The center of those vortices moved to near the thermally insulated walls upside and downside with an increase of Ha.The patterns of P were circular having a low pressure region in the enclosure center at Ha = 0, but tended to become

Steady Flows (Ra = 10 4 )
The natural convection was steady at Pr = 0.025, κ = 0.5, and Ra = 10 4 .Figure 8 shows the contour figures of Ψ, P, Θ, Ψ c , and Ψ e at Ha = 10, 100, 200, and 1000.The contour patterns at Ha = 0 were almost the same as those at Ha = 10.Therefore, those at Ha = 0 were omitted.The natural convection was damped by the effect of the Lorentz force at high Hartmann numbers.As a result, the pattern of Θ was closed to the heat conduction state with increasing the value of Ha.In regard to the vortex of Ψ and Ψ e , Energies 2020, 13, 2896 13 of 23 the center of those vortices moved from the center of the duct to near the cold wall, and the shape of those vortices changed from a circular to a flat shape.The patterns of Ψ c were generally anti-symmetric with respect to the horizontal center line.The center of those vortices moved to near the thermally insulated walls upside and downside with an increase of Ha.The patterns of P were circular having a low pressure region in the enclosure center at Ha = 0, but tended to become anti-symmetry similar to Ψ c due to the generation of a high pressure region at Ha = 1000.It is interesting that the distributions of Ψ and those of Ψ e were similar to each other.This tendency was observed in all the cases of Ha.Ψ is zero on the wall by a no-slip boundary condition of the velocity, but Ψ e does not have such a restriction so that Ψ e is not constant on the wall.Furthermore, the patterns of the contour figures also showed the similarity between P and Ψ c at high Hartmann numbers.Both distributions were very similar to each other at Ha = 1000.The reason for these phenomena is described in Section 4.5.As regards a transition process, the flow was the one-period oscillation in Ha ≤ 50, but appeared a two-period oscillation at Ha = 100, and became the one-period oscillation at Ha = 200 again.It was suspected to occur by a grid dependence.However, even investigating with the uniform grid of N = 160 at Ha = 0, the one-period oscillation was obtained, and therefore the grid dependence was not detected.With an increase of Ha, at least at κ = 0.5, it is concluded that the time development characteristic changed in the order of one-period oscillation, two-period oscillation, one-period oscillation, and a steady state.Destabilization of a flow promotes period doubling bifurcation.Such a transition process has been known as Feigenbaum's scenario [27].It is interesting that the two-period oscillation was observed only at Ha = 100, where the stabilization effect by imposing a magnetic field should be stronger than at Ha = 0 or 50. Figure 10 shows the contour figures of time mean variables, where the symbol noted by < > represents time mean value.The patterns changed in the same tendency as that at Ra = 10 4 with respect to Ha, but asymmetry appeared especially near the corners.The similarity between Ψ and Ψe was also recognized at all Ha for Ra = 5 × 10 5 , and that between P and Ψc was also recognized at high Hartmann numbers as shown in the patterns in Ha ≥ 1000.As regards a transition process, the flow was the one-period oscillation in Ha ≤ 50, but appeared a two-period oscillation at Ha = 100, and became the one-period oscillation at Ha = 200 again.It was suspected to occur by a grid dependence.However, even investigating with the uniform grid of N = 160 at Ha = 0, the one-period oscillation was obtained, and therefore the grid dependence was not detected.With an increase of Ha, at least at κ = 0.5, it is concluded that the time development characteristic changed in the order of one-period oscillation, two-period oscillation, one-period oscillation, and a steady state.Destabilization of a flow promotes period doubling bifurcation.Such a transition process has been known as Feigenbaum's scenario [27].It is interesting that the two-period oscillation was observed only at Ha = 100, where the stabilization effect by imposing a magnetic field should be stronger than at Ha = 0 or 50.
Figure 10 shows the contour figures of time mean variables, where the symbol noted by < > represents time mean value.The patterns changed in the same tendency as that at Ra = 10 4 with respect to Ha, but asymmetry appeared especially near the corners.The similarity between Ψ and Ψ e was also recognized at all Ha for Ra = 5 × 10 5 , and that between P and Ψ c was also recognized at high Hartmann numbers as shown in the patterns in Ha ≥ 1000.

Turbulence (Ra = 10 7 )
The natural convection changed into a turbulent regime oscillating irregularly at low Hartmann numbers for Ra = 10 7 .The turbulence was damped by applying the magnetic field at high Hartmann numbers.Figure 11 shows the temporal evolution of Nuh.Certainly, the oscillation amplitude of natural convection was suppressed in the range of Ha ≤ 1000 but the Nuh was not influenced significantly.The higher Ha, the smaller the amplitude of Nuh.It became steady in Ha ≥ 5000.The irregular turbulent state was kept without a transition to a periodic oscillation.The natural convection changed into a turbulent regime oscillating irregularly at low Hartmann numbers for Ra = 10 7 .The turbulence was damped by applying the magnetic field at high Hartmann numbers.Figure 11 shows the temporal evolution of Nu h .Certainly, the oscillation amplitude of natural convection was suppressed in the range of Ha ≤ 1000 but the Nu h was not influenced significantly.The higher Ha, the smaller the amplitude of Nu h .It became steady in Ha ≥ 5000.The irregular turbulent state was kept without a transition to a periodic oscillation.
The contour figures of the time mean variables are shown in Figure 12.There was little difference in the patterns in Ha ≤ 100.The maximum value and/or the minimal value of the physical values other than temperature moved from the enclosure center to near the cold wall, which was common among Ra = 10 4 , 5 × 10 5 , and 10 7 .The symmetry or the anti-symmetry did not appear when the natural convection oscillated irregularly at Ra = 10 7 .The symmetry of the patterns may have a relationship to the irregularly of the flow oscillation.The similarity between Ψ and Ψ e was recognized at all Ha even when the natural convection was turbulent, and between P and Ψ c was also recognized in Ha ≥ 2000.In regard to the conditions where P and Ψ c were similar, it was common among the three Rayleigh numbers that the flow was steady or quasi-steady state at high Hartmann numbers.The contour figures of the time mean variables are shown in Figure 12.There was little difference in the patterns in Ha ≤ 100.The maximum value and/or the minimal value of the physical values other than temperature moved from the enclosure center to near the cold wall, which was common among Ra = 10 4 , 5 × 10 5 , and 10 7 .The symmetry or the anti-symmetry did not appear when the natural convection oscillated irregularly at Ra = 10 7 .The symmetry of the patterns may have a relationship to the irregularly of the flow oscillation.The similarity between Ψ and Ψe was recognized at all Ha even when the natural convection was turbulent, and between P and Ψc was also recognized in Ha ≥ 2000.In regard to the conditions where P and Ψc were similar, it was common among the three Rayleigh numbers that the flow was steady or quasi-steady state at high Hartmann numbers.

Variation of Nusselt Number
The time-mean average Nusselt numbers changed with respect to Ha as shown in Figure 13, which were plotted on a double logarithmic graph with the vertical axis as <Nuh-1>.The distribution curve was decomposed into three regions with regard to the value of Rayleigh number.<Nuh-1> was almost constant at low Hartmann numbers, and was proportional to about square of Ha so that it converged as Nuh → 1 at high Hartmann numbers.<Nuh-1> where the power region started depended on Ra.The transition from the low Ha to the high Ha was smooth, which was shown on the double logarithmic graph (Figure 14).These lines on a double logarithmic graph were smoothly connected decreasing to Ha at an increasing tempo.Furthermore, the distribution curves of <Nuh-1>

Variation of Nusselt Number
The time-mean average Nusselt numbers changed with respect to Ha as shown in Figure 13, which were plotted on a double logarithmic graph with the vertical axis as <Nu h -1>.The distribution curve was decomposed into three regions with regard to the value of Rayleigh number.<Nu h -1> was almost constant at low Hartmann numbers, and was proportional to about square of Ha so that it converged as Nu h → 1 at high Hartmann numbers.<Nu h -1> where the power region started depended on Ra.The transition from the low Ha to the high Ha was smooth, which was shown on the double logarithmic graph (Figure 14).These lines on a double logarithmic graph were smoothly connected decreasing to Ha at an increasing tempo.Furthermore, the distribution curves of <Nu h -1> were similarity relation with respect to Ra, in which the length and the starting/end point of the connection region and the power slope were slightly different.The distribution curve at Ra = 10 7 should have a similarity relation to the other curves.Figure 14 shows a re-arrangement of Figure 13 by modifying scales of vertical and horizontal axes.The vertical axis was modified so that the Nusselt number without the magnetic field takes the value of unity, i.e., <Nuh-1>/<Nuh (Ha = 0)-1> was used.While, the horizontal axis shows a combination of the Rayleigh and Hartmann numbers, Ha/Ra 1/2 .The three different curves are almost coincident by using <Nuh-1>/<Nuh (Ha = 0)-1> and Ha/Ra 1/2 .In [4], a combination of the Grashof number and Hartmann number, Ha/Gr 1/3 was used and most of experimental data plotted are coincident around a line.On the other hand, in [6,7], a combination of Ha/Gr 1/2 was proposed for the case that the applied magnetic field direction is parallel to the main convection flow.In the latter case [6,7], the magnetic damping effect is much weaker than the former case [4].Therefore, the present phenomena indicate that the magnetic damping effect is quite weak.
When a flow of natural convection is turbulent, the maximum of the Nusselt number may appear at Ha ≠ 0 because the turbulence tends to be suppressed by effect of the magnetic field.This unfamiliar phenomenon that takes maximum Nusselt number at Ha ≠ 0 has been already reported   Figure 14 shows a re-arrangement of Figure 13 by modifying scales of vertical and horizontal axes.The vertical axis was modified so that the Nusselt number without the magnetic field takes the value of unity, i.e., <Nuh-1>/<Nuh (Ha = 0)-1> was used.While, the horizontal axis shows a combination of the Rayleigh and Hartmann numbers, Ha/Ra 1/2 .The three different curves are almost coincident by using <Nuh-1>/<Nuh (Ha = 0)-1> and Ha/Ra 1/2 .In [4], a combination of the Grashof number and Hartmann number, Ha/Gr 1/3 was used and most of experimental data plotted are coincident around a line.On the other hand, in [6,7], a combination of Ha/Gr 1/2 was proposed for the case that the applied magnetic field direction is parallel to the main convection flow.In the latter case [6,7], the magnetic damping effect is much weaker than the former case [4].Therefore, the present phenomena indicate that the magnetic damping effect is quite weak.
When a flow of natural convection is turbulent, the maximum of the Nusselt number may appear at Ha ≠ 0 because the turbulence tends to be suppressed by effect of the magnetic field.This unfamiliar phenomenon that takes maximum Nusselt number at Ha ≠ 0 has been already reported Figure 14 shows a re-arrangement of Figure 13 by modifying scales of vertical and horizontal axes.The vertical axis was modified so that the Nusselt number without the magnetic field takes the value of unity, i.e., <Nu h -1>/<Nu h (Ha = 0)-1> was used.While, the horizontal axis shows a combination of the Rayleigh and Hartmann numbers, Ha/Ra 1/2 .The three different curves are almost coincident by using <Nu h -1>/<Nu h (Ha = 0)-1> and Ha/Ra 1/2 .In [4], a combination of the Grashof number and Hartmann number, Ha/Gr 1/3 was used and most of experimental data plotted are coincident around a line.On the other hand, in [6,7], a combination of Ha/Gr 1/2 was proposed for the case that the applied magnetic field direction is parallel to the main convection flow.In the latter case [6,7], the magnetic damping effect is much weaker than the former case [4].Therefore, the present phenomena indicate that the magnetic damping effect is quite weak.
When a flow of natural convection is turbulent, the maximum of the Nusselt number may appear at Ha 0 because the turbulence tends to be suppressed by effect of the magnetic field.This unfamiliar phenomenon that takes maximum Nusselt number at Ha 0 has been already reported in several references of the natural convection of liquid metal in rectangular enclosures [4,6,28].This phenomenon should also happen in the present annular duct.In Figure 13, it seems that the Nusselt number is nearly constant in the range of Ha < 1000 at Ra = 10 7 .However, when the diagram is enlarged for the case of Ra = 10 7 as shown in Figure 15, it is observed that the Nusselt number takes its maximum at around Ha = 500.The rate of increase was about 5% compared to the Nusselt number for Ha = 0.The appearance of a maximum Nusselt number at non-zero Hartmann number deserves a more detailed discussion.
Energies 2020, 13, x FOR PEER REVIEW 20 of 25 The present numerical computation assuming the axial symmetry cannot capture the three-dimensional turbulent flow as it should be.Nevertheless, axially symmetric temporally irregular oscillatory flows have turbulent components.By applying azimuthal magnetic field, convection is hardly suppressed, but turbulent components are suppressed more efficiently, and as a result, it is considered that the flow becomes smoother and the heat transfer coefficient increases.If the three-dimensional computation for the natural convection in the annular enclosure was performed, the turbulent component would be larger, so it could be expected that a more remarkable increase in the heat transfer coefficient could be observed.This is a future issue.

Similarity
As shown in the contour maps, there are similarities between the Stokes stream function of the velocity Ψ and the electric potential Ψe for all Hartmann numbers, and those between the pressure P and the Stokes stream function of electric current density Ψc for the case of high Hartmann numbers.In this section, we consider this reason, where the governing equations are assumed as uniform in the azimuthal direction.
First, let us explain a simple consideration for the similarity between Ψ and Ψe by using Ohm's law (Equation ( 17)).It is assumed that each component of electric current density is approximately equal to zero.From Equation (17), it is recognized that Ψ and Ψe are similar to each other.
On the other hand, the similarity between P and Ψc can be derived from the momentum Equations ( 12) and ( 14) in which the inertial, viscous, and buoyancy terms are neglected due to the dominance of electromagnetic force.Since we are discussing the time-averaged field, the time-derivative term was also ignored.
If there was no R −2 factor in the last term, ∂P / ∂R and ∂ Ψc/ ∂ R would be directly proportional to each other, assuring similar spatial structure.The presence of the R −2 causes different radial variation, but in this domain, h < rave, so the difference is small.On the other hand, near the origin, the similarity would be weak, since R −2 is steep there.
In Equation ( 46), JR and JZ approximated zero, but the order of JR and JZ is 1 at most.In order to explain the similarity more generally, we attempted to derive Poisson equations.If Equation ( 17) is substituted into Equation ( 16), we can obtain The present numerical computation assuming the axial symmetry cannot capture the three-dimensional turbulent flow as it should be.Nevertheless, axially symmetric temporally irregular oscillatory flows have turbulent components.By applying azimuthal magnetic field, convection is hardly suppressed, but turbulent components are suppressed more efficiently, and as a result, it is considered that the flow becomes smoother and the heat transfer coefficient increases.If the three-dimensional computation for the natural convection in the annular enclosure was performed, the turbulent component would be larger, so it could be expected that a more remarkable increase in the heat transfer coefficient could be observed.This is a future issue.

Similarity
As shown in the contour maps, there are similarities between the Stokes stream function of the velocity Ψ and the electric potential Ψ e for all Hartmann numbers, and those between the pressure P and the Stokes stream function of electric current density Ψ c for the case of high Hartmann numbers.In this section, we consider this reason, where the governing equations are assumed as uniform in the azimuthal direction.
First, let us explain a simple consideration for the similarity between Ψ and Ψ e by using Ohm's law (Equation ( 17)).It is assumed that each component of electric current density is approximately equal to zero.From Equation (17), it is recognized that Ψ and Ψ e are similar to each other.
On the other hand, the similarity between P and Ψ c can be derived from the momentum Equations ( 12) and ( 14) in which the inertial, viscous, and buoyancy terms are neglected due to the dominance of electromagnetic force.Since we are discussing the time-averaged field, the time-derivative term was also ignored.
If there was no R −2 factor in the last term, ∂P/ ∂R and ∂Ψ c / ∂R would be directly proportional to each other, assuring similar spatial structure.The presence of the R −2 causes different radial variation, but in this domain, h < r ave , so the difference is small.On the other hand, near the origin, the similarity would be weak, since R −2 is steep there.
In Equation (46), J R and J Z approximated zero, but the order of J R and J Z is 1 at most.In order to explain the similarity more generally, we attempted to derive Poisson equations.If Equation ( 17) is substituted into Equation ( 16), we can obtain It is recognized that there is similarity between Ψ and Ψ e without explicitly using J R and J Z .The Poisson equation for the pressure can be also derived from the momentum equations excluding the inertial, viscous, and buoyancy terms.The time-derivative terms are discretized by the Euler explicit method.
These equations are substituted into Equation ( 11), and we get the Poisson equation.
As long as the continuity of mass is satisfied, the first term of the right-hand-side becomes zero.The similarity between P and Ψ c were also recognized as follows.

Consideration on Joule Heat
The Joule heat is estimated to verify the MHD assumption that the Joule heat should have been small enough.The Joule heat is considered in the dimensionless energy equation as follows.
The dimensionless Joule heat Φ em is defined as where Eckert number Ec = 10 −8 and J 2 = J R 2 + J θ 2 + J Z 2 .The dimensionless spatio-temporal average Joule heat <Φ em > m is defined in 2D as follows.

onclusions
In regard to the natural convection of liquid metal in the annular enclosure under the influen zimuthal magnetic field, the two-dimensional axisymmetric time-marching computations we formed.The following conclusions were obtained for the Prandtl number, Pr = 0.025, the radi o, κ = 0.5, the Rayleigh number, Ra = 10 4 , 5 × 10 5 , and 10 7 , and the Hartmann numbers, 0 ≤ H ,000.The contour map between the stream function of the velocity, Ψ, and the electric potential, Ψe, artmann numbers was similar, and that between the pressure, P, and the stream function of t tric current density, Ψc, at high Hartmann numbers was similar too.The period-doubling bifurcation of the flow was generally promoted by the decrease in t tmann number, but an exceptional case that did not follow this general tendency was observe a = 5 × 10 5 , for example, the natural convection was the one-period oscillation in Ha ≤ 200 but t -period oscillation only at Ha = 100.The distribution curve of the time-averaged Nusselt number, <Nuh-1>, versus Ha was divid three regions on a double logarithmic graph.<Nuh-1> was almost constant at low Hartma bers, and it was inversely proportional to about Ha 2 and converged to <Nuh> → 1.The tw ions were connected smoothly in the middle Hartmann numbers.
At high Rayleigh numbers where the natural convection was turbulent, the maximum of t e-averaged Nusselt number appeared at Ha ≠ 0 due to the suppression of the turbulence by gnetic field.It had already been reported in the case of cubic enclosure.With respect to t ymmetric natural convection in the present annular duct, the time-averaged Nusselt numb

Conclusions
In regard to the natural convection of liquid metal in the annular enclosure under the influence of azimuthal magnetic field, the two-dimensional axisymmetric time-marching computations were performed.The following conclusions were obtained for the Prandtl number, Pr = 0.025, the radius ratio, κ = 0.5, the Rayleigh number, Ra = 10 4 , 5 × 10 5 , and 10 7 , and the Hartmann numbers, 0 ≤ Ha ≤ 100,000.
The contour map between the stream function of the velocity, Ψ, and the electric potential, Ψ e , at all Hartmann numbers was similar, and that between the pressure, P, and the stream function of the electric current density, Ψ c , at high Hartmann numbers was similar too.
The period-doubling bifurcation of the flow was generally promoted by the decrease in the Hartmann number, but an exceptional case that did not follow this general tendency was observed.At Ra = 5 × 10 5 , for example, the natural convection was the one-period oscillation in Ha ≤ 200 but the two-period oscillation only at Ha = 100.
The distribution curve of the time-averaged Nusselt number, <Nu h -1>, versus Ha was divided into three regions on a double logarithmic graph.<Nu h -1> was almost constant at low Hartmann numbers, and it was inversely proportional to about Ha 2 and converged to <Nu h > → 1.The two regions were connected smoothly in the middle Hartmann numbers.
At high Rayleigh numbers where the natural convection was turbulent, the maximum of the time-averaged Nusselt number appeared at Ha 0 due to the suppression of the turbulence by a magnetic field.It had already been reported in the case of cubic enclosure.With respect to the axisymmetric natural convection in the present annular duct, the time-averaged Nusselt number took its maximum at Ha = 500 for Ra = 10 7 , which slightly increased by about 5% to Ha = 0.
Further development of this research includes linear stability analysis as well as three-dimensional numerical simulation.By performing the linear stability analysis, the critical condition for the onset of three-dimensional convection would be elucidated depending on the input parameters.By carrying out the three-dimensional numerical simulation, enhancement of heat transfer at a moderate Hartmann number could be observed in the present configuration.

6 Figure 3 .
Figure 3. Dependence of the error, ε, to the grid number, N at Pr = 0.71.

J
, radial component of electric current density J R , and axial component of electric current density J Z at Pr = 0.025, κ = 0.5, and Ra = 10 4 .Among the four figures, the contour map of Ψ c is probably most understandable what is occurring in the cross-sectional area in the enclosure.The upper blue contour lines indicate clockwise circulation of electric current density, while the lower red ones indicate contour-clockwise circulation of electric current density.Since all the walls are electrically insulated, the radially inward current in the core region must be recirculated through the boundary layers formed near the top and bottom walls.
Energies 2020, 13, x FOR PEER REVIEW 13 of 25 is probably most understandable what is occurring in the cross-sectional area in the enclosure.The upper blue contour lines indicate clockwise circulation of electric current density, while the lower red ones indicate contour-clockwise circulation of electric current density.Since all the walls are electrically insulated, the radially inward current in the core region must be recirculated through the boundary layers formed near the top and bottom walls.

Figure 7 .
Figure 7.The contour map of streamline of electric current density Ψc, electric current density vectors J 

Figure 7 .
Figure 7.The contour map of streamline of electric current density Ψ c , electric current density vectors → J , radial component of electric current density J R , and axial component of electric current density J Z at Pr = 0.025, κ = 0.5, and Ra = 10 4 .

Figure 8 .Figure 9
Figure 8. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 10 4 .From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.

Figure 8 .
Figure 8. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 10 4 .From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.

4. 2 .Figure 9
Figure9shows the temporal evolution of Nu h at Pr = 0.025, κ = 0.5, and Ra = 5 × 10 5 for the various Hartmann numbers.The computations were carried out in the range of τ = 0-2.At Ha = 200, a longer time was needed to get the periodic state among the Hartmann numbers investigated.As a magnetic field contributes to stabilize the natural convection, the higher Ha, the smaller the amplitude of oscillation.The flow exhibited oscillation in Ha ≤ 200 and was steady in Ha ≥ 500.The amplitude changed little in Ha ≤ 50 but it became smaller in Ha ≥ 100.Comparing the cases of Ha = 50 and Ha = 100, the waveforms were clearly different, but the time mean values of Nu h were not different even by 1%.

Figure 10 .
Figure 10.Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 5 × 10 5 .From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.

Figure 10 .
Figure 10.Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 5 × 10 5 .From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.

Figure 12 .
Figure 12.Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 10 7 .From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.

Energies 2020 ,
13, x FOR PEER REVIEW 19 of 25 were similarity relation with respect to Ra, in which the length and the starting/end point of the connection region and the power slope were slightly different.The distribution curve at Ra = 10 7 should have a similarity relation to the other curves.

Figure 13 .
Figure 13.Diagram of <Nuh-1> to Ha at κ = 0.5.For convenience, the results at Ha = 0 were plotted on the double logarithmic graph, and so Ha = 0 is exactly zero.

Figure 13 .
Figure 13.Diagram of <Nu h -1> to Ha at κ = 0.5.For convenience, the results at Ha = 0 were plotted on the double logarithmic graph, and so Ha = 0 is exactly zero.

Figure 13 .
Figure 13.Diagram of <Nuh-1> to Ha at κ = 0.5.For convenience, the results at Ha = 0 were plotted on the double logarithmic graph, and so Ha = 0 is exactly zero.

Figure 16 .
Figure16.Dimensionless spatio-temporal average Joule heat <Φ em > m , which was small enough to be ignored.

Table 1 .
Grid numbers and time steps.