Next Article in Journal
The Energy Performances of a Ground-to-Air Heat Exchanger: A Comparison Among Köppen Climatic Areas
Next Article in Special Issue
Energy and Exergy Analysis of Using Turbulator in a Parabolic Trough Solar Collector Filled with Mesoporous Silica Modified with Copper Nanoparticles Hybrid Nanofluid
Previous Article in Journal
Simple Power Quality Compensation with Bidirectional Battery Charger for Electric Vehicles in Single-Phase Three-Wire Distribution Feeders
Previous Article in Special Issue
Study on Surface Condensate Water Removal and Heat Transfer Performance of a Minichannel Heat Exchanger
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Influence of Azimuthal Magnetic Field

Department of Aeronautics and Astronautics, Tokyo Metropolitan University, Hino 191-0065, Japan
*
Author to whom correspondence should be addressed.
Energies 2020, 13(11), 2896; https://doi.org/10.3390/en13112896
Submission received: 27 April 2020 / Revised: 26 May 2020 / Accepted: 3 June 2020 / Published: 5 June 2020
(This article belongs to the Special Issue Convection Process and Entropy Generation in Different Fluids)

Abstract

:
Natural convection of liquid metal in an annular enclosure under the influence of azimuthal static magnetic field 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 field 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 104, 5 × 105 and 107, 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.

1. 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 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].
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.

2. Methods

2.1. 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 κ = 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
r a v e = r i n + r o u t 2 = 1 + κ 2 ( 1 κ ) h ,
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
b θ ( r ) = μ m I 2 π r .
The reference value of magnetic field b0 was defined at the location of averaged radius.
b 0 ( r = r a v e ) = μ m I 2 π r a v e = μ m I π ( r i n + r o u t )
where, I is the constant electric current through the electric coil placed at the annulus center, and μm is the magnetic permeability of fluid.

2.2. Governing Equations and Dimensionless Numbers

In considering the governing equations for magnetohydrodynamic (MHD) natural convection in the enclosure, the following assumptions were made.
  • 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.
  • 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.
  • 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.
  • The effect of viscous dissipation and Joule heat are neglected compared to the externally given temperature from the walls.
The governing equations employed in this study are summarized as follows. Where, T0 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.
u = 0
u t + ( u ) u = 1 ρ 0 p + ν 2 u + 1 ρ 0 j × b + g β ( T T 0 ) e z
T t + ( u ) T = α 2 T
j = 0
j = σ ( ψ e + u × b )
In Equation (8), the static electric potential ψe was employed to satisfy Faraday’s law. The non-dimensional variables are as follows:
R , Z = r , z h , τ = t h / u 0 , U , V , W = u , v , w u 0 , P = p ρ 0 u 0 2 , Θ = T T 0 T h T c B θ = b θ b 0 = 1 + κ 2 ( 1 κ ) 1 R , J R , J θ , J Z = j r , j θ , j z σ u 0 b 0 , Ψ e = ψ e u 0 b 0 h .
where, the characteristic velocity, u0, was defined by the characteristic length, h, and the thermal diffusivity, α, as follows:
u 0 = α h .
The dimensionless governing equations expressed in the three-dimensional cylindrical coordinate system are summarized as follows.
U R + U R + 1 R V θ + W Z = 0
U τ + U U R + V R U θ + W U Z V 2 R = P R + P r ( 2 U R 2 + 1 R U R U R 2 + 1 R 2 2 U θ 2 2 R 2 V θ + 2 U Z 2 ) 1 + κ 2 ( 1 κ ) P r H a 2 R J Z
V τ + U V R + V R V θ + W V Z + U V R = 1 R P θ + P r ( 2 V R 2 + 1 R V R V R 2 + 1 R 2 2 V θ 2 + 2 R 2 U θ + 2 V Z 2 )
W τ + U W R + V R W θ + W W Z = P Z + P r ( 2 W R 2 + 1 R W R + 1 R 2 2 W θ 2 + 2 W Z 2 ) + 1 + κ 2 ( 1 κ ) P r H a 2 R J R + R a P r Θ
Θ τ + U Θ R + V R Θ θ + W Θ Z = 2 Θ R 2 + 1 R Θ R + 1 R 2 2 Θ θ 2 + 2 Θ Z 2
J R R + J R R + 1 R J θ θ + J Z Z = 0
J R = Ψ e R 1 + κ 2 ( 1 κ ) W R , J θ = 1 R Ψ e θ , J Z = Ψ e Z + 1 + κ 2 ( 1 κ ) U R
Non-dimensional numbers appearing in this study are the Prandtl number, Grashof number, Rayleigh number, and Hartmann number.
P r = ν α , G r = g β ( T h T c ) h 3 ν 2 , R a = g β ( T h T c ) h 3 α ν = G r P r , H a = b 0 h σ ρ 0 ν
The boundary conditions for the velocity, pressure, electric current density, electric potential, and temperature are as follows.
Hot wall ( R = κ 1 κ ) : U = V = W = P R = J R = Ψ e R = 0 , Θ = 0.5 Cold wall ( R = 1 1 κ ) : U = V = W = P R = J R = Ψ e R = 0 , Θ = 0.5 Adiabatic walls ( Z = 0 , 1 ) : U = V = W = P Z = J Z = Ψ e Z = Θ Z = 0
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.

