3D MHD analysis of prototypical manifold for liquid metal blankets

The water-cooled lead lithium (WCLL) and the dual-cooled lead lithium (DCLL) are two of the breeding blanket (BB) concepts that the EUROfusion consortium is pursuing in the framework of the development of the fusion reactor industrial demonstrator DEMO. Both involve the use of a liquid metal (LM) as working fluid, the lead-lithium eutectic alloy (PbLi), due to its excellent thermal properties and the possibility to serve as both the blanket coolant and tritium breeder and carrier. Unfortunately, due to the high electrical conductivity of LMs, their motion is influenced by the magnetic field used in the reactor to confine the plasma, generating a complex phenomenology which is studied by magnetohydrodynamics (MHD). In this work, a representative prototypical manifold of a BB bottom feeder is investigated for different configurations with the custom phiFoam solver, capable of simulating unsteady, incompressible and isothermal MHD flow. The aim of this study is to investigate which configuration minimizes the flow imbalance in the manifold for the WCLL or in the poloidal breeding zone channels for the DCLL. The distribution of the flow rate among the channels is strongly influenced by the position of the feeding pipe (FP) and by the development of the MHD internal layer near the expansion, which generates important jets close to the lower plate and the upper one, where the channels are attached. The channel aligned with the FP is the one carrying most of the flow, from 55% to 82%, while in the more distant one the flow is almost stagnant, carrying from 17% to 6% of the total flow rate. The total pressure loss is also estimated and its functional dependence on the manifold configuration is discussed.


Introduction
The breeding blanket (BB) is a critical reactor component which has the fundamental tasks of breeding the tritium * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. required for the reactor operation and of extracting the thermal power generated by the fusion reactions. Some of the most promising concepts involve the use of a liquid metal (LM) as working fluid, the lead-lithium eutectic alloy (PbLi), due to its excellent thermal properties and the possibility to be used as both coolant and tritium breeder/carrier. The EUROfusion consortium is pursuing two LM BB concepts for its demonstrator reactor DEMO: a so-called driver BB, the watercooled lead lithium (WCLL) [1], characterized by a lower technological risk and to be deployed in the first phase of operation, and a prospective BB, the dual-cooled lead lithium (DCLL) [2], with attractive features and higher performance for deployment in first-of-a-kind commercial fusion reactors.
One of the main concerns for the BB design is the magnetohydrodynamic (MHD) pressure drop, caused by induced Lorentz forces, that can exceed by several orders of magnitude the hydrodynamic loss [3]. A reliable estimate of pressure drop is essential for the development of LM BB concepts. MHD pressure loss for fully developed flow in straight channels is well understood but accurate correlations are still lacking for concentrated (or 3D) pressure losses in complex geometrical elements, i.e. a bend or sudden cross-section variations. These are associated to the generation of 3D currents and have a functional dependence more challenging to discern, being strongly affected by the specific geometry considered.
Contractions and expansions were probably the first 3D geometries considered as proven by the seminal experimental studies of Branover et al [4] in the late sixties. They analyzed the distribution of pressure along the duct and in proximity of an expansion considering a transverse magnetic field and liquid mercury. Regarding bends, Stieglitz et al [5] experimentally studied the flow in a 90 • bend that moves the flow from the direction perpendicular to the magnetic field to the parallel one. A striking feature of this case is the generation of internal layers (ILs) aligned with the magnetic field and propagating from the corner. This configuration has the highest pressure loss whereas, as observed analytically by Walker [6] and numerically by Bühler [7], if the curve occurs on a plane perpendicular to the field, pressure loss is much lower. This behavior is found also for cross-section variation in the plane perpendicular to the imposed magnetic field, as shown by Molokov [8], with regard to duct expansions parallel to magnetic field.
Analytical studies are conducted in the asymptotic limit of high magnetic field, B → ∞, notable early works being [9,10], and are useful to gain insights on physical phenomena, despite not being representative of reactor conditions. For instance, Bühler [11] investigated the 3D MHD flow in a expansion and the effect of the expansion length ranging from a smooth diffuser to a sudden aperture. Studies that consider more realistic ranges of characteristic parameters are possible only by direct numerical simulations that resolve the MHD layers. Mistrangelo [12] investigated a sudden expansion with finite wall conductivity showing how, even for a high magnetic field and low velocities, the inertial effects are important. Feng et al [13] considered the effect of expansion length and ratio and highlighted that the 3D pressure drop increases with this latter parameter and approach asymptotically a maximum value above 4.
Manifolds often consist of an expansion and distribution in an array of channels. Some early studies by Moon et al [14] and Tillack et al [15] have shown that the redistribution of the flow in the IL can produce a mass flow rate imbalance among the channels. The geometry taken as a reference is usually composed of a single channel that supplies sub-channels stacked in the magnetic field direction through a symmetrical expansion, also aligned with the field. Morley et al [16] investigated an insulated manifold of this type showing how, even if the flow imbalance tends to become uniform as the ratio between the magnetic and inertial forces increases, the subchannels aligned with the inlet always tend to catch more flow, since also the parabolic profile of the jet that forms in the side layers (SLs) is centered with the inlet channel. Rhodes et al [17,18] characterized the main source of 3D MHD pressure loss, sudden cross-section variation from feeding pipe (FP), for an insulated manifold, and he also constructed a correlation for pressure drop for different data range.
In this paper, a prototypical manifold representing a BB bottom feeder is investigated through phiFoam, a custom OpenFOAM solver based on the icoFoam application. As reference, figure 1 shows the in-vessel PbLi loop of the WCLL that, referring to the outboard (OB) segment, is composed of: a FP, which supplies independently every segment of the blanket with 'fresh' PbLi coming from the ex-vessel loop; a bottom feeder, which is the object of this study and is located in the bottom part of the blanket, receives the PbLi exiting the FP and conveys the flow toward the spinal manifold; a spinal manifold, placed between the water coolant manifold and the breeding zone (BZ), distributes and retrieves the PbLi among the poloidally stacked breeding cells; the BZ, where the PbLi flows mostly in radial direction on horizontal planes; the top collector, which is located in the top part of the blanket and conveys the PbLi from the PbLi manifold to the draining pipe (DP); the DP, attached at about 2/3 of the height of the OB segment, which extracts the tritium-rich PbLi from the blanket and closes the cycle returning the LM to the ex-vessel loop. A more detailed description of this system can be found in [19].
We consider the case of a compact component in which the LM undergoes sudden asymmetrical toroidal expansion, like in the bottom feeder showed in figure 1, and then is immediately conveyed through a right angle bend to poloidal subchannels. The walls bounding the manifold are assumed to have zero electrical conductivity. A geometry similar to the one analyzed in this paper is the one studied by Mistrangelo and Bühler [20] and by Chen et al [21], where, however, the expansion is symmetrical and the distribution in the poloidal sub-channels does not take place immediately after the expansion but after an entrance length. The configuration considered in this study, shown in figure 2, is often encountered in BB designs, see [1,2], but is seldomly treated in numerical simulations in favor of ones with more limited geometrical complexity, like in [17,18].
A parametric analysis is performed to assess the effect of inlet position on the flow distribution and the pressure loss, finding that the configuration with a single central feeding channel, compared to the ones with two lateral channels, is the configuration that maximizes both the non-uniformity of the flow rate among the channels and the pressure drop. Regarding the first point, it is emphasized how obtaining a uniform flow rate among the channels is important both for the transport of the tritium, reducing the possibility of its accumulation in the stagnation zones and, above all, for the BBs in which the LM is used as a coolant, as the DCLL. Figure 2(a) shows the geometry considered by our numerical model, which is representative of the basic component layout envisioned by BB designs, even if actual complexity may vary greatly across them. LM path can be schematized as an abrupt expansion from a small FP (characterized by a toroidal size 2 a, a poloidal one h and an axial length L in = 4a) into a main chamber (e, h, and c) with a simultaneously distribution in poloidal sub-channels (3 in figure 2(a), with sizes 2 d, L ch = 6a and c). Among the parameters listed in table 1, an important one is b i , i.e. the toroidal distance between inlet pipe center and furthest chamber wall which governs the main flow features through the position ratio r p = b i /e. As shown in figure 2, we consider three component layouts (A, B, and C) characterized with different values of the b i parameter.

