1. Introduction
A blanket is an indispensable component of a tokamak fusion reactor. It surrounds the reaction chamber and serves the triple purpose of shielding the magnets from the neutrons generated by the reaction, converting the neutrons’ energy into heat to be diverted into an external power cycle and breeding tritium fuel via interaction of neutrons with lithium atoms in the blanket. The most promising concept is that of a liquid metal blanket, in which a Li-containing fluid (most likely alloy PbLi) flows continuously through a network of ducts. The flow is affected by a very strong (up to 10–12 T) steady magnetic field and very strong internal heating due to absorption of neutrons (tens of
$\text{MW}~\text{m}^{-3}$
near the wall facing the reaction chamber).
The first of these effects has been a subject of rather extensive studies guided primarily by two considerations: that laminarization of the flow by the strong magnetic field results in poor mixing and transport properties and that the Lorentz force and magnetohydrodynamic (MHD) boundary layers lead to the flow resistance being much stronger than the resistance in analogous flows without magnetic fields (see, e.g., Molokov, Moreau & Moffatt Reference Molokov, Moreau and Moffatt2007).
Curiously, the effect of non-uniform internal heating on the flow was, for long time, largely ignored. This changed in recent laboratory experiments, such as those on mercury flows in pipes and ducts (see, e.g., Belyaev et al. Reference Belyaev, Listratov, Razuvanov and Sviridov2011, Reference Belyaev, Genin, Listratov, Melnikov, Sviridov, Sviridov, Ivochkin, Razuvanov and Shpansky2013; Genin et al. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011a ,Reference Genin, Zhilin, Ivochkin, Razuvanov, Sviridov, Shestakov and Sviridov b ; Melnikov et al. Reference Melnikov, Sviridov, Sviridov and Razuvanov2014; Poddubnyi et al. Reference Poddubnyi, Razuvanov, Sviridov and Ivochkin2014), and computations (Mistrangelo & Bühler Reference Mistrangelo and Bühler2013; Vetcha et al. Reference Vetcha, Smolentsev, Abdou and Moreau2013; Zikanov, Listratov & Sviridov Reference Zikanov, Listratov and Sviridov2013; Lv & Zikanov Reference Lv and Zikanov2014; Zhang & Zikanov Reference Zhang and Zikanov2014; Liu & Zikanov Reference Liu and Zikanov2015). While limited to moderately high values of Hartmann and Grashof numbers (experiments and accurate computations at the parameters of a real fusion reactor are as of yet unfeasible) and, therefore, not amounting to an unambiguous proof, the studies produced convincing data indicating that the structure and transport properties of blanket flows are changed profoundly by thermal convection. In many configurations, in particular in flows through vertical (Vetcha et al. Reference Vetcha, Smolentsev, Abdou and Moreau2013; Melnikov et al. Reference Melnikov, Sviridov, Sviridov and Razuvanov2014; Poddubnyi et al. Reference Poddubnyi, Razuvanov, Sviridov and Ivochkin2014; Liu & Zikanov Reference Liu and Zikanov2015) and horizontal (Genin et al. Reference Genin, Zhilin, Ivochkin, Razuvanov, Belyaev, Listratov and Sviridov2011a ; Lv & Zikanov Reference Lv and Zikanov2014; Zhang & Zikanov Reference Zhang and Zikanov2014) ducts and in boxes with conducting walls (Mistrangelo & Bühler Reference Mistrangelo and Bühler2013), convection leads to large-amplitude fluctuations of temperature or to the formation of hot and cold spots in the walls. Should similar effects appear in an actual blanket, they would result in strong and possibly unsteady thermal stresses in the walls causing rapid deterioration of wall material and even loss of the structural integrity of the blanket. In this sense, the convection effect is highly undesirable and may even invalidate the currently pursued blanket designs, such as dual coolant lead lithium (DCLL; Smolentsev, Moreau & Abdou Reference Smolentsev, Moreau and Abdou2008) or helium cooled lead lithium (HCLL; Mistrangelo & Bühler Reference Mistrangelo and Bühler2009).
In this paper, we argue and illustrate by numerical simulations that the effect of convection is not necessarily detrimental, but can, in fact, be exploited as an efficient mechanism of mixing and heat transfer within a blanket. A flow in a long duct aligned with a strong uniform magnetic field is analysed. The system can be considered as an idealized model of a toroidal (parallel to the main component of the magnetic field) duct of a blanket. Such designs were considered in the early days of fusion technology (Hunt & Hancox Reference Hunt and Hancox1971; Smith et al. Reference Smith, Baker, Sze, Morgan, Abdou, Piet, Schultz, Moir and Gordon1985) for self-cooled (both breeding and cooling carried out by liquid metal) blankets. The concept was, probably prematurely, abandoned in the 1990s, after theoretical and numerical analysis of laminar steady-state flows showed complexity and poor predictability related to strong sensitivity to variations of the magnetic field and wall material properties, details of the design, etc. (see, for example, the chapter by Bühler in Molokov et al. Reference Molokov, Moreau and Moffatt2007).
The focus of this work is on the theoretical possibility of a blanket with toroidal separately cooled (heat transfer is mostly carried out by external heat exchangers) ducts. Apart from their importance for fusion technology, the results are, in our view, of general interest and beauty. A rare example of two-dimensional turbulence in a realistic system is presented.
2. Physical model and numerical method
A flow of an incompressible electrically conducting fluid (a liquid metal) in a horizontal duct of square cross-section is considered (see figure 1
a). There is non-uniform internal heating of volumetric rate
$q_{0}q(x)$
, while the walls of the duct are maintained at constant temperature
$T_{0}$
. There is no mean flow along the duct. The system can be considered as a model of a blanket, where cooling is carried out primarily by an auxiliary (e.g. pressurized He) circuit built into walls and the flow of liquid metal along the duct is only needed for tritium extraction and, therefore, is negligibly slow (
$1~\text{mm}~\text{s}^{-1}$
or less).
Figure 1. (a) Geometry of the flow and coordinate system. For the three-dimensional flows considered in § 3.2, the sketch represents a cross-section of the flow domain that has the form of a long duct. (b) The non-dimensional internal heating rate
$q(x)$
(solid line) and the temperature distribution (dashed line) that would develop in response to this heating in a purely conductive regime without the influence of the walls at
$y=\pm 1$
.
As an approximation, we neglect the poloidal component of the magnetic field (approximately 5 % of the total field strength in actual reactors) and curvature effects. The magnetic field is assumed to be uniform, horizontal and perfectly aligned with the axis of the duct. We also assume that the duct is long (tens of hydraulic diameters), so that the effects of its ends can be neglected and an asymptotic model with periodic inlet–exit conditions can be considered in three-dimensional simulations. Furthermore, in such a system, a sufficiently strong magnetic field completely suppresses axial velocity gradients and makes the flow two-dimensional (see, e.g., Davidson Reference Davidson2001). The electric currents and Lorentz force become zero in such a flow. The two-dimensional state is stable or unstable to three-dimensional perturbations depending on the strength of the magnetic field, as measured by the Stuart number,
$N=Ha^{2}Re^{-1}$
, where
$Ha$
and
$Re$
are properly defined Hartmann and Reynolds numbers (see the discussion below for the definition relevant to our case), and on the axial wavelength of the perturbations (see, e.g., Thess & Zikanov Reference Thess and Zikanov2007). The Stuart number corresponding to the typical blanket condition is of the order of
$10^{2}$
, which indicates stability. We assume stability at first and perform tests of the robustness of this assumption later.
As typical scales, we use the duct half-width
$d$
for length, typical magnitude
$q_{0}$
for internal heating rate,
${\rm\Delta}T=q_{0}d^{2}{\it\kappa}^{-1}$
for temperature, free-fall speed
$U=\sqrt{{\it\beta}g{\rm\Delta}Td}$
for velocity,
$d/U$
for time and
${\it\rho}U^{2}$
for pressure, and obtain, with the Boussinesq approximation,
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+\frac{1}{\sqrt{Gr}}{\rm\nabla}^{2}\boldsymbol{u}+T\boldsymbol{e}_{y}+\boldsymbol{F}_{L}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T=\frac{1}{Pr\sqrt{Gr}}({\rm\nabla}^{2}T+q), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}=T=0\quad \text{at}~x=\pm 1,y=\pm 1, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}=(u,v,w)$
and
$T$
are the non-dimensional velocity and temperature deviation from
$T_{0}$
, and 
$$\begin{eqnarray}q=\exp -(x+1)\end{eqnarray}$$
is an approximation of the non-dimensional rate of internal heating due to absorption of neutrons proposed by Smolentsev, Morley & Abdou (Reference Smolentsev, Morley and Abdou2006) and shown in figure 1(b) (the wall at
$x=-1$
is nearest to the reaction chamber). The Lorentz force
$\boldsymbol{F}_{L}$
is zero in two-dimensional (
$z$
-independent) flow. In three-dimensional flows, it is determined as
$$\begin{eqnarray}\boldsymbol{F}_{L}=Ha^{2}Re^{-1}\boldsymbol{j}\times \boldsymbol{e}_{z}.\end{eqnarray}$$
The electric current
$\boldsymbol{j}$
is determined by Ohm’s law,
$$\begin{eqnarray}\boldsymbol{j}=-\boldsymbol{{\rm\nabla}}{\it\phi}+\boldsymbol{u}\times \boldsymbol{e}_{z},\end{eqnarray}$$
where the electric potential
${\it\phi}$
is a solution of the Poisson equation expressing the instantaneous electric neutrality of the fluid assumed in the quasistatic approximation:
$$\begin{eqnarray}{\rm\nabla}^{2}{\it\phi}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{u}\times \boldsymbol{e}_{z}).\end{eqnarray}$$
The boundary conditions for
${\it\phi}$
are those of perfect electric insulation,
$$\begin{eqnarray}\frac{\partial {\it\phi}}{\partial n}=0\quad \text{at}~x=\pm 1,y=\pm 1.\end{eqnarray}$$
Periodicity of the electric potential
${\it\phi}$
, velocity
$\boldsymbol{u}$
, temperature fluctuations
$T$
and pressure fluctuations
$p$
at
$z=0,L$
, where
$L$
is the length of the computation domain, is assumed in the case of three-dimensional flows. The non-dimensional parameters are the Prandtl number
$Pr={\it\nu}/{\it\chi}$
, the Grashof number
$$\begin{eqnarray}Gr=\frac{g{\it\beta}q_{0}d^{5}}{{\it\kappa}{\it\nu}^{2}},\end{eqnarray}$$
which is a square of the Reynolds number
$Re=Ud/{\it\nu}$
, and the Hartmann number
$$\begin{eqnarray}Ha=Bd\left(\frac{{\it\sigma}}{{\it\rho}{\it\nu}}\right)^{1/2}.\end{eqnarray}$$
In the formulae above,
${\it\kappa}$
,
${\it\beta}$
,
${\it\nu}$
,
${\it\chi}$
and
${\it\sigma}$
are respectively the thermal conductivity, thermal expansion coefficient, kinematic viscosity, temperature diffusivity and electric conductivity of the fluid. Applicability of the Boussinesq approximation to convection flows within fusion reactor blankets is, generally, questionable because temperature variations can be large (see Gray & Giorgini Reference Gray and Giorgini1976 for a discussion of application limits). Nevertheless, the Boussinesq approximation is typically used in blanket studies (see, e.g., Vetcha, Smolentsev & Abdou Reference Vetcha, Smolentsev and Abdou2009; Mistrangelo & Bühler Reference Mistrangelo and Bühler2011; Belyaev et al.
Reference Belyaev, Genin, Listratov, Melnikov, Sviridov, Sviridov, Ivochkin, Razuvanov and Shpansky2013; Lv & Zikanov Reference Lv and Zikanov2014; Zhang & Zikanov Reference Zhang and Zikanov2014). As indicated by Ni et al. (Reference Ni, Qiu, Su, Tian and Wu2012), the non-Boussinesq effects can be significant quantitatively, but are unlikely to have profound qualitative effects on the flow. Our approach in this study is to use the approximation, partially justifying it by the fact that the particularly interesting flows for us are turbulent and, thus, temperature non-uniformity is reduced, and postpone considering its accuracy for the future.
The problem is solved numerically in primitive variables using the finite-difference scheme developed by Krasnov, Zikanov & Boeck (Reference Krasnov, Zikanov and Boeck2011) and modified to include thermal convection and implicit treatment of viscosity and conduction terms by Zikanov et al. (Reference Zikanov, Listratov and Sviridov2013) and Zhang & Zikanov (Reference Zhang and Zikanov2014). The spatial and time discretizations are of the second order. The grid is uniform in the axial direction (in three-dimensional computations) and clustered towards walls according to the coordinate transformation
$$\begin{eqnarray}x=A_{x}\sin ({\rm\pi}{\it\xi}/2)+(1-A_{x}){\it\xi},\quad y=A_{y}\sin ({\rm\pi}{\it\eta}/2)+(1-A_{y}){\it\eta},\end{eqnarray}$$
where
$-1\leqslant {\it\xi}\leqslant 1$
and
$-1\leqslant {\it\eta}\leqslant 1$
are the transformed coordinates in which the grid is uniform and the blending coefficients
$A_{x}=A_{y}=0.96$
are used. The same combination of clustering and finite-difference discretization was successfully used by Krasnov, Zikanov & Boeck (Reference Krasnov, Zikanov and Boeck2012) for simulations of MHD duct flows at high Reynolds and Hartmann numbers. Each simulation is conducted for a long time to obtain a fully developed state and then accumulate data over an interval sufficient for accurate time averaging. The Nusselt number is defined as
$$\begin{eqnarray}Nu=\frac{Q}{\overline{T}},\end{eqnarray}$$
where
$Q=\int _{A}q\,\text{d}A$
is the total heating rate per unit length of the duct and
$\overline{T}=A^{-1}\int _{A}T\,\text{d}A$
is the mean temperature. In addition to that, the following integral characteristics are computed in the two-dimensional flows (see, e.g., Clercx & van Heijst Reference Clercx and van Heijst2009, for a discussion of their meaning): kinetic energy, enstrophy and angular momentum with respect to the centre of the duct,
$$\begin{eqnarray}\displaystyle & \displaystyle E=\frac{1}{A}\int _{A}|\boldsymbol{u}|^{2}\,\text{d}A, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\Omega}=\frac{1}{2A}\int _{A}{\it\omega}^{2}\,\text{d}A, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle L=\frac{1}{A}\int _{A}(xv-yu)\,\text{d}A=-\frac{1}{2A}\int _{A}r^{2}{\it\omega}\,\text{d}A. & \displaystyle\end{eqnarray}$$
$W=(1/A)\int _{A}{\it\omega}\,\text{d}A$
, where
${\it\omega}=\partial v/\partial x-\partial u/\partial y$
. It is easy to verify by direct integration of
${\it\omega}$
and use of velocity boundary conditions that
$W$
should be zero for all flows.
In the simulations of three-dimensional flows, our main attention is given to the perturbations with respect to the two-dimensional solutions:
$$\begin{eqnarray}\boldsymbol{u}^{\prime }=\boldsymbol{u}-\overline{\boldsymbol{u}},\quad T^{\prime }=T-\overline{T},\end{eqnarray}$$
where
$\overline{\boldsymbol{u}}$
and
$\overline{T}$
are the instantaneous averages in the axial direction. In particular, we follow the deviation from two-dimensionality using the perturbation energies
$$\begin{eqnarray}E^{\prime }=\langle \,f^{\prime 2}\rangle ,\end{eqnarray}$$
where
$\langle \cdots \rangle$
stands for volume averaging and
$f$
stands for perturbations of a velocity component or temperature.
A grid sensitivity study has been performed to identify, for every value of
$Gr$
, the size of the grid such that further increase does not significantly change the integral characteristics of the two-dimensional flow (see table 1). We also require that the average vorticity
$W$
obtained by the finite-difference differentiation of velocity components and the Simpson-rule integration of
${\it\omega}$
on the non-uniform grid is sufficiently close to zero. Table 1 shows that the error is small in comparison to the r.m.s. vorticity, which is approximately one in all the cases. In three-dimensional analysis, the grid step in the axial direction
${\rm\Delta}z$
is always kept as approximately 0.1, which is 100 times smaller than the typical wavelength of the unstable modes observed.
Table 1. Grid sensitivity study: the parameters of the computational grids and the integral characteristics of the two-dimensional solutions are shown. The characteristics are obtained by time averaging at
$Gr\geqslant 10^{7}$
. The grids used in the simulations of two-dimensional flows are underlined.
Results are presented for
$10^{6}\leqslant Gr\leqslant 10^{11}$
(see table 1),
$Pr=0.0321$
(PbLi alloy at approximately 570 K) and
$800\leqslant Ha\leqslant 10^{4}$
. This can be compared with the typical values expected in full-scale fusion reactors:
$Gr$
up to
$10^{12}$
and
$Ha$
up to
$10^{4}$
(see Molokov et al.
Reference Molokov, Moreau and Moffatt2007). The axial length of the computational domain in the three-dimensional simulation is
$4{\rm\pi}\leqslant L_{z}\leqslant 30{\rm\pi}$
.
3. Results
3.1. Two-dimensional flows
The non-uniform temperature distribution produced by internal heating causes thermal convection flows in all of our cases. The structure and behaviour of the developed flows are illustrated in figures 2–4. We show instantaneous distributions of the streamfunction, temperature and amplitude of the velocity in figure 2, time evolution of the Nusselt number and average kinetic energy in figure 3 and spectral decompositions of the average kinetic energy signal
$E(t)$
in figure 4. At
$Gr=10^{6}$
, the flow is steady state and consists of two circulation rolls (see figure 2
a–c). Topologically similar but unsteady flows are obtained at
$Gr=10^{7}$
and
$10^{8}$
(see figure 2
d–f). At
$Gr=10^{7}$
, the unsteady fluctuations have three isolated frequencies and low amplitudes in comparison to the mean (see figure 4
a). At
$Gr=10^{8}$
(see figures 2
d–f and 4
b), the amplitude is higher, and the Fourier analysis of the signal gives one strong frequency and weaker frequencies distributed continuously over a small interval.
Figure 2. Instantaneous distributions of the streamfunction
${\it\Psi}$
(solid lines indicate counterclockwise motion while dashed lines indicate clockwise motion), temperature
$T$
and amplitude of the velocity
$|\boldsymbol{u}|$
in two-dimensional flows at
$Gr=10^{6}$
(a–c),
$Gr=10^{8}$
(d–f) and
$Gr=10^{11}$
(g–i). Note that the isolevels are different at different values of
$Gr$
.
Figure 3. The Nusselt number (a) and average kinetic energy (b) in two-dimensional flows at
$Gr=10^{8}$
(dashed lines) and
$Gr=10^{11}$
(solid lines). Only parts of actual simulations are shown. Note that different scales are used at
$Gr=10^{8}$
and
$Gr=10^{11}$
.
Figure 4. Spectral decompositions of the average kinetic energy signal
$E(t)$
in two-dimensional flows at
$Gr=10^{7}$
(a),
$Gr=10^{8}$
(b) and
$Gr=10^{11}$
(c) (in log scale).
In the most interesting cases for us of high
$Gr$
(at
$10^{9}\leqslant Gr\leqslant 10^{11}$
), the flow is turbulent (see figure 2
g–i and the curves for
$Gr=10^{11}$
in figures 3 and 4). The Fourier analysis produces continuous spectra in wide ranges of frequencies (see figure 4
c). The flow fields show features typical for forced two-dimensional turbulence in a confined domain (see, e.g., Clercx & van Heijst Reference Clercx and van Heijst2009): continuous instability and breakdown of large-scale structures, their recreation by forcing, thin shear layers and vorticity filaments at the walls. We note that an inverse energy cascade is not to be expected in our system, since the convection forcing acts on a length scale of the same order as the size of the flow domain.
The integral properties (2.13)–(2.16) are listed in table 1. We see that the average enstrophy
${\it\Omega}$
changes little with
$Gr$
. With the kinetic energy balance in a statistically steady state being
$-2Gr^{-1/2}{\it\Omega}+(1/A)\int _{A}Tv\,\text{d}A\approx 0$
, this implies that the energy production by convection forcing satisfies
$$\begin{eqnarray}\int _{A}Tv\,\text{d}A\sim Gr^{-1/2}.\end{eqnarray}$$
As illustrated in figure 5, the other integral properties change significantly with
$Gr$
. The Nusselt number grows as
$$\begin{eqnarray}Nu\sim Gr^{0.2}\end{eqnarray}$$
in turbulent regimes. Similar scaling was recently observed by Goluskin & Spiegel (Reference Goluskin and Spiegel2012) in simulations of two-dimensional convection in a horizontal layer with uniform internal heating. The average kinetic energy
$E$
and the amplitude of the average angular momentum
$-L$
decrease by orders of magnitude apparently as a result of redistribution of kinetic energy from large to small length scales. Interestingly, the scaling coefficients at large
$Gr$
are nearly the same for the two quantities and are equal to the negative of the scaling coefficient for the Nusselt number. In all of our flows,
$L<0$
, which reflects the asymmetry introduced by the forcing.
Figure 5. Integral time-averaged (at
$Gr\geqslant 10^{7}$
) characteristics of the two-dimensional flows listed in table 1. Approximate slopes at
$10^{9}\leqslant Gr\leqslant 10^{11}$
are indicated.
For the design of a liquid metal blanket with toroidal ducts imitated by our model (with very low flow rate and cooling by heat exchangers built into the walls), it is important to know how the heat flux from the interior is distributed over the walls. The non-dimensional heat flux
$$\begin{eqnarray}q_{w}=\left.-\frac{\text{d}T}{\text{d}n}\right|_{wall},\end{eqnarray}$$
where
$\boldsymbol{n}$
is an outward-facing normal, is shown in figure 6 for the three flow states of figure 2. We do not see potentially dangerous localized zones of particularly strong or particularly weak flux. At
$Gr=10^{6}$
and
$10^{8}$
, the walls at
$x=-1$
and
$y=1$
are heated somewhat more strongly than the other two walls. In the turbulent flow at
$Gr=10^{11}$
, intense mixing results in a more uniform heat flux distribution.
Figure 6. Distributions of the time-averaged wall heat flux corresponding to the two-dimensional flow states in figure 2. The averaging is performed at the stage of fully developed flow over 200 time units. At
$Gr=10^{11}$
, this corresponds to approximately eight convective time units
$(E/2)^{-1/2}$
. Note that different scales are used for
$q_{w}$
at
$Gr=10^{6}$
(a),
$Gr=10^{8}$
(b) and
$Gr=10^{11}$
(c).
3.2. Stability of two-dimensional flow regimes to three-dimensional perturbations
Three-dimensional transient analysis is performed to test the robustness of the two-dimensionality assumption. We take a two-dimensional flow as an initial state, add small-amplitude three-dimensional random perturbations and compute the flow evolution for at least 200 time units. The two-dimensional state is considered to be stable or unstable if the perturbation energies (2.18) behave as illustrated in figure 7(a) or (b) respectively, i.e. show continuous nearly exponential decay or grow and then saturate at a higher amplitude. The flows are computed in a domain with periodic inlet/exit conditions, length
$L_{z}=4{\rm\pi}$
,
$10{\rm\pi}$
or
$30{\rm\pi}$
and at
$800\leqslant Ha\leqslant 10^{4}$
. The domain length
$L_{z}$
is an important parameter of the analysis. It determines the maximum of the axial wavelength
${\it\lambda}_{z}$
of the perturbations allowed by the model. This affects the strength of suppression of perturbations by the magnetic field, since the rate of Joule dissipation of the electric currents induced by a flow structure is proportional to
$Ha^{2}Re^{-1}{\it\lambda}_{z}^{-2}$
. At any strength of the magnetic field one can, therefore, select a sufficiently long domain such that perturbations with
${\it\lambda}_{z}\sim L_{z}$
survive the suppression. Our choice
$L_{z}\leqslant 30{\rm\pi}$
can be related to existing designs of liquid metal blankets for fusion reactors, where duct widths tend to be of the order of a few cm and duct lengths rarely exceed 2 m. The non-dimensional length
$L_{z}$
is, thus, not larger than 100. The parameters of the three-dimensional simulations conducted are listed in table 2. The most extensive computations have been conducted at
$Gr=10^{9}$
, while at
$Gr=10^{10}$
and
$Gr=10^{11}$
just a few runs have been performed to confirm that the situation is qualitatively the same. In order to make three-dimensional computations feasible, the size of the computational grid in the
$x$
–
$y$
-plane has been reduced in comparison to the two-dimensional simulations (see tables 1 and 2).
Figure 7. The evolution of the three-dimensional perturbation energy
$E^{\prime }$
is shown for
$Ha=800$
,
$Gr=10^{9}$
,
$L_{z}=4{\rm\pi}$
(a) and
$Ha=2000$
,
$Gr=10^{9}$
,
$L_{z}=4{\rm\pi}$
(b).
Table 2. Results of the three-dimensional simulations. The computational parameters, the energy of the three-dimensional perturbations
$E^{\prime }$
and the ratio between
$E^{\prime }$
and the kinetic energy of the corresponding two-dimensional flows
$E_{2D}$
are shown.
The main results are illustrated by the case
$Gr=10^{9}$
in figure 8. For every
$L_{z}$
, there is a critical
$Ha$
above which the perturbation energies decrease and the flow maintains two-dimensionality. When
$Ha$
is lower than the critical value, transition to three-dimensionality occurs. Calculations conducted at
$Gr=10^{10}$
and
$10^{11}$
produce qualitatively similar results and show that the two-dimensional states become less stable at higher
$Gr$
– the stability requires larger
$Ha$
or smaller
$L_{z}$
. An example of a developed three-dimensional flow is illustrated in figure 9 by the case
$Ha=10^{4}$
,
$Gr=10^{11}$
,
$L_{z}=10{\rm\pi}$
. The Hartmann number in this case is smaller than the critical value, and the perturbation energies increase in a manner similar to figure 7(b). The magnitude of the perturbation energies at the saturation level is small, of the order of
$10^{-5}$
. We can estimate that the perturbation velocity is only approximately 2 % of the velocity of the two-dimensional flow (see table 2). This indicates that the flow at
$Ha=10^{4}$
,
$Gr=10^{11}$
and
$L_{z}=10{\rm\pi}$
remains nearly two-dimensional. This is confirmed by the two-point correlation coefficient shown in figure 9(e). The coefficient is calculated using the vertical velocity component as
$$\begin{eqnarray}R(L)=\frac{\langle v(\boldsymbol{x},t)v(\boldsymbol{x}+L\boldsymbol{e}_{z},t)\rangle }{\langle v^{2}(\boldsymbol{x},t)\rangle },\end{eqnarray}$$
where
$\boldsymbol{x}$
is the point
$(0,0,z)$
and the averaging is performed over the streamwise direction. No time averaging is applied in (3.4). We see that, although there are some three-dimensional structures, the flow does not vary much in the axial direction, with the flow coefficient remaining above 0.99. Apart from weak three-dimensional perturbations, the flow remains very similar both qualitatively and quantitatively to the flow obtained in the two-dimensional computations (compare figures 9
d and 2
g–i). It can also be observed in figure 9 that the three-dimensional structures are located near the bottom of the duct and have a typical axial wavelength of approximately 12 unit lengths. The bottom location suggests that the transition to three-dimensionality in this flow may be caused by instabilities forming in the thin shear layers near the duct bottom in the two-dimensional solution (see figure 9
d).
Figure 8. Stability diagram for the two-dimensional solutions at
$Gr=10^{9}$
. Crosses indicate the stability thresholds determined in the computations.
Figure 9. Results of the three-dimensional simulations. Instantaneous distributions of the vertical velocity
$u_{y}$
, temperature
$T$
and spanwise velocity
$u_{x}$
at
$Ha=10^{4}$
,
$Gr=10^{11}$
,
$L_{z}=10{\rm\pi}$
are shown in the cross-section
$x=0$
(a–c), and of
$T$
and the velocity vectors in the transverse cross-section
$z=L_{z}/2$
(d). The two-point correlation coefficient
$R$
is shown in (e). The scales of the
$y$
and
$z$
coordinates in (a–c) are related as 1:3.
The strength of the three-dimensional structures is largely dependent on how far away the parameters are from the critical line. For example, at
$Ha=10^{4}$
,
$Gr=10^{11}$
and
$L_{z}=4{\rm\pi}$
, which can be compared with the
$L_{z}=10{\rm\pi}$
in figure 9, we still observe a transition to three-dimensionality, but visual inspection of the plots analogous to figure 9(a–c) yields only a slight tilt of the isolines. The energies of the three-dimensional velocity perturbations saturate at approximately
$10^{-6}$
, i.e. at a level significantly lower than for
$L_{z}=10{\rm\pi}$
(see table 2).
On the contrary, at parameters far from the critical line, the three-dimensionality becomes clearly visible. This is illustrated in figure 10 by the case
$Gr=10^{10}$
,
$Ha=2000$
,
$L_{z}=10{\rm\pi}$
. The correlation coefficient drops to approximately 0.5 in this case. The three-dimensional flow structures are clearly visible, have a range of active wavelengths and occupy the entire duct.
Figure 10. Results of the three-dimensional simulations. Instantaneous distributions of the vertical velocity
$u_{y}$
, temperature
$T$
and spanwise velocity
$u_{x}$
at
$Ha=2000$
,
$Gr=10^{10}$
,
$L_{z}=10{\rm\pi}$
are shown in the cross-section
$x=0$
(a–c), and of
$T$
and the velocity vectors in the transverse cross-section
$z=L_{z}/2$
(d). The two-point correlation coefficient
$R$
is shown in (e). The scales of the
$y$
and
$z$
coordinates in (a–c) are related as 1:3.
We have explored the effect of three-dimensionality on heat transfer and found that the Nusselt number is increased as a result of the transition. As an illustration, figure 11 shows
$Nu$
in the three-dimensional flow at
$Gr=10^{9}$
,
$Ha=800$
and
$L_{z}=4{\rm\pi}$
approximately 8 % higher than in its two-dimensional counterpart.
Figure 11. The Nusselt number
$Nu$
at
$Gr=10^{9}$
as computed for the two-dimensional regime (black, solid) and the three-dimensional regime at
$Ha=800$
,
$L_{z}=4{\rm\pi}$
(red, dashed).
One can see a certain analogy between our flow and the inertial (flywheel) convection found in low-Prandtl Rayleigh–Bénard or Marangoni systems, often with periodic or free-slip boundaries (see, e.g., Jones, Moore & Weiss Reference Jones, Moore and Weiss1976; Clever & Busse Reference Clever and Busse1981; Boeck & Thess Reference Boeck and Thess1997). The inertial convection has the form of strong and nearly circular two-dimensional rolls with nearly uniform vorticity and weak viscous dissipation in the core. Considering the structure of our flows in figure 2, we see that the analogy is, to some degree, valid at low
$Gr$
. At high
$Gr$
, the vorticity filaments generated at the solid walls penetrate into the core flow leading to breakdown of large rolls and two-dimensional turbulence. The analogy becomes invalid.
Our last point of concern is the phenomenon of large-scale intermittency discovered in the MHD flow in a channel with a spanwise magnetic field by Boeck et al. (Reference Boeck, Krasnov, Thess and Zikanov2008) and further explored by Dey & Zikanov (Reference Dey and Zikanov2012). It has been found for this flow that in a substantial range of
$Ha$
, the magnetic field is not strong enough to prevent three-dimensional secondary instability of the spanwise Tollmien–Schlichting rolls, but is sufficiently strong to suppress the turbulence after it develops. As a result, the flow experiences nearly periodic evolution through cycles of growth of two-dimensional rolls, their three-dimensional instability leading to a turbulent burst and decay to nearly a laminar Poiseuille flow.
There appears to be no evident reason why a similar intermittent behaviour cannot exist in our case or in the case of inertial convection, the Tollmien–Schlichting rolls being replaced by the two-dimensional convection structures described in § 3.1. Our computations at
$Gr=10^{9}$
,
$Gr=10^{10}$
and
$Gr=10^{11}$
conducted near the three-dimensional instability thresholds for long times (up to 3000 time units) have not produced a clear intermittent behaviour. In all runs at supercritical parameters, the flows remain three-dimensional indefinitely after the initial growth and saturation of perturbations. A possible explanation is that, unlike the channel flow case of Boeck et al. (Reference Boeck, Krasnov, Thess and Zikanov2008) and Dey & Zikanov (Reference Dey and Zikanov2012), the three-dimensional perturbations remain weak in our case and do not destroy the underlying two-dimensional flow structures, and thus do not remove the mechanisms leading to the three-dimensional instability. At the same time, we have observed fluctuations of
$E^{\prime }$
by nearly an order of magnitude in several simulations at parameters very close to the threshold (for example at
$Gr=10^{9}$
,
$Ha=2500$
,
$L_{z}=20{\rm\pi}$
and
$Gr=10^{9}$
,
$Ha=3700$
,
$L_{z}=30{\rm\pi}$
). Furthermore, intermittent behaviour may occur at lower
$Gr$
. These questions deserve future consideration.
4. Concluding remarks
We have considered convection caused by internal heating in a horizontal duct with a longitudinal magnetic field. The flow can be considered as a model of liquid metal flows in a conceptual blanket for a tokamak fusion reactor. In such a blanket, liquid metal moves slowly through toroidal ducts. Cooling and transport of heat to an external power generation facility is carried out by a conventional (e.g. water or compressed helium) heat exchanger built into the walls, while the liquid metal serves the breeding and shielding purposes.
It has been first assumed in our analysis that the magnetic field is sufficiently strong to make the flow two-dimensional and reduce the problem to that of natural convection in a square.
Using, as an example, the duct half-width
$d=5$
cm and physical properties of LiPb at 570 K (Schulz Reference Schulz1991), we can evaluate the combination
$g{\it\beta}d^{5}/{\it\kappa}{\it\nu}^{2}$
as approximately
$536~\text{m}^{3}~\text{W}^{-1}$
. This means that the range
$10^{9}\leqslant Gr\leqslant 10^{11}$
, in which we have found two-dimensional turbulence, corresponds to a typical heating rate
$q_{0}$
from approximately 1.9 to
$190~\text{MW}~\text{m}^{-3}$
. This range entirely covers the variety of heating rates expected in both experimental (such as ITER) and future production fusion reactors.
Three-dimensional computations have shown that the flow maintains its two-dimensional form at
$Ha\sim 10^{4}$
and
$10^{9}\leqslant Gr\leqslant 10^{11}$
typical for reactor conditions in ducts of fairly large length. Even when transition to three-dimensionality occurs, the energy of the three-dimensional perturbations is small in comparison to the energy of the two-dimensional flow components and is limited to large wavelengths (see figures 7
b and 9 for an example of such a flow at typical values of
$Gr$
and
$Ha$
). The three-dimensionality has little impact on the flow structure and little to moderate impact on the heat transfer, so the key flow parameters can be reasonably accurately evaluated in a two-dimensional analysis. The situation may change when realistic inlet and exit conditions and accompanying three-dimensional MHD effects are included. This deserves to be addressed in future studies.
Accepting that the flow of liquid metal in the ducts of our hypothetical blanket is turbulent and has structure and properties similar to those obtained for
$Gr=10^{11}$
in this work, we immediately arrive at the conclusion that the concept is promising. Turbulent two-dimensional convection provides an effective mechanism of mixing and transport resulting in a reasonably uniform heat transfer into cooled walls. Unlike the currently pursued blanket designs (Smolentsev et al.
Reference Smolentsev, Moreau and Abdou2008; Mistrangelo & Bühler Reference Mistrangelo and Bühler2009), the effect of convection is not detrimental, but is, in fact, exploited as a beneficial mechanism of the blanket’s operation.
Although the results are theoretical at this point, they show potential new advantages of the old concept of a blanket with toroidal ducts. This inspires us to conclude that further studies are necessary, in which one would consider real-life complexities, such as the finite length of the duct and deviations from two-dimensionality, as well as features of a specific design of the blanket.
Acknowledgements
The authors are grateful to Y. Listratov for a careful reading of the manuscript and valuable suggestions. Financial support was provided by the US NSF (Grant CBET 1232851).