2.3. 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
R i + 1 / 2 R i n = 1 2 { 1 + 1 a tanh [ ( i 1 N 1 2 ) ln ( 1 + a 1 a ) ] } ,
where Rin 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 Ri+1/2Rin = 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].
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.
C C L i , j n + 1 = J R i + 1 / 2 , j n + 1 J R i 1 / 2 , j n + 1 Δ R i + J R i + 1 / 2 , j n + 1 + J R i 1 / 2 , j n + 1 2 R i + J Z i , j + 1 / 2 n + 1 J Z i , j 1 / 2 n + 1 Δ Z j
J R i + 1 / 2 , j n + 1 = ( Ψ e i + 1 , j Ψ e i , j Δ R i + 1 / 2 1 + κ 2 ( 1 κ ) W c n t R i + 1 / 2 ) n + 1
J Z i , j + 1 / 2 n + 1 = ( Ψ e i , j + 1 Ψ e i , j Δ Z j + 1 / 2 + 1 + κ 2 ( 1 κ ) U c n t R i ) n + 1
W c n t = Δ R i + 1 ( Δ Z j + 1 / 2 W i , j 1 / 2 + Δ Z j 1 / 2 W i , j + 1 / 2 ) + Δ R i ( Δ Z j + 1 / 2 W i + 1 , j 1 / 2 + Δ Z j 1 / 2 W i + 1 , j + 1 / 2 ) ( Δ R i + Δ R i + 1 ) ( Δ Z j 1 / 2 + Δ Z j + 1 / 2 )
U c n t = Δ Z j + 1 ( Δ R i + 1 / 2 U i 1 / 2 , j + Δ R i 1 / 2 U i 1 / 2 , j + 1 ) + Δ Z j ( Δ R i + 1 / 2 U i + 1 / 2 , j + Δ R i 1 / 2 U i + 1 / 2 , j + 1 ) ( Δ Z j + Δ Z j + 1 ) ( Δ R i 1 / 2 + Δ R i + 1 / 2 )
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 Section 3.1, Section 3.2 and Section 3.3. The largest time step at which the computation did not diverge was chosen from 1 × 10n, 2 × 10n, and 5 × 10n, where n is an integer.