Governing equations, dimensionless groups and physical assumptions
The computational fluid dynamics code used for this study is the phiFoam solver [22], a custom OpenFOAM tool for the resolution of unsteady, incompressible and isothermal MHD flows that employs the ϕ-formulation and the current density (j) conservative scheme developed by Ni et al [23] in order to ensure charge conservation (∇ · j = 0). The implemented governing equations are reported in dimensionless form below [24]: whereũ = u/u 0 is the scaled velocity,p = p/(ρu 2 0 ) the scaled pressure,B = B/B 0 the scaled magnetic field,j = j/(σu 0 B 0 ) the scaled electric current density andφ = ϕ /(L c u 0 B 0 ) the scaled electric potential field. The quantities u 0 , ρ, B 0 , σ and L c = a are, respectively, the inlet average velocity, the density, the imposed magnetic field, the electrical conductivity and the characteristic length.
The interaction parameters N = σB 0 2 L c /(ρu 0 ) gives the ratio of electromagnetic to inertial forces, whereas the Hartmann number Ha = √ NRe = L c B 0 √ σ/(ρν), with Re = u 0 L c /ν the Reynolds number and ν the kinematic viscosity, represents the ratio of electromagnetic to viscous forces. It also characterizes the thickness of boundary layers appearing close to walls normal to B 0 (figure 2(a), the so-called Hartmann   walls) and to walls parallel to B 0 (side walls). In the first case, the Hartmann layer (HL) width is δ HL = L c /Ha. For the SL instead δ SL = L c / √ Ha [24]. As previously mentioned, abrupt geometric variations are associated with the generation of 3D currents which, instead of closing in the walls of the cross-section as for a fully developed flow, may also close axially inside of the fluid. Furthermore, as shown by studies of Hunt and Leibovich [25], an MHD IL, or Ludford layer, can be formed in the vicinity of a sharp change in wall curvature and can causes significant flow redistribution and pressure drop. Indeed, even in the presence of a flow that is globally inviscid (Ha ≫ 1) and inertialess (N ≫ 1), the IL can still be affected by inertial and viscous forces and its thickness δ IL is controlled by both Ha and N. This leads to three possible regimes and distinguished by the value of the regime ratio R r = N/Ha 3/2 : an inertial-electromagnetic (IE) regime for R r ≪ 1, a viscous-electromagnetic (VE) regime for R r ≫ 1 and an IE-VE regime for R r ≈ 1 [24,25].
An adiabatic, isothermal and pressure-driven PbLi flow, which physical properties are collected in table 2, has been considered to study the magneto-hydraulic features of the bottom feeder. The imposed magnetic field B 0 is toroidal (figure 2(a)), stationary and uniform, and the walls are perfectly insulated. As can be seen from table 3, all the analyzed cases are characterized by governing parameters placing them below the threshold that mark the transition from laminar to quasi-two dimensional (Q2D) turbulent regime (Re = 65Ha 1/2 , figure 9 in [3]), therefore the flow is assumed to be laminar.
The regime ratio considered in our model is equal to 2.055 × 10 −2 and falls in the IE regime (R r ≪ 1), which is a typical value for separately cooled blankets and it is of particular importance since it cannot be studied relying on asymptotic methods [22]. This parameter is kept constant and, as the Hartmann number varies, we obtain an average velocity u 0 to be imposed on the inlet, as collected in table 3. Other relevant parameters are presented in the same table for the three Hartmann number considered in this study. To provide a frame of reference, those of the two LM BB concepts developed in the European fusion programme, the WCLL and DCLL, are also shown [22].

