1. Introduction
Breeding blankets are crucial systems projected in future nuclear fusion reactors by magnetic confinement whose main purpose is to regenerate the tritium burnt in the plasma. They are designed to absorb the energetic neutrons created in the fusion reactions (). In the core of the blanket, the neutrons produce tritium by transmuting the lithium contained in the blanket. The amount of tritium bred this way has to be extracted from the reactor for being processed, stored and eventually injected in the plasma. When including a neutron multiplier in the blanket, the tritium bred in the blanket can be higher than the tritium burnt in the plasma which allows maintaining the desired tritium self-sufficiency of the plant.
Within the framework of the EUROfusion blanket project [
1,
2], the WCLL blanket concept [
3] has been selected as one of the two candidates for driver blanket of the European DEMO. This concept is based on an eutectic alloy of PbLi (Pb acting as neutron multiplier and Li acting as breeder) flowing at low velocities. The blanket is directly exposed to highly energetic particle fluxes coming from the blanket. For this reason, liquid water pressurized at 155 bar is used for cooling both the PbLi and the steel structures. A network of water tubes immersed in the liquid metal is included in the design for this purpose.
The PbLi is a good electrical conductor immersed in the strong magnetic field used for the plasma confinement. Therefore, electric currents will be induced inside the flow producing Lorentz forces that alter its dynamics. This interaction known as Magnetohydrodynamics (MHD) is dominant in most of the WCLL PbLi flowpath. Nevertheless, in the regions closer to the first wall (FW), the high particle fluxes from the plasma are very energetic which are translated into a very high non-homogeneous heat source. This together with the arrangement of the cold water tubes immersed in the PbLi flow can produce temperature differences inside the bulk of the fluid of hundreds of degrees. In these areas, the buoyancy forces are expected to be comparable or even higher than the Lorentz ones.
In this work, magneto-convective simulations of the WCLL blanket are presented. Those are focused on a very particular region of the design: the frontal part of the central outboard elementary cell (
Figure 1). This region is located very close to the plasma and it is characterized by a very high heat flux and a consequent dense arrangement of water tubes. This location is of the special interest from the blanket design perspective since it presents a very high tritium production. In its vicinity, it is expected that buoyancy forces play a very important role on PbLi dynamics.
The WCLL elementary cell is composed of six parallel circuits along the toroidal direction. Each of them is fed by a rear PbLi manifold and presents a U-shape flowpath along the radial direction. The analyses of this work are focused on the frontal part of one of the six parallel circuits; in the area where the flow describes a 180
turn, close to the FW of the blanket. It is worth noting that there is not a physical separation in between the frontal regions of the six parallel circuits. Therefore, in a central circuit, such as the one considered in this analysis, the lateral radial-poloidal walls only cover the rear part of the domain (the so-called radial channels
Figure 2).
The numerically analyzed geometry has been simplified. In particular, the shape of the third row of tubes has been straightened, making them completely toroidaly oriented. This way, the studied geometry presents a good symmetry along the toroidal direction which allows using structured meshes in the analyses. Indeed, structured meshes provide better stability in MHD simulations. In the ANSYS-Fluent MHD solver, unstructured meshes can lead to violations of charge conservation forcing extremely small time steps which would rise the computational cost.
The simplification performed does not affect the first and second rows of tubes which have the same shape and position than in the WCLL design. These tubes are closer to the FW where buoyancy forces are expected to be stronger. The square section FW cooling channels are not directly included in the computational domain. It only covers the part of the FW in between the cooling channels and the PbLi.
2. Numerical Model
In the present analyses, the dynamics of the PbLi flow is ruled by the Navier–Stokes equation in which the Lorentz and buoyancy terms are added as a volumetric sources:
where the PbLi density and dynamic viscosity are denoted by
and
, respectively. The vector field
is the velocity of the fluid. The Lorentz force is the cross product of the current density (
) and the external magnetic field (
) which is assumed to be constant.
represents the gravity vector field and
T the temperature field.
Under the inductionless approximation (the magnetic field created by the induced currents is negligible in comparison with the external one), the electric potential (
) is well defined and the generalized Ohm’s law can be written as follows:
where
is the electrical conductivity of the media. The MHD problem can be closed by applying the divergence to (
2) and assuming charge conservation (
). Therefore, a Poisson’s equation for the electric potential has to be solved coupled with (
1). This equation is solved also in the solid domains. The electric potential and electric currents are continuous across the solid-fluid interfaces.
Concerning the buoyancy forces in the fluid, they are treated using the Boussinesq approximation:
where
is the characteristic temperature of the problem,
is the PbLi density at that temperature and
is the thermal expansion coefficient. The rest of the fluid and solid properties are assumed to be constant.
The temperature field evolves following the energy conservation equation:
where
is the specific heat capacity of the material and
its thermal conductivity. The volumetric heat source
Q represents the effect of the neutrons and photons (nuclear heating). As a first approximation, this source only depends on the radial coordinate. The shape of this kind of functions are derived from the neutronics analyses performed for the WCLL (e.g., [
4]). In this work, the following piecewise functions have been used. These functions provide a steeper power density in the region close to the FW in comparison with the exponential functions used in previous studies which are only accurate in the tail of the curve (e.g., [
5,
6]).
The origin of the radial coordinate is located in the contact plane between the FW and the PbLi. The
r coordinate decreases towards the PbLi inlet and grows towards the FW. The shape of the volumetric heat source is plotted in
Figure 3.
The magnitude of the nuclear heating is of the order of MW/m. The Ohmic or Joule heating (∼) is expected to be several orders of magnitude smaller than the nuclear heating. Indeed, in the WCLL conditions: mm/s, and m), which means that: A/m. Thus, the Joule heating (∼ W/m) has been neglected in the present analysis.
In a dimensionless analysis of (
1), the square of the Hartmann number (
) represents the ratio between Lorentz forces and Viscous forces:
where
a is the channel semilength along the field direction (toroidal direction). In the central outboard WCLL conditions, the Hartmann number in the frontal part of the blanket is estimated to be approximately 8800. This means that viscous forces will be strongly suppressed by the magnetic field.
Complementary, the Grashof number (
) represents the ratio between the buoyancy and viscous forces:
where
b is the characteristic length of the heat transfer process,
is the characteristic temperature difference and
g is the gravitational acceleration. Estimating
in the different regions of the WCLL is not a straightforward task. Indeed, the heat source is extended along the complete computational domain (in a non-homogeneous way) and there are multiple cooling regions. Therefore, both the characteristic length and the temperature difference are not known before performing an specific thermal calculation. It is sometimes found in the literature (e.g., [
7,
8]) an estimation of the characteristic temperature as
, using the semilength along the radial direction as the characteristic length. In the case of study, this strategy would provide an unrealistic temperature difference of several thousands of degrees and
or more. This assumption is only valid (although very conservative) when the cooling comes from the sides of the system (e.g., a cooling plate of the Dual Coolant Lithium Lead blanket (DCLL)). Indeed, the expression does not take into account the cooling effect of the internal tubes which will keep temperatures differences in more moderate values, reducing
as well. More realistic estimations can be made with this formula considering only the PbLi volume contained in between two rows of tubes. With this strategy,
K and
.
In purely natural or free magneto-convection problems, buoyant forces are balanced with Lorentz forces. The characteristic velocity scale is given by the ratio between the
and
:
[
9]. Alternatively, in mixed-convection problems, the velocity scale is given by the characteristic velocity scale of the pressure driven flow (
). In these scenarios, the ratio between
and
weights the relation between buoyant and Lorentz forces [
10]. Indeed, normalizing
,
,
p,
,
and
T by
,
,
,
,
and
, respectively, Equation (
1) can be written as follows:
where the symbol
represents a normalized variable (e.g.,
). Inertial effects are weighted to the Lorentz forces with the usual interaction parameter (
) or by the square of the so-called Lykoudis number (
) in natural convection problems:
. In such problems, only for
of the order of
, inertial effects are significant.
In natural convection situations, the Péclet number () can be related also with and : . Therefore, the advective heat transfer is strongly suppressed by the magnetic field. This effect is also true on mixed-convection problems although a sufficiently high can as usual make advection a significant heat transfer mechanism.
Both natural and mixed-convection scenarios have been studied in the past for different combinations of
,
and
under the influence of a transversal magnetic field. MHD flows in vertical channels with a non-homogeneous heat source and insulated walls have been studied within the framework of the DCLL blanket [
11,
12]. Above a critical value of
, the magnetic field is able to stabilize both the buoyancy-assisted (upward) flow [
13] and buoyancy-opposed (downward) flow [
14]. The latter requires more intense magnetic fields for the stabilization.
In the case of horizontal ducts, both natural convection [
15] and mixed-convection scenarios [
10] have been studied considering heated surfaces. Volumetric heat sources in horizontal channels have been studied for natural-convection regimes within the framework of the Helium Cooled Lithium Lead (HCLL) blanket concept analyses [
16]. In this case, the cooling plates and the spatially varying heat source originate the flow movement. More recently, studies dedicated to the WCLL blanket show that the presence of the cooling tubes immersed in the PbLi also plays an important role on determining the buoyant recirculation patterns [
17,
18].
The strong magnetic field and heat source present in the frontal part of the WCLL blanket implies a very challenging problem from the computational point of view. Being able to resolve the MHD boundary layers while capturing the dynamics of the buoyant vortexes requires very small mesh sizes and very small time steps in a complex geometry. To relax the computational requirements, this work considers reduced values. The final objective would be to gradually increase until reproducing the real blanket conditions. In this paper, two different situations are compared: and . In each case, the external magnetic field has been tuned accordingly.
To keep the ratio between Lorentz and buoyant forces as similar as possible to the WCLL conditions, the gravity field has been reduced as well keeping the ratio between
and
constant. For this purpose, the reduced gravity field (
) and reduced magnetic field (
) follow the following scaling rule:
Equation (
10) will only keep the ratio
constant if the characteristic temperature difference of the reduced case (
) is equal to real
. Since the heat source is the same, this is expected to be approximately true. Results obtained in
Section 3 for different values of
and
are consistent with this approximation.
Simulation Conditions
The material properties of the PbLi [
19] and the steel (EUROFER [
20]) at
K are exposed in
Table 1.
Additionally, the geometrical inputs taken from the WCLL frontal region are exposed in
Figure 4. The toroidal length of the domain is 243 mm.
Regarding the boundary conditions, convective ones (
) have been used in the internal surfaces of the tubes and in the external side of the FW to represent the cooling effect of the water. The water stream temperatures (
) of both circuits are assumed at 311.5 °C. This value is the average temperature of the water thermal cycle [
18]. The heat transfer coefficients (
h) are 11,175 W/m
K for the tubes and 22,012 W/m
K for the FW. The values are obtained using the Dittus-Boelter correlation taking into account the design water velocities of both circuits. Electrically insulated boundary conditions (
) are applied in these surfaces as well. In the inlet channel, the PbLi enters at a constant temperature
K and with a constant velocity
mm/s.
Periodic boundary conditions are applied in the top and bottom steel plates to reproduce the presence of an analogous WCLL cell at the top and at the bottom of the studied one (translational symmetry). Likewise, periodic boundary conditions are applied in the lateral radial-poloidal surfaces since the PbLi parallel circuits stacked along the toroidal direction are supposed to be thermally similar and circulations between parallel circuits are expected to be negligible.
In the internal PbLi/steel surfaces, continuity of the temperature, heat flux, electric potential and electric currents are considered.
Concerning the initial conditions, a pure conductive heat transfer model (treating the PbLi as a solid material) has been used to obtain an initial temperature map. Using this map as the initial condition reduces the time needed to reach relevant conditions, which reduces significantly the computational time. This map is exposed in
Figure 5. In agreement with previous studies [
21], the hottest region is located at the end of the radial channels while the cooling tubes keep the frontal region at moderate temperatures (600–700 K).
As mentioned in
Section 1, the geometry has been simplified in order to employ a structured mesh. A multi-block structured mesh has been designed for this purpose using ICEM-CFD. O-grid structures are included in the vicinity of the tubes. After some testing, it was preferred to extend the radial elements of the tube walls towards the fluid domain in order to ensure good resolution in the boundary layers close to the tubes. For the same reasons, a hyperbolic clustering of cells has been applied towards the walls.
Figure 6 depicts a zoomed view of the computational mesh used for the
Ha = 1000 case.
3. Computational Results
Transient magneto-convective simulations have been performed using the MHD solver integrated in the ANSYS-Fluent platform. This solver has been tested against the experimental results in a NaK loop [
22] under pure MHD conditions. Moreover, a benchmarking exercise with another 4 MHD codes was successfully conducted under magneto-convective conditions [
23].
Results are presented for the
Ha = 1000 case and for the
Ha = 2000 case after 200 s and 100 s, respectively. The gravity field has been scaled according with (
10). In both cases, the 3D solutions obtained present a quite good symmetry along the toroidal direction. This is expected not only because of the toroidal symmetry of the geometry but also because the magnetic field tends to align the convective vortexes along its direction. This kind of behavior is called quasi-two-dimensional (Q2D) turbulence [
24]. For simplicity, results are presented only in the central radial-poloidal plane. Small deviation from these results can be found in other radial-poloidal planes but they are not significant. Weak temperature toroidal dependence has been found as well in hydrodynamics works in the regions where the tubes are toroidally oriented [
21,
25].
Figure 7 depicts a comparison between the temperature distribution of the purely conducting model (or the initial condition) and the results for
Ha = 1000 and
Ha = 2000 cases in the frontal part of the domain. Convective heat-transfer distorts the conductive temperature map. However, the effect is quite moderate which points to a heat transfer scenario dominated by conduction.
In any case, the flow motion boosts heat transfer near the cooling tubes decreasing peak and average temperatures. High temperature regions are slightly displaced towards the upper part of the domain due to the buoyancy force. There are relatively small but appreciable differences between both magneto convective cases. This indicates that even when keeping the overall ratio between buoyancy and Lorentz forces, local effects play a role on heat transfer. For example, higher velocity jets are developed next to the conducting walls at higher numbers.
The computed velocity vector field is exposed in
Figure 8. In the first and second row of tubes, medium sized vortexes, (size comparable with the tube diameter) appear at both sides of each tube. Similar structures also appear next to the FW as a result of being the only cooled wall of the system. The vector field is qualitatively rather similar in both cases. Quantitatively, the velocity scale is of the order of ∼
m/s (
) in both cases as well. However, peak velocities are higher (factor 2) in the
case as a result of the electrical interaction between the fluid and conducting tube walls.
Results also show a weak connection between the pressure driven flow in the radial channels and the rest of the domain. Indeed, a big recirculation region is developed in between the third row of tubes and the radial channels blocking the flow. As a result, the majority of the flow that comes from the inlet channel goes to the outlet channel through a narrow region close to the tip of the baffle plate. This effect was observed in pure MHD flows in simpler geometries that consider no tubes [
26]. Recirculations regions are also observed in both the inlet and outlet radial channels. To confirm their apparation, the recirculation in the radial channels needs to be further investigated. In the present analyses, the proximity of the PbLi inlet and outlet faces might be introducing unrealistic effects. Other studies of the WCLL [
17] predict recirculations in the radial channels as well. However, those seem to be less pronounced. The dissimilarities might be caused by differences in the heat source or the geometrical approximations employed in both works. Recirculations in the radial channels have been also computed for the HCLL TBM [
16]. Nevertheless, the latter work considers a natural magneto-convection flow influenced by the HCLL horizontal cooling plates. These conditions are significantly different than in the WCLL radial channels.
Figure 9 depicts the electric potential distribution in both cases. Electric potential differences arises at both sides of the tubes. Comparing the potential contours with the vector velocity field, it is observed that local maxima and minima of potential are located in the center of the vortical structures. The vortexes extend along the whole toroidal direction of the computational domain. The alignment of the vortexes with the magnetic field direction is a characteristic of Q2D flows. In this flows, the contours of the electric potential coincide approximately with the stream lines of the flow. Indeed, assuming that the flow can be described in 2D, the electric potential is proportional to the stream function (
):
3D potential iso-surfaces are exposed in
Figure 10. The iso-surfaces are a good representation of the Q2D flow whose vortical structures are elongated along the magnetic field direction. This behavior is allowed and enhanced by symmetry of the geometry along this direction. This kind of elongated potential isosurfaces and their relation with the velocity stream-lines have been also obtained for the former EU-HCLL TBM conditions where the cooling of the PbLi was made by some radial-toroidal cooling plates [
16].
The sizes of the potential isosurfaces are dependent on the ratio
. Increasing the magnetic field while keeping
constant is expected to reduce the size of the vortexes [
11].
4. Estimation of the Nusselt Number
Estimations of the Nusselt number (
) are of great interest for design purposes. Knowing
or the heat transfer coefficient
h that appears in its definition allows performing system level heat transfer analyses.
represents the ratio between the convective and conductive processes in the fluid-solid interface:
The heat transfer coefficient can be estimated from the computational results by evaluating the average heat flux at the surface (
[W/m
]) and the superficial temperature (
) of each tube (and the FW) of the system:
The temperature is defined as a temperature sufficiently far away from the interface. Before computing, it is always necessary to specify the way in which the characteristic length (L) and the temperature far away () have been defined. In this case, has been defined as the average temperature of a cylindrical surface (a plane in the case of the FW) at a distance from the tube of interest. The tube surface and the cylindrical surface form a hollow cylinder (a prism in the case of the FW). The characteristic length is defined as the ratio between the volume and the area of external surfaces ().
The distance between the tube and cylindrical surface () has been picked using the purely conducting model. Indeed, it is defined as the minimum length needed to obtain = 1 with the conductive model results. This way, the definition is consistent with the extreme case of purely conductive heat transfer scenario. As a consequence of this definition, the characteristic length is slightly different for each tube but it is around 10 mm in every case which is roughly half the distance between tubes.
Table 2 depicts the heat transfer coefficients and
computed from the simulations results. The tubes of the different rows are labeled from the top to the bottom of the WCLL cell. This way the tube 1 of the row 1 is the tube located in the top left corner of
Figure 4 and tube 4 of row 3 the one located in the bottom right corner.
In agreement with the qualitative conclusions deduced from the obtained temperature maps (
Figure 7), the quantitative analyses confirm that heat transfer is mostly ruled by conduction in the studied region. Indeed,
is very close to unity for all interfaces. There are small differences between the tubes and it is clear that the third row of tubes is slightly more affected by convection than the others. The third row of tubes is closer to the radial channels and therefore it is more influenced by the PbLi flow coming from them. In the other tubes, the flow is essentially moving by natural or free convection and unaffected by the flow of the channels.
It can also be deduced that is around 5–10% higher when increasing (and ). This small difference is related with the increase of the velocity jets in the vicinity of the conducting walls with . In the real blanket conditions ( and the real gravity field) the heat transfer by convection might play a more important role but it is not expected to be comparable with the conducting mechanism.
5. Estimation of the Grashof Number
The value of the Grashof number is dependent on the characteristic length (
b) and temperature differences (
) picked in the definition (
8). In the case of study this definition is not trivial since there are multiple heat sinks (each tube and the cooled FW) and the heat source is extended along the whole domain in a non-homogeneous way.
A conservative approach is considering the maximum temperature difference inside the PbLi (approximately 180 K) and the total radial length of the frontal cavity as characteristic length. With this strategy it is obtained that in the case and in WCLL conditions. If the system were ruled by pure natural convection heat transfer, the Reynolds number would be which is of the order of the average obtained in the simulations.
Alternatively, a local definition can be used based on the temperature differences between some previously defined toroidal-poloidal planes located at different radial positions. This definition is motivated by tubes disposition in rows and by the radial dependence of the heat source.
The planes are defined in the middle of each tube row and in between them. This means that distances between two adjacent planes are approximately the size of the Q2D vortexes obtained in the simulations. Multiple values of
have been derived considering the differences between the planes average temperatures (
) and the distance (
b) between them. The
computed this way are presented in
Figure 11 for the
case.
The values presented are calculated for a reduced gravity field (
10). For obtaining the results in the real WCLL conditions,
has to be multiplied by the square of the ratio between target
(8800) and the reduced
(2000), in other words
. After the re-scaling, the values obtained in the real conditions vary between
and
, depending on the planes considered. Most likely, the most relevant
are the ones computed using two adjacent planes since this is the size of the convective vortexes associated to the problem. This implies a maximum value of
in between the third row of tubes and the radial channels and a minimum value of
in between the FW plane and the first plane. It is worth noting that the
obtained in between two rows of tubes is of the order of
which is one order of magnitude below the initial estimations.
6. Conclusions
This work presents a magneto-convective simulation of the EU-WCLL blanket central outboard elementary cell. The analyses are focused on the frontal region which is close to the FW and subjected to a very high volumetric heat source and to a very high magnetic field. Geometrically, this region is characterized by a group of cooling tubes that crosses the fluid domain mostly in the perpendicular direction to the flow (parallel to the magnetic field). Simulations have been implemented using the MHD solver of ANSYS-Fluent.
According to the results, the dynamics of this region is driven mostly by natural convection. Little influence of the colder radial PbLi flow was found in the frontal region. Indeed, most of the pressure driven flow goes from the inlet to the outlet channel very close to the tip of the baffle plate interacting weakly with the rest of the domain. Medium size vortexes (dimension of the order of the tube diameter) appear in the frontal part of the domain. Having vortexes isolated in the frontal region can have important implications for tritium transport since high concentrations can arise in these regions.
Heat transfer is mostly driven by conduction in the system. Indeed, temperature profiles are relatively similar to the purely conductive solution used as an initial condition. Hot spots are a bit distorted and displaced towards the top part of the domain but not in a strong way. This effect is a consequence of the magnetic field that suppresses the heat transfer via advection. The suppression of the advective heat transfer mechanism by the magnetic field was an expected result in agreement with previous experiences. Despite the conductive nature of the heat transfer, the dense network of water tubes are able to keep the PbLi and more importantly the structural steel below critical temperatures (550 °C).
Nu numbers and heat transfer coefficients have been computed for each cooling tube and for the FW. Values very close to unity have been obtained for in every case. numbers have been also estimated using the simulations outcomes. A global of ∼ has been obtained from the maximum temperature difference obtained in the WCLL cell. Local estimations of in between rows of tubes (∼size of the vortexes) provides a value two or three orders of magnitude smaller: .
Future developments should be focused on reaching the actual values of the magnetic field and gravity field in the blanket. Moreover, the influence of the tube curvature should be evaluated since they will break the toroidal symmetry and possibly the Q2D structures. In fact, it is expected that curvature of the third row of tubes will break the toroidal extension of the bigger vortexes next to the baffle plate. Moreover, some of this tubes penetrate in the PbLi radial channels which will mitigate the hot spot located next to the end of the baffle plate. The effects could significantly affect the recirculation found in the radial channels. This circumstance has to be further investigated.
To include the tube curvature, unstructured tetrahedral elements will be most likely unavoidable. An investigation of the numerical stability of magneto-convective problems with this kind of meshes should be performed. Finally, different orientations of the gravity and the magnetic field should be analyzed. This would allow studying other scenarios in WCLL cells different from the central one. These have variable orientations with respect to the horizontal plane. Besides, the poloidal component of the magnetic field might play also a significant role on magneto-convective results.