2.4. 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, JR and JZ, 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 JR, JZ, and Ψe is explained, where the superscript on the left indicates the iteration number. At the first step, JR and JZ were obtained temporarily.
J 1 R i + 1 / 2 , j n + 1 = Ψ 0 e i + 1 , j n + 1 Ψ 0 e i , j n + 1 Δ R i + 1 / 2 1 + κ 2 ( 1 κ ) W c n t n + 1 R i + 1 / 2
J 1 Z i , j + 1 / 2 n + 1 = Ψ 0 e i , j + 1 n + 1 Ψ 0 e i , j n + 1 Δ Z j + 1 / 2 + 1 + κ 2 ( 1 κ ) U c n t n + 1 R i
Next, Ψe is corrected.
Ψ m + 1 e i , j n + 1 = Ψ m e i , j n + 1 + δ m Ψ e i , j n + 1 = Ψ m e i , j n + 1 C m C L i , j n + 1 ( C C L i , j / Ψ e i , j ) m n + 1 = Ψ m e i , j n + 1 J m R i + 1 / 2 , j n + 1 J m R i 1 / 2 , j n + 1 Δ R i + J m R i + 1 / 2 , j n + 1 + J m R i 1 / 2 , j n + 1 2 R i + J m Z i , j + 1 / 2 n + 1 J m Z i , j 1 / 2 n + 1 Δ Z j 1 Δ R i Δ R i 1 / 2 + 1 Δ R i Δ R i + 1 / 2 + 1 Δ Z j Δ Z j 1 / 2 + 1 Δ Z j Δ Z j + 1 / 2
Then, JR and JZ are regarded as a function of Ψe and are corrected by Equations (29)–(32).
J m + 1 R i + 1 / 2 , j n + 1 = J m R i + 1 / 2 , j n + 1 + ( J R i + 1 / 2 , j Ψ e i , j ) n + 1 m δ m Ψ e i , j n + 1 = J m R i + 1 / 2 , j n + 1 + 1 Δ R i + 1 / 2 δ m Ψ e i , j n + 1
J m + 1 R i 1 / 2 , j n + 1 = J m R i 1 / 2 , j n + 1 1 Δ R i 1 / 2 δ m Ψ e i , j n + 1
J m + 1 Z i , j + 1 / 2 n + 1 = J m Z i , j + 1 / 2 n + 1 + 1 Δ Z j + 1 / 2 δ m Ψ e i , j n + 1
J m + 1 Z i , j 1 / 2 n + 1 = J m Z i , j 1 / 2 n + 1 1 Δ Z j 1 / 2 δ m Ψ e i , j n + 1
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 L2 norm of the energy equation was set as 10−8 or less of the L2 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.

2.5. Definition of Nusselt Number

In cylindrical coordinate, the ratio of local heat fluxes, Q = Qconv/Qcond, was defined, where Qconv was the local heat flux by the convective heat transfer as Equation (33) and Qcond 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. Qcond is in inversely proportion to the local radius R.
Q c o n v ( R , Z ) = U Φ Φ / R
Q c o n d ( R ) = 1 R ln κ
The definition of the average Nusselt number conformed to that by Davis [21].
N u R = 0 1 Q ( R , Z ) d Z
N u m = 0 1 N u R d R
N u v = N u R ( R = R a v e )
N u h = N u R ( R = R i n )
To calculate Qconv, 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,
Φ i + 1 / 2 Δ R i + 1 Φ i + Δ R i Φ i + 1 Δ R i + Δ R i + 1 .
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:
( Φ R ) i + 1 / 2 Φ i + Φ i + 1 Δ R i + 1 / 2 .
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
( Φ R ) i + 1 / 2 a 3 Φ i 1 b 3 Φ i + b 3 Φ i + 1 a 3 Φ i + 2 2 a b ( b a ) ( b + a ) ,
where
a = Δ R i + 1 / 2 2 , b = Δ R i 1 / 2 + Δ R i + 1 / 2 + Δ R i + 3 / 2 2 .
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, NuR. In this formula, to obtain an integration value, the functional value at a midpoint on a grid side, Qconv i+1/2,j, was multiplied by the length of a grid side, ΔZj.
N u R = 0 1 Q c o n v d Z j Q c o n v i + 1 / 2 , j Δ 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 Qconv i+1/2,j on a midpoint must be obtained with the second order.

3. Verification