Domain discretization and numerical model
The governing equations (1)-(4) have been discretized on a hexahedral, structured and non-uniform mesh with greater refinement next to the walls and in proximity of the expansion and the channels contractions to resolve the MHD layers (HL, SL, IL). The number of divisions in these ones, following the result of the mesh sensitivity study reported in section 3 for a similar three-dimensional case, is set equal, respectively, to 3, 8 and 7. Given the symmetry of the geometrical component, as shown in figure 1 for the WCLL bottom feeder, and of the expected flow solution, only half manifold has been simulated.
Regarding the boundary conditions (BCs) shown in figure 2(a), the fully developed velocity profile (u FD ), calculated by a dedicated simulation on a channel with the same geometry and governing parameters, is imposed at the inlet whereas, at the three outlet channels, an outflow BC with zero reference pressure is imposed. The noSlip BC has been applied to all the walls, except for the symmetry patch (see figure 2(a)), and, for all boundaries, the zero-flux BC is imposed for the electric potential [22].
The temporal discretization was performed with the Crank-Nicolson scheme while the gradient and Laplacian of all the variables were spatially discretized, respectively, with the cell-limited second order and with the Gaussian integration for central difference scheme. For the advection scheme, the Gaussian integration with the linear-upwind divergence scheme was used and a linear interpolation was used for all the variables. About algebraic solvers, the preconditioned conjugate gradient solver with the simplified diagonal-based incomplete Cholesky preconditioner has been employed for both the pressure and electric potential equations, while the preconditioned bi-conjugate gradient (PBiCGStab) solver with the simplified diagonal-based incomplete LU preconditioner has been used for the velocity equations. Three PISO cycles and three external electric potential cycles were performed for each time-step [22], which value was set to guarantee the Courant-Friedrichs-Lewy condition [27], considering a maximum Courant number equal to 0.9. The simulation was performed until reaching steady state solution, determined by the negligible change of some control quantities, such as the pressure drop, the output flow of the channels and the velocities in some significant points of the domain. Regarding the internal convergence criteria, a tolerance from 1 × 10 −6 to 1 × 10 −9 was chosen for all the variables.

Verification and validation
Before continuing with the presentation and discussion of the manifold analyses, the main outcomes of the 3D validation carried out on the phiFoam solver is reported. Further details, as well as the 2D benchmark for Hartmann number up to 5000, can be found in [22]. The 3D benchmark aims to test the performance of the solver with complex geometries such as those present in the blanket. Local velocity and pressure profiles are considered for the purpose of mesh sensitivity but the focus is on assessing the reliability of the code for the prediction of integral quantities (flow distribution and 3D pressure drop), which are the ones directly impacting the blanket design. Figure 3 shows the geometrical model representing the experimental section considered in [28], which consists of a liquid mercury flow through an expansion and then its distribution in a set of three parallel rectangular channels stacked in the direction of a uniform magnetic field B 0 . The entire section is made of acrylic and was therefore considered perfectly insulated. Since the expansion region is the one in which 3D MHD effects predominantly exist, to better characterize these effects, the half-width b (see figure 3) and the average velocity u 0 at the expansion have been chosen, respectively, as the characteristic length and velocity in the calculation of the dimensionless groups. The geometrical dimensions and physical properties used can be found in [22]. Since the 3D currents couple the flow at the expansion with that at the inlet, it is important to consider the expansion ratio r ex = b/d, which in this case is equal to 4. The computational grid was generated with the same methodology described in section 2.3, shown in figure 4, as well the numerical setup and BCs, that are the same in section 2.3 without the symmetry BC since, for this case, was considered the full geometry shown in figure 3.

Mesh sensitivity
To evaluate the independence of the results from the resolution of the computational grid, a mesh sensitivity study was carried out considering five meshes with different resolution for the Hartmann, side and ILs, as shown in table 4, and considering as control parameters the distribution of the flow rate among the three channels, the pressure drop and the value of the local velocity detected near the expansion (point p in figure 3) which characterizes the IL. The pressure drop is calculated as the surface integral average in the inlet section, since a zero pressure is imposed at the outlet, while the flow rate is obtained by the integration of axial velocity in channels cross-sections. As can be seen in table 4, already with the distribution of nodes in the layers used in M2, it is possible to correctly characterize the flow, presenting a maximum difference relative to   [17], an IL is generated upon expansion which conveys the flow near the side walls, generating an M-shape velocity profile both in a direction perpendicular to the magnetic field and parallel to it, the latter shown in figure 5(a). These 3D MHD phenomena occur also at the contractions due by the three channels, as can be seen from the jets generated near the internal wall of the external channel (peak at about 0.044 in figure 5(b)) and near both the walls of the internal channel (peak about 0.034). In the middle of the channels, the flow becomes uniform as the classic slug profile for the insulated channels, with a slightly higher core velocity for the central channel (figure 5(c)). Subsequently, the expansion due to the outlet box produces ILs at the edges of the channels (figure 5(d)), similar to preceding expansion but less intense, and then the flow tends again to its fully developed state away from this last perturbation (figure 5(e)). As regards the pressure distribution along the manifold, shown in figure 5(f ), the local maxima and minima occur near the geometric discontinuities, where additional Lorentz forces are generated due, in turn, by the generation of the 3D currents. Due to the different direction of circulation of the currents, the Lorentz force generates minima in the pressure at expansions (at about 0.22 and 0.5 in figure 5(f )) and a local maximum at contraction (at about 0.35).

