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
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.
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.
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,
u0, 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.
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.
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
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/2 −
Rin = 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.
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 × 10
n, 2 × 10
n, and 5 × 10
n, 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.
Then,
JR and
JZ 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 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.
The definition of the average Nusselt number conformed to that by Davis [
21].
To calculate
Qconv, the temperature,
Φ, and the temperature gradient,
, 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,
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.
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.
4. Results
4.1. Steady Flows (Ra = 104)
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 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 × 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
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 = 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.
4.3. Turbulence (Ra = 107)
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 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.
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 = 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/
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 = 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.
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.
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, / 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
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.
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.
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.
The dimensionless Joule heat
Φem is defined as
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.
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.