3.1. 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, 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 = 103, 104, 105, and 106, and then the error, ε, was obtained.
ε = N u p r e s e n t N u D a v i s N u D a v i s
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 = 106. The error was generally less than 1% with N = 20, 40, 80, and 160 at Ra = 103, 104, 105, and 106, respectively.

3.2. Unsteady Characteristics with Uniform Grids

Time development characteristics were verified at Pr = 0.025 and Ra = 5 × 105. 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 × 105, 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 × 105.

3.3. 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 = 107, 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 = 107. 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 = 107. 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 = 107 for convenience of the computational time.

3.4. 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.
U = 1 R Ψ Z , W = 1 R Ψ R , J R = 1 R Ψ c Z , J Z = 1 R Ψ c R
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 JR 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 density vectors J , radial component of electric current density JR, and axial component of electric current density JZ at Pr = 0.025, κ = 0.5, and Ra = 104. 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.

4. Results

4.1. Steady Flows (Ra = 104)

The natural convection was steady at Pr = 0.025, κ = 0.5, and Ra = 104. 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 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.

4.2. Periodic Oscillation (Ra = 5 × 105)

Figure 9 shows the temporal evolution of Nuh at Pr = 0.025, κ = 0.5, and Ra = 5 × 105 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 Nuh were not different even by 1%.
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 = 104 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 × 105, and that between P and Ψc was also recognized at high Hartmann numbers as shown in the patterns in Ha ≥ 1000.

4.3. Turbulence (Ra = 107)

The natural convection changed into a turbulent regime oscillating irregularly at low Hartmann numbers for Ra = 107. 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 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 = 104, 5 × 105, and 107. The symmetry or the anti-symmetry did not appear when the natural convection oscillated irregularly at Ra = 107. 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.

4.4. 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> 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 = 107 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/Ra1/2. The three different curves are almost coincident by using <Nuh–1>/<Nuh (Ha = 0)–1> and Ha/Ra1/2. In [4], a combination of the Grashof number and Hartmann number, Ha/Gr1/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/Gr1/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 = 107. However, when the diagram is enlarged for the case of Ra = 107 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.
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.

4.5. 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.
J R = Ψ e R + 1 + κ 2 ( 1 κ ) 1 R 2 Ψ R 0 , J Z = Ψ e Z + 1 + κ 2 ( 1 κ ) 1 R 2 Ψ Z 0
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.
0 = P R + 1 + κ 2 ( 1 κ ) P r H a 2 R 2 Ψ c R , 0 = P Z + 1 + κ 2 ( 1 κ ) P r H a 2 R 2 Ψ c Z
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
2 Ψ e R 2 + 1 R Ψ e R + 2 Ψ e Z 2 = 1 + κ 2 ( 1 κ ) 1 R 2 ( 2 Ψ R 2 1 R Ψ R + 2 Ψ Z 2 ) .
It is recognized that there is similarity between Ψ and Ψe without explicitly using JR and JZ.
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.
U n + 1 = U n + Δ τ [ P R + 1 + κ 2 ( 1 κ ) P r H a 2 R 2 Ψ c R ] n W n + 1 = W n + Δ τ [ P Z + 1 + κ 2 ( 1 κ ) P r H a 2 R 2 Ψ c Z ] n
These equations are substituted into Equation (11), and we get the Poisson equation.
2 P R 2 + 1 R P R + 2 P Z 2 = 1 Δ τ ( U n R + U n R + W n Z ) + 1 + κ 2 ( 1 κ ) P r H a 2 R 2 ( 2 Ψ c R 2 1 R Ψ c R + 2 Ψ c Z 2 )
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.
2 P R 2 + 1 R P R + 2 P Z 2 = 1 + κ 2 ( 1 κ ) P r H a 2 R 2 ( 2 Ψ c R 2 1 R Ψ c R + 2 Ψ c Z 2 )