Benchmark results
First, the flow rate distribution between the three channels is compared to the one measured by Messadek and Abdou [28]. In [28], several experimental tests are reported, with an Hartmann number and Reynolds number ranging between 421 < Ha < 2429 and 4000 < Re < 120 000, almost all in the Q2D MHD regime. Since the use of a DNS mesh has been deemed to be too computationally expensive and a dedicated model for the Q2D regime has not yet been implemented in phiFoam, the data set at Ha = 1503 and Re = 6400 has been selected for validation purposes. It should be noted that this is the only test case available in laminar conditions, even if relatively close to the laminar-Q2D transition. In order to obtain an estimate of the flow rate imbalance, the data plotted in figure  3 in [28] has been digitized and a linear interpolation was performed [22].
The benchmark results are reported in table 5. As can be seen, there is a consistent error on the left channel flow rate, about 12%, and a smaller value relating to the other two channels (≃5%). The solver returns an equivalent flow distribution between the three channels which is, in theory, the expected result according to what is reported in [28]. This is also confirmed by the numerical analyses reported in [17], where in the range of parameters 1000 < Ha < 2000, 50 < Re < 3750, r ex = 4 and 0.5 < L ex /b < 3 (L ex /b = 1 in our case), the difference in percentage of the flow in the channels was calculated to be at most on order 1 × 10 −3 %. This leads us to believe that the discrepancy found with the experimental data must be attributed more to the uncertainty of the digitalization procedure and interpolation to extract the data from the original figure, rather than being associated to a fault in the physical prediction by the code. A  discrepancy of ≈5% is therefore taken as a representative behavior of our model in the prediction of the flow rate repartition.
To quantify how much the transition to a Q2D regime may affect the prediction of the flow distribution in our numerical model, simulations were performed in a parameter space consistent with the initially disregarded experimental data from [28]. We observed deviations with respect to the flow rate distribution among the channels from 20% to 30% or even the outright divergence of the solution.
In the second part of the benchmark, the control parameter chosen is the 3D MHD pressure drop due only to the expansion, reported in [17], where the authors performed several numerical simulations in the range of parameters of 4 < r ex < 12, 1000 < Ha < 6570 and 50 < Re < 2500. By subtracting from the total numerical pressure drop the 2D pressure drop calculated using the Shercliff formula [29], they obtained the increase due to the 3D losses generated at the expansion and then these data were sorted in groups with equal expansion ratio obtaining an analytical formula for the prediction of 3D losses for this geometry and regime, finding that the 3D pressure drop scales linearly with ρu 0 N 2/3 for the IE regime, in excellent agreement with the predictions previously made by Molokov [8].
Following the same methodology, the 3D pressure drop for only the expansion, shown in figure 6, has been evaluated for Ha = 1503 and for Re = 942, founding an underestimation of ≃8.6% with respect to the Rhodes formula [17]. This discrepancy can be attributed to the expansion ratio of the case used (r ex = 4) which, even though it falls within the range considered in [17], is in any case the lower limit and therefore the one that could be most affected by the interpolation error.  Figure 7(a) shows the distribution of the scaled electric potential (ϕ 0 = L c u 0 B 0 ) in the yz-plane which cuts the expansion in half (x = b A , see figure 2(a)), whereas figure 7(b) shows the topology of the currents by means of streamlines of the scaled electric density current (j 0 = σu 0 B 0 ) generated near the inlet, the expansion and the outlet channels at Ha = 2000 for configuration A. As expected, the electric potential varies axially in correspondence of the expansion, where we expect a local velocity gradient caused by the enlargement of the flow area. Three-dimensional electric currents close inside the fluid propagating forward and backward from the expansion in the toroidal direction (x-direction), but also partly into the vertical channels, electrically connecting the three fluid zones. Indeed, the currents in the inlet channel, despite having imposed a Shercliff velocity distribution, are not perfectly 2D, given the proximity of this section to the area of influence of the 3D currents. The same happens for the outlet section where the greater distance from the interface with the feeder makes this phenomenon less conspicuous. Figure 8 shows the scaled (F L0 = σu 0 B 2 0 ) z-and ycomponent Lorentz force distribution superimposed to the velocity streamlines in the xz-plane at the middle of the inlet channel (y = h/2, figure 8(a)) and in the yz-plane at x = b A ( figure 8(b)) for Ha = 2000. The y-direction currents interact with the magnetic field (x-direction) generating a stream-wise Lorentz force ( figure 8(a)) that brakes the fluid in the center of the channel, while the axial currents (z-direction) generate the F L,y component which pushes the flow toward the side walls, as shown in figure 8(b). This generates an IL in proximity of the expansion which propagates parallel to the magnetic field lines (figure 8(a)) in agreement with the Ludford's theory [25]. This distribution of the Lorentz force causes the PbLi to avoid the center of the channel, being instead conveyed near the side walls, and strongly influences the velocity distribution downstream of the expansion, as shown in figure 8(b).

Configuration A
The scaled stream-wise velocity distribution for the Ha = 2000 case is shown in six cross-sections shown in figure 9. At the beginning of the inlet channel (cross-section A in figure 9) the velocity distribution matched the imposed Shercliff profile ( figure 10) showing how, even if the induced currents at that location are slightly affected by the 3D currents generated upon expansion as shown in figure 7(b), their distortion from a purely 2D distribution has no relevant impact on the velocity distribution.
By approaching the expansion the flow moves toward the SLs, as already shown with the velocity streamlines in figure 8, and when it reaches the expansion (cross-section B in figure 9) it separates in two ILs that form on the expansion wall, as shown in figure 11 where is reported the scaled velocity x-component (aligned with the magnetic field) on the two sampling lines shown in figure 9. As can be seen from the latter, the asymmetry of the expansion causes the IL that is generated near the wall of the short side of the case (red line) to be more intense than the one that is generated on the long side (blueline ), in agreement with the distribution of the Lorentz force in figure 8(a).
The concentration of the flow in the ILs generates two highvelocity regions, as shown in figure 12 that refers to crosssection B shown in figure 9. The first one is found close to the lower plate of the case (red arrow in the figure), whereas the second one is observed at the channel connections (blue line). Due to the presence of the immediate connection of the outlet channels to the expansion, the velocity distributions differs greatly between the lower half of the case and the upper one, which leads to a more complex one in the latter region. Near the lower plate (red line in the plot of figure 12(a)), which is less affected by the presence of outlet channels, there is the generation of two jets that have a shape similar to that identified by Rhodes et al [17], Mistrangelo and Bühler [20] and Chen et al [21] for an insulated symmetrical case without channels or with some distance between channels and expansion, but with a imperfect M-shape due to the asymmetry of the expansion. On the other hand, the upper part is strongly influenced by the presence of the channels (blue line in the plot of figure 12(a)) and presents a an inflection point at the boundary between channels C3 and C1 (figure 9), separating the flow in a higher velocity part near the channel C3 and a lower on, near channel C2.  Halfway through the expansion (cross-section C in figure 9), the velocity assumes an almost parabolic shape near the lower plate only on one side of the expansion (red line in the plot of figure 13), instead of across the entire section, as it happens for a symmetrical expansion and without an immediate channel distribution [13,17,20,21]. In the upper part, which corresponds to the position in the middle of the inlet section of the outlet channels, the radial component of the velocity is close to zero while the poloidal component, shown in figure 14(a) considering cross-section D in figure 9, assumes a very complex distribution, which is dramatically influenced by the presence of the wall that separates the channels C3 and C2 ( figure 9).
As shown in 14(a) (cross-sections D, E and F in figure 9), this splits the jet into two asymmetrical sections, as the separating wall is not located in the center of the inlet channel, where the velocity reaches high values in the portion of space inside the cross-section of the inlet channel (from about 0.55 to 0.9 of the red plot in 14(a)) while there is a counter flow in the external portion. In addition to the velocity field asymmetry in the x-direction, the bend of the channel directly attached to the expansion generates a z-direction asymmetry, where the velocity near the internal wall (red line in plot of 14(a)) is considerably higher than the one at the external wall (blue line).
The velocity distribution tends to become uniform as the flow continues along the channels (cross-section E in figure 9), at least in the x-direction as shown in figure 14(b), while the non-uniformity in the z-direction is still considerable. At the channel outlet (cross-section F in figure 9) the velocity distribution returns similar to the classic fully developed slug-type distribution ( figure 14(c)), albeit with a considerable mass flow rate imbalance among the channels.
As expected, having imposed a fixed regime ratio, when the intensity of the magnetic field is reduced the fundamental characteristics of the flow remain unchanged. The similarity of the velocity distribution in the Hartmann number range considered is clearly visible in figures 15(a)-(c), where is reported the scaled velocity y-component (u y /u 0 ) profile in the middle of the channel C3 (see figure 9) along a line that divides it in two in the z-direction (from z = L in to z = L in + c, see figure 2(a)) and at three different y-positions: at the inlet (y = h), in the middle (y = h + L ch /2) and at the outlet (y = h + L ch ). Focusing the attention on the lowest part of the channel ( figure 15(a)), the velocity jet located closer to the expansion is about 10% higher than the one that forms on the opposite wall for all the Ha numbers. A deformed paraboliclike core region, clearly affected by inertial forces, is formed between about 40% and 80% of the channel depth and is characterized by an almost constant velocity across the three cases. On this position, very close to the expansion, the jets have an intensity that increases with Ha but, progressing along the channel (figures 15(b) and (c)), this difference is reduced. At the same time, the central core tends toward the classical slug flow solution for an electrically insulated channel. Figure 16 shows the variation of the scaled pressure (p 0 = 1/2ρu 2 0 ) along the coordinate s, reported in the detail of the figure, which does not consider the length parallel to the magnetic field where the pressure variation is negligible. The pressure has a strong decrease as we approach the expansion and, immediately after, there is recovery due to the axial Lorentz force locally accelerating the flow ( figure 8(a)), in agreement with the results in [17,20]. The pressure continues to increase until the inlet of the channels (786 mm). There some nonlinear behavior is observed due to the 3D effects caused by the contraction of the inlet channels, before the pressure profile decreasing linearly in rest of the channel, signaling the recovery of the fully developed flow. It should be noted that for all simulations, regardless of the configuration and the magnetic field intensity considered, low-frequency oscillations in time of the variables have been observed with respect to an average central value, as shown in figures 17(a)-(c), due to the generation of vortex structures aligned with the magnetic field in the vicinity of the expansion, as shown in figure 18(a) by means of the isosurface at 0.1 of the Q-criterion [30]. Figure 17(a) shows the time evolution of the scaled velocity y-component of the jet peak near the expansion of the channel C3 (red circle in figure 15(a)). This location is the one in which the oscillations have the largest amplitude and has been taken as a reference to evaluate the convergence of the solution, which was considered reached when the oscillations deviate a maximum ±5% with respect to the average value of the solution from t > 100 s. The integral variables, such as the pressure drop ( figure 17(b)) and flow rate in the channels (figure 17(c)), exhibit oscillation with smaller amplitude (±1%), as shown in figures 17(b) and (c) and reach a statistically steady behavior much earlier (t ⩽ 40 s) than the local quantities (t ⩾ 200 s). A similar approach to evaluate the convergence of the solution was also used for configurations B and C, where the considered peak velocity is located in the carrier channel, i.e. the one more closely aligned with the inlet, the channel C2 for the configuration B and the C1 for the configuration C [22]. Figure 17(a) shows that the amplitude and frequency of the inlet oscillations increase with Ha. The explanation for this phenomenon is related to the structure of the wall jets developing at the inlet of the poloidal channel side walls due to the interplay between ILs and the flow incoming from the manifold tank. As it can seen in figure 15 these structures are reminiscent of the wall jets typical of the M-shaped velocity profile occurring for a MHD flow in an electrically conductive duct. In particular, it can be observed that an inflection point is present at the connection with the developing core flow and that the peak velocity is proportional to a power law of Ha. These features are known to be a source for instability due to the development of progressively larger shear stress between the slow-moving core and the thinner and quicker wall jet [31]. In the present case, velocity oscillations in the jet are only weakly suppressed due to the electrically insulated walls and favorable orientation with the magnetic field direction, hence the observed behavior [24]. Fluctuations tend to be dampened as the flow travels through the poloidal channels and regains its fully developed state.
Concerning the expansion vortexes, the largest one is located at the corner formed by the lower plate and the external wall of the case and in proximity to the attachment of the channels, in correspondence with the wall of the expansion. Similar structures are also highlighted in some cases analyzed by Rhodes et al [17]. These vortexes promote a mixing in the expansion, redistributing the flow between the channels as a function of the entrance position of the fluid particles in the inlet channel, as shown in figure 18(b), which represent a map of the velocity streamlines originating from different positions in the inlet channel approximately midway between the inlet and the expansion. The fluid particles along the horizontal midline of the inlet channel (black line in figure 18(b)) are conveyed almost exclusively to the channel C3 (see figure 9), as well as the fluid particles near the upper side wall (blue line in figure 18(b)) which, after having fed the vortex located at the expansion, enter channel C3 adjacent to the wall closer to the inlet channel. The particles that come from the lower side wall (orange in figure 18(b)) penetrate further into the case, gyrate in the vortex located at the corner, and tend to flow predominantly up the channel C1. On the other hand, considering the map in y-direction, the particles entering the symmetry axis of the inlet channel (green line in figure 18(b)) feed mainly the channel C3, but those in the lower area that are able to penetrate more into the expansion end up in the other two channels. The fluid that enters near the Hartmann walls feeds the vortexes that form close to the corners, between these walls and the expansion (see figure 18(a)). The particles that enter near the external Hartmann layer (violet in figure 18(b)) mainly feed the channel C3, while the ones that enter near the internal HL (cyan in figure 18(b)), feed mainly the C2. Channel C1 is the one in which the fewest streamlines arrive and, as will be seen below, is the one in which there is the smallest fraction of the total flow rate.

Configuration B and C
We now briefly consider the flow features of the other two configurations, focusing above all on the differences due to the different position of the inlet channel and referring to [32] for a more complete discussion: configuration B, in which the inlet channel is positioned in the center of our model ( figure 2(b)) and configuration C, in which the inlet channel is positioned on the right-side of our model ( figure 2(c)), forming a symmetrical expansion when considering the whole toroidal extension of the manifold, like in the DCLL bottom feeder.
It must be noted the latter configuration is conceptually different from the others, as the inlet channel of the complete geometry is only one and therefore the available cross-section for LM ingress in the manifold is halved. In analyzing this configuration, it has been decided to keep constant the average velocity imposed at the inlet in the other configurations, as described in section 2.2 and shown in table 3, instead of doubling it to have the same total flow rate circulating in the manifold. This choice has two justifications: the first is that doubling the flow rate would have pushed the MHD flow above the laminar-Q2D threshold, currently not treatable with the reference code as pointed out in section 3, and the second is that this maintains the same balance between the IEviscous force in the ILs as described by R r (see section 2.3). Furthermore, it favors a certain similarity in the MHD flow features that is beneficial to highlight the differences arising due to the inlet position. It is reasonable to expect that the flow rate distribution for configuration C should not change considerably even with double the flow rate currently adopted, since the condition R r ≪ 1 is maintained in our simulation and, as well, the relevance of inertial forces. As such, we expect our estimate to provide some guidance for the flow rate distribution in a parameter space closer to DCLL conditions.
The distribution of the electric potential and of the electric currents are quite similar with respect to the ones of configuration A, with the generation of a complex 3D topology of the currents that couples the inlet channel, the expansion and the outlet channels. This is reflected on the Lorentz force distribution that, also in these cases, tends to pushes the flow toward the side walls, concentrating it in the ILs [22].
With reference to the planes shown in figure 9 for configuration A, completely similar to those considered for the other two configurations, it can be seen from figure 19(a) that for configuration B the velocity distribution generated by the expansion is almost symmetrical, due to the fact that the inlet is positioned in the middle of the model. Similarly to configuration A, the velocity at the lower plate (red line in figure 19(a)) assumes the classic M-shape profile due to the acceleration of the flow in the ILs, while the velocity near the channels (blue line) is influenced by the latter but to a lesser extent if compared with case A (cfr figures 19(a) and 12(a)), due to the aforementioned central position of the inlet channel.
In the middle of the expansion ( figure 19(b)) the jet that forms near the lower plate (red line) in the half of the expansion is almost symmetrical with an almost parabolic shape similar to that shown in [17,20]. The upper jet, on the other hand, is strongly influenced by the presence of the walls that separate the outlet channels, as shown in figure 20(a) (red line). Indeed, the inlet channel is wider than the outlet channel (table 1), and therefore the strongly accelerated flow by the Lorentz force impacts with the edges formed by the attachment of the channel C2. This reflects the jet largely in the channel C2, unlike case A where the jet was split centrally and divided between channels C3 and C2, generating a high velocity zone on the inner wall of channel C2 and reverse flow zones on the inner walls of the other two channels (red line in figure 20(a)). This strongly disturbs the velocity distribution in channel C2, which is peaked close to the wall bordering the expansion also near the outlet while in the other two channels a uniform configuration is achieved, as shown in figure 20(b).
Considering configuration C, the fact that the flow coming from the inlet channel does not meet a separating wall, but can enter channel C1 almost undisturbed, means that the distribution of the expansion velocity near the channels is quite regular, as shown in figure 21(a) (blue line), while near the lower plate we find the M-shape profile (red line). The symmetry of the expansion leads to the formation of a localized almost parabolic jet in the lower part at the middle of the expansion ( figure 21(b)). The upper jet, as already mentioned, is not broken and deflected by any wall in this configuration and practically enters the channel C1 aligned with the inlet entirely, as shown in figure 22(a).
In this configuration are clearly visible the jets that form on the wall opposite to the expansion (blue line in figure 22(a)). Indeed, the Lorentz force concentrates the flow in the lower plate and, when the flow impacts with the aforementioned wall, it rise upwards and collides with the walls separating the channels, forming the jets. The flow remains highly nonuniform in channel C1 up to the exit of channel, as shown in figure 22(b), while in the other channels the velocity becomes almost uniform at a very low average velocity, creating a strong flow rate imbalance among the channels.

Flow rate distribution
As already mentioned in the previous sections, the morphology of the upper jet greatly influences the distribution of the mass flow rate in the outlet channels, preferentially conveying it to the outlet channel aligned with the inlet one. Figures 23(a)-(c) shows the flow rate distribution among the channels for, respectively, configuration A, B and C. For configuration A, the one more similar to the WCLL bottom manifold (see figure 1), the distribution remains roughly fixed as the Hartmann number varies, as shown in figure 23(a), and which shows a strong displacement of the flow towards the channel C3 (see figure 9), that carries ≃54% of the flow rate, to the detriment of channel C1, the furthest away from the inlet, that carries less than ≃10% for high Hartmann numbers. For comparison, the hydrodynamic laminar flow repartition (Ha = 0) was also simulated using the same imposed flow rate as Ha = 1000 (see table 3) for the three configurations, showing a much more balanced repartition. Therefore, the 3D MHD effects due to the expansion greatly influences the imbalance of the flow rate among the channels and must be carefully considered in the evaluations of the tritium transport, that could concentrate in the almost stagnant channel.
Also for the configuration B, shown in figure 23(b), the flow rate is strongly unbalanced in favor of the outlet channel aligned with the inlet (channel C2), which passes from ≃61% to ≃65% from Ha = 1000 to Ha = 2000. Regarding the other two channels, the situation is slightly better as there is no heavily penalized channel. Indeed, for this configuration, the other two channels are equally spaced from the inlet (see figure 2(b) and carry about ≃18% of the total flow rate. However, also for this configuration the MHD impact on the flow distribution is evident.
The worst imbalance occurs for configuration C, as shown in figure 23(c), the reference layout for the DCLL, where the channel C1 carries most of the flow rate, from ≃78% to ≃82% passing from Ha = 1000 to Ha = 2000, and in the channel furthest from the expansion (channel C3) the flow is practically stagnant, since only ≃5% of the flow rate is able to flow into it. Even if the high flow rate places the DCLL in the Q2D turbulence regime that seems to smooth the flow [28], this flow rate is probably too low for acceptable heat transfer.

Pressure drop
As regards the pressure loss, table 6 collects the total (∆p), two-dimensional (∆p 2D ) and three-dimensional (∆p 3D ) pressure loss for all the cases analyzed in this study. For the following considerations, it is useful to consider the expansion ratio for the different configurations, defined as r ex = b i /a for configuration A and B and as r ex = e/a for configuration C. Indeed, although in this model there is no variation in the width of the downstream and upstream ducts before expansion as in the classic cases for which this parameter is defined [13,17,33], the variation in the position of the inlet actually modifies the expansion, increasing the aspect ratio through as the input moves towards the center of the expansion (configuration C). As shown in table 6 and as will be analyzed below, an increase in ∆p 3D correlates with an increase in r ex , in agreement with the literature [13,17,33].
For the configuration A, the total pressure drop increases by ≃4 times going from Ha = 1000 to Ha = 2000, combining both the increase in the intensity of the magnetic field and that of the average velocity, the latter an intrinsic consequence of the assumption of constant ratio between N and Ha 3/2 (table 3). The 2D loss, similarly to what was done in [17] and explained in section 3 for the 3D benchmark, was evaluated using the Shercliff formula [29], considering each section separately (inlet channel, expansion, outlet channel) and adding the three contributions. The 3D loss results from the difference with the total one, computed as the area-average pressure on the inlet since at the outlet has been imposed the zero-reference value (see figure 2(a)). To compare the 3D losses generated at the expansion with the 2D ones, the equivalent length L eq 3D was defined as the length of the outlet channel which generates a ∆p 2D equal to the 3D one given the same flow conditions, scaled with the characteristic length of the outlet channel, i.e. its half-width aligned with the magnetic field (d in table 1). As can be seen from table 6, the weight of the 3D loss on the total one substantially increases with Ha, corresponding to L eq 3D ≃ 70 for Ha = 1000 to L eq 3D 105 for Ha = 2000. Figure 24 shows the scale law of the 3D scaled pressure drop (∆p 0 = 1/2ρu 2 0 ) as a function of N. According with the predictions of the IL theory for the IE regime [24,25], the ∆p 3D scales linearly with ρu 2 0 N 2/3 . Using the same nomenclature employed by Rhodes et al [17], the coefficients k ie and d ie take on the values of, respectively, 1.2400 and −9.5478. These are different with respect to the Rhodes model due to the geometric difference between the two manifold layouts. To extrapolate the 3D pressure drop for the configuration A under the WCLL operative magnetic field (OMF) condition, we assume that the relationship shown in figure 24 remains linear as N increases. This is consistent with what reported in [17] and it is deemed reasonable since R r is held constant and that the flow features are well delineated for Ha = 1000. Given these assumptions, we estimate a ∆p OMF ≃ 414 Pa.
From the 3D pressure loss values, it is possible to evaluate the generic 3D loss coefficient k, collected in table 6, by the relations [3] ∆p 3D = ζ 1 2 ρu 2 0 (5) where ζ is a local MHD resistance coefficient that is function of the N and Ha, other than the problem geometry and magnetic field orientation, and L 3D is the characteristic length scale for the 3D flow [29].
Applying the relation and coefficients k ie and d ie developed by Rhodes et al [17], that is relative to the 3D pressure drop for the expansion only, we obtain a coefficient k lower than ≃20%. This difference can be explained with the difference in the considered manifold layout. The presence of a right angle turn perpendicular to the field, the contraction at the poloidal channels attachment point and separating walls between these channels, which are obstacles that prevent a free expansion of the flow, all contribute to a significant more complex flow than the one observed in [17]. It is likely that the particular arrangement of the poloidal channels has the most significant effect since a bend in a plane perpendicular to the magnetic field should not generate large 3D losses, as it is known from the observations of Walker [6], Bühler [7] and Molokov [8]. Similarly, Rhodes et al [17,18] stated that, when the sum of the cross-sections of the outflow channels does not differ much from the inflow cross-section of the expansion, as happens in this case, the 3D loss due to the following contraction can be neglected and the additional loss being mostly controlled by the parameters of the primary expansion. However, it should be noted that, in this case, the bend occurs simultaneously with the expansion and therefore the velocity, as shown in figure 12, turns on an inclined plane with respect to the magnetic field. It is difficult to reduce our model to the ideal conditions previously described and it seems reasonable that the higher 3D losses observed in our model must be a consequence of the complex interplay of all these effects that perturbate the current distribution.
Similarly, also for the configuration B (table 6) passing from Ha = 1000 to Ha = 2000 the increase in total pressure drop is ≃4 times, while the loss is lower for each Ha at configuration A, in accordance with the lower expansion ratio (table 1) for this case [13,17,33]. Having obtained the 2D and 3D losses with the same methodology presented for the configuration A, also in this case the 3D losses are responsible for almost all the total ones, where also in this case the equivalent length moves in a range of ≃63 ÷ 96, and the scale law of ∆p 3D as a function of N scales linearly with ρu 2 0 N 2/3 [22]. With the same considerations made for configuration A, the ∆p 3D under OMF condition is ≃374 Pa, associated with a k coefficient equal to ≃0.036, as shown in table 6. Applying the relation developed by Rhodes et al [17], we obtain a coefficient k, relative to the expansion only, ≃25% lower than that of the manifold, of the same order with the difference found for configuration A. It must be emphasized that the expansion ratio of this configuration, equal to 2.82, is out of the range in which the semi-empirical relationship proposed by Rhodes was evaluated (4 ÷ 12) and therefore the uncertainty about its applicability is added to the difference caused by the presence of the curve and the channels.
Finally, considering the configuration C, also in this case passing from Ha = 1000 to Ha = 2000 the increase in total pressure drop is ≃4 times. As expected from the higher value of the expansion ratio (see table 1), this is the configuration that presents the largest pressure loss for all Hartmann numbers. The linear tendency of ∆p 3D with respect to ρu 2 0 N 2/3 is confirmed [22]. Considering that this configuration is strictly related to the DCLL bottom manifold, using the same methodology described for the other configurations and the parameters in table 3, the total pressure drop at OMF condition considering is ≃7579 Pa, a significantly higher value than the estimates for the WCLL as the flow rate in the DCLL is much higher. It should be emphasized that DCLL bottom feeder presents some differences with respect to the geometric model used in these analyses. Indeed, the distribution channels are not immediately attached to the expansion but after a certain length, which will result in a lower ∆p 3D compared with that predicted to our model due to the larger pressure recovery allowed by that configuration [11]. However, the expansion, which as seen is the geometric element that causes the greatest pressure drop, is well represented by the numerical model. Considering that the DCLL uses ceramic walls this estimate can be considered representative of the component. Indeed, even if given the high flow rate places the DCLL in the Q2D turbulence regime, Burr et al [34] have shown that, since the Q2D vortexes are aligned with the magnetic field, they have a negligible Joule and viscous dissipation compared to the pressure drop caused by the IL. Table 6. Total, 2D, 3D pressure drop and equivalent length Leq 3D for all the analyzed cases and estimation at operating magnetic field (OMF) conditions for the configuration A, B and C. The value of k at OMF conditions for the configuration A is estimated with the relation presented in figure 24 and using the data in table 3 and similarly for configurations B and C [22]. The expansion ratio is rex = b i /a and rex = e/a for the configuration C).

Conclusions and future works
In this paper, we have simulated the MHD flow in a prototypical BB bottom feeder made from a perfectly insulating material using a custom-made OpenFOAM solver. The 3D analysis has been focused on Ha = 1000, 1500, 2000 and for geometry with different inlet channel positions. The flow features are strongly influenced by the generation of axial currents in the expansion, which in turn generate a Lorentz force with both a poloidal and an axial component. This significantly impacts the velocity distribution by, for instance, generating high velocity jets at the expansion. The sub-channel walls interact with the upper one by creating a complex phenomenology, dependent on the stiffening plate geometry, that affects the balance of the flow rate among the channels. Generally, the channel aligned with the inlet is favored, not so much for its privileged position, as seen for the hydrodynamic case where the imbalance is present but much lower with respect the MHD cases, but because of the morphology of the upper jet which carries the bulk of the flow. For all the analyzed cases, the channels more distant from the inlet are reached by a flow rate that is very low (<10%), likely insufficient to guarantee the efficient re-circulation of the breeder and acceptable heat transfer coefficients, which calls for counter-measures to be devised. An expansion length could be employed to minimize the three-dimensional effects and the inlet position should be chosen to minimize flow imbalance [17,20]. However, the former strategy could result in a reduction of the available blanket thickness for the BZ [1]. Integration and mechanical stability of the modified manifold are other points that need to be evaluated. Another strategy may rely on the positioning of calibrated orifices at the beginning of the channels, a solution which however increases the pressure drop of the component. Finally, the effect of the conductivity of the wall must be considered which, as shown by Tillack and Morley [15] and Chen et al [35], has a beneficial effect on the distribution of the flow rate, at the expense of an increase in pressure drop. The latter situation will be investigated in future numerical campaigns.
As regards the pressure drop, for all the cases analyzed, more than 90% is due to the 3D losses due to the expansion, the bend and the entrance of the channels. The lower pressure drop occurs for configuration B which has the lowest expansion ratio, with a loss for Ha = 2000 that is approximately 23% lower than configuration C, which has the greatest loss as well as the worst distribution of the flow rate. Configuration B, on the other hand, is also the one that presents the best balance between the courses. Therefore, a configuration with two inlet channels, considering the whole manifold, positioned in such a way as to divide the case into the most equal possible sections is certainly the best performing one and should be preferred if possible in the envelope of interface constraints of the component with the other tokamak systems, i.e. PbLi loop, and without affecting the blanket thermo-mechanical stability.
The presented OpenFOAM solver is planned to be further developed by implementing support for finite wall conductivity (solid-fluid coupling considering a segregated multiregion approach) and a complete MHD heat transfer model (including buoyancy). Furthermore, it is planned to extend the MHD implementation on the interIsoFoam solver, a twophase incompressible, isothermal solver able to simulate highdensity ratio mixtures [32].