4.6. 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.
Θ τ + U Θ R + V R Θ θ + W Θ Z = 2 Θ R 2 + 1 R Θ R + 1 R 2 2 Θ θ 2 + 2 Θ Z 2 + Φ e m
The dimensionless Joule heat Φem is defined as
Φ e m = P r H a 2 E c J 2 ,
where Eckert number Ec = 10−8 and J2 = JR2 + Jθ2 + JZ2. The dimensionless spatio-temporal average Joule heat <Φem>m is defined in 2D as follows.
Φ e m m = Φ e m d τ d r d z d τ d r d z
Since ∂Θ/∂τ~Φem, <Φem>m means the Joule heat per unit time in the enclosure. All the values examined for <Φem>m were less than 10−2 as shown in Figure 16. The Joule heat was hence ignored.

5. 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 = 104, 5 × 105, and 107, 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 × 105, 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, <Nuh–1>, versus Ha was divided into three regions on a double logarithmic graph. <Nuh–1> was almost constant at low Hartmann numbers, and it was inversely proportional to about Ha2 and converged to <Nuh> → 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 = 107, 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.

Author Contributions

Conceptualization, T.T.; methodology, T.M. and T.T.; software, T.M.; validation, T.M. and T.T.; formal analysis, T.M.; investigation, T.M.; data curation, T.M.; writing—original draft preparation, T.M.; writing—review and editing, T.T.; visualization, T.M.; supervision, T.T.; project administration, T.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

b magnetic flux density = (0, bθ, 0) (T)
B dimensionless magnetic flux density = (0, Bθ, 0) (-)
b0value of magnetic flux density at r = rave (T)
e z unit vector in the axial direction
ggravitational acceleration (m/s2)
GrGrashof number (-)
hwidth or height length of duct (m)
HaHartmann number (-)
Iconstant electric current (A)
j electric current density = (jr, 0, jz) (A/m2)
J dimensionless electric current density = (JR, 0, JZ) (-)
Nnumber of grids in the R- or Z-direction (-)
Nuhaverage Nusselt number at r = rin (-)
Numaverage Nusselt number within whole domain (-)
Nuvaverage Nusselt number at r = rave. (-)
NuRaverage Nusselt number at r = const. (-)
ppressure (Pa)
Pdimensionless pressure (-)
PrPrandtl number (-)
Qlocal heat flux (J/(m2·s))
Qcondlocal conductive heat flux (J/(m2·s))
Qconvlocal convective heat flux (J/(m2·s))
rr coordinate (m)
Rdimensionless r coordinate (-)
RaRayleigh number (-)
raveaverage radius (m)
rinradius at inner wall (m)
routradius at outer wall (m)
ttime (s)
Ttemperature (K)
Tctemperature at outer (cold) wall (K)
Thtemperature at inner (hot) wall (K)
T0reference temperature = (Th + Tc)/2 (K)
ur-directional velocity component (m/s)
Udimensionless r-directional velocity component (-)
u0reference velocity (m/s)
vθ-directional velocity component (m/s)
Vdimensionless θ-directional velocity component (-)
wz-directional velocity component (m/s)
Wdimensionless z-directional velocity component (-)
zz coordinate (m)
Zdimensionless z coordinate (-)
Greek symbols
αthermal diffusivity (m2/s)
βvolumetric coefficient of expansion at T0 (1/K)
Δdiscrete variant
θθ coordinate (rad)
Θdimensionless temperature (-)
κradius ratio (-)
μmmagnetic permeability (H/m)
νkinematic viscosity (m2/s)
ρdensity (kg/m3)
σelectric conductivity (1/(Ω·m))
τdimensionless time (-)
Φdimensionless temperature = Θ + 0.5 (-)
ψeelectric potential (V)
Ψdimensionless stream function of velocity (-)
Ψcdimensionless stream function of electric current density (-)
Ψedimensionless electric potential (-)
Subscripts or superscripts
cntlinear interpolation
i, jgrid point
mnumber of iterations
ntime step
Brackets
< >time-mean

References

  1. Molokov, S.; Moreau, R.; Moffatt, H.K. Fluid Mechanics and Its Application. In Magnetohydrodynamics: Historical Evolution and Trends; Springer: Berlin/Heidelberg, Germany, 2007; Volume 80. [Google Scholar]
  2. Smolentsev, S.; Moreau, R.; Bühler, L.; Mistrangelo, C. MHD thermofluid issues of liquid-metal blankets: Phenomena and advances. Fusion Eng. Des. 2010, 85, 1196–1205. [Google Scholar] [CrossRef]
  3. Fumizawa, M. Natural convection experiment with liquid NaK under transverse magnetic field. J. Nucl. Sci. Tech. 1980, 17, 98–105. [Google Scholar] [CrossRef]
  4. Okada, K.; Ozoe, H. Experimental heat transfer rates of natural convection of molten gallium suppressed under an external magnetic field in either the X, Y, or Z direction. J. Heat Transf. 1992, 114, 107–114. [Google Scholar] [CrossRef]
  5. Ozoe, H.; Okada, K. The effect of the direction of the external magnetic field on the three-dimensional natural convection in a cubical enclosure. Int. J. Heat Mass Transf. 1989, 32, 1939–1954. [Google Scholar] [CrossRef]
  6. Tagawa, T.; Ozoe, H. Enhanced heat transfer rate measured for natural convection in liquid gallium in a cubical enclosure under a static magnetic field. J. Heat Transf. 1998, 120, 1027–1032. [Google Scholar] [CrossRef]
  7. Tagawa, T.; Ozoe, H. Enhancement of heat transfer rate by application of a static magnetic field during natural convection of liquid metal in a cube. J. Heat Transf. 1997, 119, 265–271. [Google Scholar] [CrossRef]
  8. Authié, G.; Tagawa, T.; Moreau, R. Buoyant flow in long vertical enclosures in the presence of a strong horizontal magnetic field. Part 2. Finite enclosures. Eur. J. Mech.—B/Fluids 2003, 22, 203–220. [Google Scholar] [CrossRef]
  9. Wang, Z.H.; Meng, X.; Ni, M.J. Liquid metal buoyancy driven convection heat transfer in a rectangular enclosure in the presence of a transverse magnetic field. Int. J. Heat Mass Transf. 2017, 113, 514–523. [Google Scholar] [CrossRef]
  10. Tagawa, T.; Ozoe, H. The natural convection of liquid metal in a cubical enclosure with various electro-conductivities of the wall under the magnetic field. Int. J. Heat Mass Transf. 1998, 41, 1917–1928. [Google Scholar] [CrossRef]
  11. Tagawa, T.; Ozoe, H. Combined Effect of the magnetic strengths and electro-conductivity of the enclosure wall on the natural convection of liquid metal in a cubical enclosure. J. Mater. Process. Manuf. Sci. 2002, 10, 137–146. [Google Scholar] [CrossRef]
  12. Di Piazza, I.; Ciofalo, M. MHD free convection in a liquid-metal filled cubic enclosure. I. Differential heating. Int. J. Heat Mass Transf. 2002, 45, 1477–1492. [Google Scholar] [CrossRef]
  13. Gajbhiye, N.L.; Eswaran, V. MHD buoyant flow in a cubical enclosure at low to high Hartmann number. Int. J. Therm. Sci. 2018, 134, 168–178. [Google Scholar] [CrossRef]
  14. Moreau, R.J. Fluid Mechanics and Its Application. In Magnetohydrodynamics; Kluwer Academic Publishers: Dordrecht, The Netherlands; Norwell, MA, USA, 1990; Volume 3. [Google Scholar]
  15. Kakarantzas, S.C.; Sarris, I.E.; Vlachos, N.S. Natural convection of liquid metal in a vertical annulus with lateral and volumetric heating in the presence of a horizontal magnetic field. Int. J. Heat Mass Transf. 2011, 54, 3347–3356. [Google Scholar] [CrossRef]
  16. Kakarantzas, S.C.; Benos, L.T.; Sarris, I.E.; Knaepen, B.; Grecos, A.P.; Vlachos, N.S. MHD liquid metal flow and heat transfer between vertical coaxial cylinders under horizontal magnetic field. Int. J. Heat Fluid Flow 2017, 65, 342–351. [Google Scholar] [CrossRef]
  17. Afrand, M.; Rostami, S.; Akbari, M.; Wongwises, S.; Esfe, M.H.; Karimipour, A. Effect of induced electric field on magneto-natural convection in a vertical cylindrical annulus filled with liquid potassium. Int. J. Heat Mass Transf. 2015, 90, 418–426. [Google Scholar] [CrossRef]
  18. Mahfoud, B.; Bessaïh, R. Magnetohydrodynamic counter-rotating flow in a cylindrical cavity. Int. J. Heat Mass Transf. 2016, 93, 175–185. [Google Scholar] [CrossRef]
  19. Wang, W.; Li, B.W.; Varghese, P.L.; Leng, X.Y.; Tian, X.Y. Numerical analysis of three-dimensional MHD natural convection flow in a short horizontal cylindrical annulus. Int. Comm. Heat Mass Trans. 2018, 98, 273–285. [Google Scholar] [CrossRef]
  20. Ashour, R.F.; Kelley, D.H.; Salas, A.; Starace, M.; Weber, N.; Weier, T. Competing forces in liquid metal electrodes and batteries. J. Power Sources 2018, 378, 301–310. [Google Scholar] [CrossRef] [Green Version]
  21. Beltrán, A. MHD natural convection flow in a liquid metal electrode. Appl. Therm. Eng. 2017, 114, 1203–1212. [Google Scholar] [CrossRef]
  22. Kelley, D.H.; Weier, T. Fluid Mechanics of Liquid Metal Batteries. Appl. Mech. Rev. 2017, 70. [Google Scholar] [CrossRef] [Green Version]
  23. Weber, N.; Nimtz, M.; Personnettaz, P.; Salas, A.; Weier, T. Electromagnetically driven convection suitable for mass transfer enhancement in liquid metal batteries. Appl. Therm. Eng. 2018, 143, 293–301. [Google Scholar] [CrossRef] [Green Version]
  24. Weber, N.; Nimtz, M.; Personnettaz, P.; Weier, T.; Sadoway, D. Numerical simulation of mass transfer enhancement in liquid metal batteries by means of electro-vortex flow. J. Power Sources Adv. 2020, 1, 100004. [Google Scholar] [CrossRef]
  25. Tagawa, T.; Ozoe, H. Effect of Prandtl number and computational schemes on the oscillatory natural convection in an enclosure. Numer. Heat Transf. Part A Appl. 1996, 30, 271–282. [Google Scholar] [CrossRef]
  26. De Vahl Davis, G. Natural convection of air in a square cavity: A bench mark numerical solution. Int. J. Numer. Methods Fluids 1983, 3, 249–264. [Google Scholar] [CrossRef]
  27. Feigenbaum, M.J. The transition to aperiodic behavior in turbulent systems. Commun. Math. Phys. 1980, 77, 65–86. [Google Scholar] [CrossRef]
  28. Hof, B.; Juel, A.; Mullin, T. Magnetohydrodynamic damping of convective flows in molten gallium. J. Fluid Mech. 2003, 482, 163–179. [Google Scholar] [CrossRef] [Green Version]
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.
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.
Energies 13 02896 g001
Figure 2. Staggered grid system.
Figure 2. Staggered grid system.
Energies 13 02896 g002
Figure 3. Dependence of the error, ε, to the grid number, N at Pr = 0.71.
Figure 3. Dependence of the error, ε, to the grid number, N at Pr = 0.71.
Energies 13 02896 g003aEnergies 13 02896 g003b
Figure 4. Time development of Nuh with uniform grids at Pr = 0.025, κ = 0.99, Ra = 5 × 105, and Ha = 0.
Figure 4. Time development of Nuh with uniform grids at Pr = 0.025, κ = 0.99, Ra = 5 × 105, and Ha = 0.
Energies 13 02896 g004
Figure 5. Time development of Nuh with non-uniform grids at Pr = 0.025, κ = 0.5, Ra = 107, and Ha = 0.
Figure 5. Time development of Nuh with non-uniform grids at Pr = 0.025, κ = 0.5, Ra = 107, and Ha = 0.
Energies 13 02896 g005
Figure 6. Dependency of time-mean averaged Nusselt numbers, <Num>, <Nuv>, and <Nuh> to grid number, N, at Pr = 0.025, κ = 0.5, Ra = 107, and Ha = 0.
Figure 6. Dependency of time-mean averaged Nusselt numbers, <Num>, <Nuv>, and <Nuh> to grid number, N, at Pr = 0.025, κ = 0.5, Ra = 107, and Ha = 0.
Energies 13 02896 g006
Figure 7. The contour map of streamline of electric current density Ψc, electric current density vectors J , radial component of electric current density JR, and axial component of electric current density JZ at Pr = 0.025, κ = 0.5, and Ra = 104.
Figure 7. The contour map of streamline of electric current density Ψc, electric current density vectors J , radial component of electric current density JR, and axial component of electric current density JZ at Pr = 0.025, κ = 0.5, and Ra = 104.
Energies 13 02896 g007
Figure 8. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 104. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Figure 8. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 104. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Energies 13 02896 g008
Figure 9. Time development of Nuh at κ = 0.5 and Ra = 5 × 105.
Figure 9. Time development of Nuh at κ = 0.5 and Ra = 5 × 105.
Energies 13 02896 g009
Figure 10. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 5 × 105. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Figure 10. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 5 × 105. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Energies 13 02896 g010
Figure 11. Time development of Nuh at κ = 0.5 and Ra = 107.
Figure 11. Time development of Nuh at κ = 0.5 and Ra = 107.
Energies 13 02896 g011
Figure 12. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 107. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Figure 12. Contour figures for several values of Hartmann number at κ = 0.5 and Ra = 107. From top to bottom, stream function of velocity, pressure, temperature, stream function of electric current density, and electric potential are shown.
Energies 13 02896 g012
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. 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.
Energies 13 02896 g013
Figure 14. Diagram of <Nuh–1>/<Nuh (Ha = 0)–1> to Ha/Ra1/2 at κ = 0.5.
Figure 14. Diagram of <Nuh–1>/<Nuh (Ha = 0)–1> to Ha/Ra1/2 at κ = 0.5.
Energies 13 02896 g014
Figure 15. Enlarged diagram of time-mean average Nusselt number versus Hartmann number at κ = 0.5 and Ra = 107.
Figure 15. Enlarged diagram of time-mean average Nusselt number versus Hartmann number at κ = 0.5 and Ra = 107.
Energies 13 02896 g015
Figure 16. Dimensionless spatio-temporal average Joule heat <Φem>m, which was small enough to be ignored.
Figure 16. Dimensionless spatio-temporal average Joule heat <Φem>m, which was small enough to be ignored.
Energies 13 02896 g016
Table 1. Grid numbers and time steps.
Table 1. Grid numbers and time steps.
RaGrid TypeNΔτ
1 × 104Uniform401 × 10−3
5 × 105Uniform1205 × 10−5
1 × 107Non-uniform801 × 10−5

Share and Cite

MDPI and ACS Style

Masuda, T.; Tagawa, T. Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Influence of Azimuthal Magnetic Field. Energies 2020, 13, 2896. https://doi.org/10.3390/en13112896

AMA Style

Masuda T, Tagawa T. Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Influence of Azimuthal Magnetic Field. Energies. 2020; 13(11):2896. https://doi.org/10.3390/en13112896

Chicago/Turabian Style

Masuda, Takuya, and Toshio Tagawa. 2020. "Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Influence of Azimuthal Magnetic Field" Energies 13, no. 11: 2896. https://doi.org/10.3390/en13112896

APA Style

Masuda, T., & Tagawa, T. (2020). Axisymmetric Natural Convection of Liquid Metal in an Annular Enclosure under the Influence of Azimuthal Magnetic Field. Energies, 13(11), 2896. https://doi.org/10.3390/en13112896

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop