Verification and Validation of COMSOL Magnetohydrodynamic Models for Liquid Metal Breeding Blankets Technologies

Abstract

Liquid metal breeding blankets are extensively studied in nuclear fusion. In the main proposed systems, the Water Cooled Lithium Lead (WCLL) and the Dual Coolant Lithium Lead (DCLL), the liquid metal flows under an intense transverse magnetic field, for which a magnetohydrodynamic (MHD) effect is produced. The result is the alteration of all the flow features and the increase in the pressure drops. Although the latter issue can be evaluated with system models, 3D MHD codes are of extreme importance both in the design phase and for safety analyses. To test the reliability of COMSOL Multiphysics for the development of MHD models, a method for verification and validation of magnetohydrodynamic codes is followed. The benchmark problems solved regard steady state, fully developed flows in rectangular ducts, non-isothermal flows, flow in a spatially varying transverse magnetic field and two different unsteady turbulent problems, quasi-two-dimensional MHD turbulent flow and 3D turbulent MHD flow entering a magnetic obstacle. The computed results show good agreement with the reference solutions for all the addressed problems, suggesting that COMSOL can be used as software to study liquid metal MHD problems under the flow regimes typical of fusion power reactors.

1. Introduction

In the nuclear fusion framework, the breeding blanket is a key system devoted to power extraction, shielding and tritium production. Among the different designs of the blanket, in the Water-Cooled Lithium-Lead (WCLL) breeding blanket of DEMO [1] and in the WCLL test blanket module of ITER, the liquid, electrically conducting LiPb is adopted as working fluid to address the above-mentioned functions. The intense magnetic field used in fusion reactors to confine the burning plasma has a strong influence on the flow behavior, producing a magnetohydrodynamic (MHD) effect. In addition, serious temperature gradients are present in the breeding blanket, giving rise to buoyancy forces. The presence of the external magnetic field produces an additional MHD pressure drop, and although rougher MHD studies can predict the ∆p [2], other phenomena, such as turbulence and buoyancy-driven convection, have a drastic impact on blanket performance and require a deeper analysis. In addition, tritium transport mechanisms [3,4,5] are influenced by the magnetic field; therefore, a detailed solution of magnetohydrodynamics is necessary, and it is obtainable with 3D multiphysics models for the breeding blanket. Smolentsev et al. [6] proposed activity for verification and validation of MHD codes for fusion applications, consisting of a series of benchmark problems whose results are known from experimental data or trusted analytical and numerical solutions. In particular, the five problems cover a wide range of magnetohydrodynamic flows which are of interest for fusion applications:

• 2D fully developed laminar steady flow;
• 3D laminar, steady developing flow in a nonuniform magnetic field;
• Buoyancy-driven flow in a square cavity;
• Quasi-two-dimensional (Q2D) turbulent flow;
• 3D turbulent flow.

In this work, a verification and validation procedure of the developed 3D codes is performed, taking as reference [6]. Verification and validation activities of MHD codes were conducted by different authors for 1D codes [7], 3D codes under the COMSOL environment [8] and other codes [9,10,11,12]. Here, in order to further extend the analysis performed by [8], two benchmark cases involving unsteady flows are solved. The Q2D turbulence case proposed by Burr [13] was solved by adopting a modified version of the k − ε model, derived by [14], while the fully 3D turbulent problem was tackled using the Large Eddy Simulation (LES) Residual Based Variational Multiscale (RBVM) method [15,16,17]. For the magnetoconvection case, the problems selected are the ones proposed by Di Piazza and Buhler [18], which are of particular interest for liquid metal breeding blanket technologies. In effect, in the first problem, the non-isothermal condition is due to differentially heated boundaries, while in the second case, temperature gradients are produced by internal heat generation. Both conditions are typical in breeding blanket systems, where strong temperature gradients and intense volumetric heat generation are present.

2. Governing Equations

The flow of liquid metal in breeding blanket conditions is characterized by a small magnetic Reynolds number $R_m=\mu \sigma LU$ (−). Here, $\mu$ (Pa s) and $\sigma$ (S m⁻¹) are, respectively, the dynamic viscosity and the electrical conductivity of the fluid, $L$ (m) and $U$ (m s⁻¹) are the characteristic length and velocity of the flow. The magnetic Reynolds number represents the ratio between induction and diffusion of the magnetic field, so, in the low-$R_m$ approximation, the magnetic field transport can be considered purely diffusive [19]. For an incompressible fluid under an imposed time-independent magnetic field with $R_m$ << 1 the MHD equations, mass conservation Equation (1), momentum conservation Equation (2), current conservation Equation (3) and Ohm’s law Equation (4) are presented.

$$\nabla ·\overrightarrow{u}=0$$

(1)

$$\rho \frac{\partial \overrightarrow{u}}{\partial t}+\rho \left(\overrightarrow{u}·\nabla \right)\overrightarrow{u}=-\nabla p+\mu {\nabla }^2\overrightarrow{u}+\overrightarrow{J}\times \overrightarrow{B}$$

(2)

$$\nabla ·\overrightarrow{J}=0$$

(3)

$$\overrightarrow{J}=\sigma \left(-\nabla \varphi +\overrightarrow{u}\times \overrightarrow{B}\right)$$

(4)

where $\overrightarrow{u}$ (m s⁻¹) is the velocity vector, $p$ (Pa) is the pressure, $\rho$ (kg m⁻³) is the density of the fluid, $\mu$ (Pa s) is the dynamic viscosity, $\overrightarrow{J}$ (A m⁻²) is the current density vector, $\overrightarrow{B}$ (T) is the magnetic flux density, $\sigma$ (S m⁻¹) is the electrical conductivity and $\varphi$ (V m⁻¹) is the electrostatic potential.

For non-isothermal problems, the buoyancy force contribution must be added to momentum conservation Equation (2), which, under the Boussinesq hypothesis, becomes:

$${\rho }_0\frac{\partial \overrightarrow{u}}{\partial t}+{\rho }_0\left(\overrightarrow{u}·\nabla \right)\overrightarrow{u}=-\nabla p+\mu {\nabla }^2\overrightarrow{u}+\overrightarrow{J}\times \overrightarrow{B}+({\rho }_0+\Delta \rho )\overrightarrow{g}$$

(5)

where ${\rho }_0$ (kg m⁻³) is the reference density, $\Delta \rho =\rho -{\rho }_0$ is the density variation with respect to the reference density ${\rho }_0$ and $\overrightarrow{g}$ (m s⁻²) is the gravity vector. $\Delta \rho$ can be further expressed as $-{\rho }_0\beta \left(T-T_0\right)$, where $\beta$ (T⁻¹) is the thermal expansion coefficient and $T_0$ (K) is the reference temperature. The general heat transfer equation must be solved simultaneously:

$$\rho c_p\left(\frac{\partial T}{\partial t}+\overrightarrow{u}·\nabla T\right)+\nabla ·\left(-\lambda \nabla T\right)=\dot{Q}$$

(6)

where $T$ (K) is the temperature, $c_p$ (J kg⁻¹ K⁻¹) is the specific heat at constant pressure, $\lambda$ (W m⁻¹ K⁻¹) is the thermal conductivity and $\dot{Q}$ (W m⁻³) is the volumetric heat generation rate.

The MHD flow is characterized, in addition to the magnetic Reynolds number, by additional nondimensional numbers. Some of the most important is the interaction parameter:

$$N=\frac{\sigma L{B}^2}{\rho U}$$

(7)

that expresses the ratio of the Lorentz force to inertia force and the Hartmann number:

$$Ha=LB{\left(\sigma /\mu \right)}^{1/2}$$

(8)

whose square represents the ratio of the Lorentz force to viscous forces. The Hartmann number and the interaction parameter are related to the Reynolds number, $Re=\rho LU=H{a}^2/N$. For the thermal problems, it is interesting to recall the Grashof number that expresses the square of the ratio between buoyant and viscous forces in the fluid:

$$Gr=\frac{L^3{\rho }^2\beta \Delta Tg}{{\mu }^2}$$

(9)

where $\Delta T$ is a characteristic temperature difference, depending on the problem considered.

3. Verification and Validation of COMSOL Code

The procedure proposed by [6] was followed in order to verify the applicability of COMSOL models, and in the next sections, the benchmark cases are presented, as well as the results of the computations. The first problem is the 2D fully developed MHD flow in a rectangular duct analytically addressed by Shercliff [20] and Hunt [21] (Section 3.1). The second is an experimental case, proposed by Picologlou et al. [22,23,24], involving the study of flows under a fringing magnetic field that investigates the transition from magnetohydrodynamics to ordinary hydrodynamics, presented in Section 3.2. Then, two non-isothermal flow problems are considered, solved numerically by Di Piazza and Bühler [18] (Section 3.3). Lastly, two experimental turbulent flow cases [13,25] are resolved, presented in Section 3.4 and Section 3.5, respectively.

The equations introduced in Section 2 were implemented in COMSOL, exploiting the single-phase flow, heat transfer and electric current modules.

3.1. Two-Dimensional Fully Developed Laminar Steady MHD Flow

The laminar, fully developed, incompressible flow of a conducting fluid driven by a pressure gradient along a rectangular duct under an imposed transverse magnetic field is considered. Shercliff [20] and Hunt [21] solved this problem analytically, using different boundary conditions. Particularly, for Shercliff’s case, the four walls of the duct are nonconducting, while for Hunt’s case, the two walls perpendicular to the magnetic field, called Hartmann walls, are conducting, while the walls parallel to $\overrightarrow{B}$, called side walls, are electrically insulated. The wall conductance ratio:

$$c_w=\frac{{\sigma }_w{t}_w}{\sigma L}$$

(10)

expresses the ratio between the electrical conductivity ${\sigma }_w$ (S m⁻¹) and thickness $t_w$ (m) of the walls and the electrical conductivity of the fluid and the characteristic length of the flow. For Shercliff’s case $c_w=0$ for all the four walls, for Hunt’s case $c_w=0$ for the side walls, while is non-null for the Hartmann walls. As suggested by [6], the wall conductance ratio considered is $c_w=0.01$. Four values of the Hartmann number were selected: $Ha=500, 5000, \mathrm{10,000}, \mathrm{15,000}$.

The problem was solved in dimensionless form, using a hexahedral, structured mesh, analog to the one proposed by Sahu et al. [8], shown in Figure 1. To minimize the computational cost, considering that the solution is symmetric with respect to $x$ and $y$ axis, just a quarter of the domain is considered, applying proper boundary conditions. Elements are generated in $x$ and $y$ direction with a geometric distribution, maximizing the number of cells in the side and Hartmann layers. The mesh, selected following a grid convergence study [5], consists of $50$ elements in the $x$-direction and $75$ in the $y$-direction for $Ha=500$, $64\times 100$ for $Ha=5000$ and $78\times 125$ for $Ha=\mathrm{10,000}$ and $\mathrm{15,000}$; eight cells were always ensured in the Hartmann and side layers. Although the problem is invariant in $z$-direction, one element is generated in the direction of the flow to consider the vector products involving the magnetic field. The velocity boundary conditions adopted in the study are non-slip at the duct walls and periodic flow conditions at the inlet and outlet of the duct. The electrical boundary conditions are electrical insulation for the side walls and thin wall condition $\overrightarrow{J}·\widehat{n}=c_w{\nabla }^2\varphi$, with $\widehat{n}$ unit vector perpendicular to the wall, which represents conservation of electric charge in the plane of the wall.

Solutions obtained are compared to those reported by Smolentsev, choosing the dimensionless flow rate $\tilde{Q}$ as comparison parameter, defined as

$$\tilde{Q}=\frac{4}{b^2}·\frac{\mu u}{\left(-\partial p/\partial x\right)}$$

(11)

where $u$ (m s⁻¹) is the $x$-component of the velocity vector, $b$ (m) is half the Hartmann wall length and $\partial p/\partial x$ (Pa m⁻¹) is the imposed pressure gradient in the direction of the flow. The relative error between COMSOL results and the analytical solutions by Shercliff and Hunt is evaluated as:

$${\epsilon }_{rel}=\left|1-\frac{\tilde{Q}}{{\tilde{Q}}_{an}}\right|$$

(12)

where $\tilde{Q}$ (−) is the COMSOL solution and ${\tilde{Q}}_{an}$ (−) is the analytical solution. The calculations are reported in

Table 1. Comparison between the results of Shercliff’s case and Hunt’s case for 2D fully developed laminar steady flow. Results in terms of $\tilde{Q}$ (−).

Shercliff Case
$\mathit{H}\mathit{a}$(−)	$\mathbf{Analytical}\tilde{\mathit{Q}}$(−)	$\mathbf{COMSOL}\tilde{\mathit{Q}}$(−)	${\mathit{\epsilon }}_{\mathit{r}\mathit{e}\mathit{l}}$(%)
500	$7.680\times {10}^{-3}$	$7.690\times {10}^{-3}$	$0.130$
5000	$7.902\times {10}^{-4}$	$7.906\times {10}^{-4}$	$0.0456$
10,000	$3.965\times {10}^{-4}$	$3.946\times {10}^{-4}$	$0.478$
15,000	$2.648\times {10}^{-4}$	$2.660\times {10}^{-4}$	$0.453$
Hunt Case
$\mathit{H}\mathit{a}$(−)	$\mathbf{Analytical}\tilde{\mathit{Q}}$(−)	$\mathbf{COMSOL}\tilde{\mathit{Q}}$(−)	${\mathit{\epsilon }}_{\mathit{r}\mathit{e}\mathit{l}}$(%)
500	$1.405\times {10}^{-3}$	$1.406\times {10}^{-3}$	$0.0356$
5000	$1.907\times {10}^{-5}$	$1.904\times {10}^{-5}$	$0.184$
10,000	$5.169\times {10}^{-6}$	$5.163\times {10}^{-6}$	$0.118$
15,000	$2.425\times {10}^{-6}$	$2.410\times {10}^{-6}$	$0.635$

. The presented table shows good agreement between analytical and numerical values.

3.2. Three-Dimensional Laminar, Steady Developing MHD Flow in a Nonuniform Magnetic Field

In the second benchmark case, a conducting fluid flows in two different ducts, with rectangular and circular cross sections, in the presence of a nonuniform magnetic field at the exit from a magnet. This case was experimentally investigated at the Argonne National Laboratory on ALEX (Argonne’s Liquid metal EXperiment) facility [22,23,24]. The system employed eutectic NaK as a working fluid in a room temperature closed loop.

In this problem, the magnetic field changes in the direction of the flow $x$, $\overrightarrow{B}=B\left(x\right)\widehat{y}$, with $\widehat{y}$ unit vector in the $y$-direction, and this requires, considering the previously analyzed 2D case, the additional discretization of the domain in the $x$-direction. The velocity boundary conditions adopted in the study are non-slip at the duct walls and imposed average velocity at the inlet. The electrical boundary condition is a thin wall condition on the walls.

3.2.1. Rectangular Duct

The symmetry of the problem is exploited, and only a quarter of the duct cross section is considered. The mesh is constituted by a symmetric distribution of elements in the direction of the flow, maximizing the number of cells in the central region, where the magnetic field is changing the most. In $y$ and $z$ directions, the mesh is analog to the one proposed in Section 3.1. The total number of elements is $2.79\times {10}^5$. The equations are solved in dimensionless form.

The parameters adopted for the study are $Ha=2900$, $N=540$ and $c_w=0.07$. The quantity selected for the comparison with the experimental results is the dimensionless axial pressure difference, that is, the pressure difference developed in the axis of the duct, scaled by $\sigma U{B}_0^2$. The results are presented in Figure 2, where the magnetic field profile scaled by $B_0$ and the axial pressure difference obtained by Picologlou et al. [22] are shown in the present work. Good agreement between the curves can be appreciated. The biggest discrepancy appears in $-5<x/L<0$, where COMSOL tends to overestimate the pressure difference. This behavior is also found in work by Sahu [8] and from the HIMAG Code calculations [26,27]. The difference between the two solutions is calculated using the integral of the curves with the following relation, called integral error index:

$${\epsilon }_{int}=\left|1-\frac{{{\displaystyle \int }}_{x_{min}}^{x_{max}}\mathsf{\Delta }p\left(x\right)\mathrm{d}x}{{{\displaystyle \int }}_{x_{min}}^{x_{max}}\mathsf{\Delta }{p}_A\left(x\right)\mathrm{d}x}\right|$$

(13)

where $\Delta p_A$ (−) and $\Delta p$ (−) are, respectively, the ALEX experiment and the present work nondimensional axial pressure difference. The integrals are computed numerically, using the trapezoidal rule, and the resulting error is $1.10\%$.

3.2.2. Circular Duct

In $x$-direction the mesh adopted is equivalent to the previous case, whereas, in the $y$ and $z$ plane, $25$ boundary layers are considered, generated from the first layer of thickness ${10}^{-6}$ m, with a growth rate of 1.3. The total number of elements is $3.03\times {10}^5$.

The parameters adopted for the study are $Ha=6600$, $N=\mathrm{10,700}$ and $c_w=0.027$. In Figure 3, the results are presented. The curves are matching very well, and the error is $0.913\%$.

3.3. Magneto-Convection

The following benchmark cases were developed with the aim to include also representative cases for the liquid metals breeding blankets, such as the WCLL, being characterized by non-isothermal conditions and internal volumetric heating.

The flow of an electrically conducting fluid in a long vertical channel of the rectangular cross section was considered [16]. The flow is promoted by buoyancy forces arising from non-isothermal conditions; hence, we refer to this as magnetoconvection. With reference to Figure 4, the imposed magnetic field is $\overrightarrow{B}=B_0\widehat{y}$ and the gravitational acceleration $\overrightarrow{g}=-g\widehat{x}$, with $\widehat{x}$ unit vector in $x$-direction, is aligned with the channel axis. Within this frame, two cases are considered: a differentially heated duct and a uniformly heated duct, solved by Di Piazza and Bühler [18] with the CFX commercial code (currently Ansys CFX). For both the problems, the Hartmann number is $Ha=100$.

The problem is 2D, and the COMSOL solution is obtained for the nondimensional problem, expressed in depth in [28]. For the cases considered, $Gr\ll H{a}^4$, from which results that the inertial term of the Navier–Stokes equation became negligible. The mesh adopted is analog to the one shown in Figure 1, where only one element is generated in the direction of $\overrightarrow{g}$. The total number of elements for the mesh selected is $5120$, with $64$ elements in the $x$-direction and $80$ in the $y$-direction; eight cells were always ensured in the Hartmann and side layers. The velocity boundary conditions adopted in the study are non-slip at the duct walls and period flow conditions at the inlet–outlet. The electrical boundary condition is a thin wall condition on the walls. The temperature boundary conditions are defined, for the two different cases, in the next sub-sections.

3.3.1. Differentially Heated Duct

The two boundaries placed in the side walls ($z=-b$ and $z=+b$) are kept at different temperatures, while the Hartmann walls are thermally insulated, and there is no internal heat generation.

A sensitivity analysis, with wall conductance ratio $c_w$ as a changing parameter, was carried out, and the nondimensional velocity profile at $y=0$ as a function of $z$ for half duct is shown in Figure 5 for $Ha=100$.

It is interesting to notice that for the lower values of $c_w$, the damping effect of magnetohydrodynamics is less evident, while in the core region, the solution is still dominated by buoyancy and Lorentz forces and exhibits a linear behavior, with a slope of $\sim Ha$ for the perfectly insulating walls case [29]. In the lower conductivity cases, jets are not present. This is due to the fact that for low values of $c_w$ the side layer becomes better conducting than the side walls, and high current jets are now present in the layers parallel to the side walls. They are also parallel to $\overrightarrow{B}$, so they do not interact with the magnetic field; therefore, the electromagnetic forces in the side layers become negligible, and the dominant effect is due to viscous dissipation.

For these cases, a comparison was made with respect to the numerical solutions [18] obtained with the CFX code. The integral of the curves is selected as a comparison parameter, and the relative difference is calculated with the integral error index:

$${\epsilon }_{int}=\left|1-\frac{{{\displaystyle \int }}_{z_{min}}^{z_{max}}u\left(0,z\right)\mathrm{d}z}{{{\displaystyle \int }}_{z_{min}}^{z_{max}}{u}_B\left(0,z\right)\mathrm{d}z}\right|$$

(14)

where $u_B\left(0,z\right)$ (−) is [18] $x$-direction dimensionless velocity profile and $u\left(0,z\right)$ (−) is the current work profile, both taken at $y=0$. The error comparison is reported in

Table 2. Comparison between COMSOL code and Di Piazza and Buhler solutions for $Ha=100$ . Results in terms of integral error index.

	Differentially Heated Duct	Uniformly Heated Duct
${\mathit{c}}_{\mathit{w}}(-)$	${\mathit{\epsilon }}_{\mathit{r}\mathit{e}\mathit{l}} \left(\mathit{\%}\right)$	${\mathit{\epsilon }}_{\mathit{r}\mathit{e}\mathit{l}} \left(\mathit{\%}\right)$
$0$	$0.957$	$1.50$
$0.01$	$0.326$	$4.78$
$0.1$	$1.77$	$0.585$
$\infty$	$1.36$	$0.770$

. As can be appreciated, for all the cases, the maximum difference is lower than 2%.

3.3.2. Uniformly Heated Duct

For the internally heated duct case, a volumetric heat generation $\dot{Q}$ is present, and the boundary at the side walls is kept at a fixed and equal temperature, while the Hartmann walls are thermally insulated.

The nondimensional velocity profiles at $y=0$ and as a function of $z$-coordinate, a result of an analog sensitivity analysis to the one presented for the differentially heated duct case, are shown in Figure 6 for $Ha=100$. For high wall conductivity ratios, the additional forces damp the velocity profile in the core region, and velocity jets are present in the side layers. For small values of $c_w$, jets are no more present, and the solution at the side layers is dominated by viscous effects.

The error comparison is reported in Table 2, where the maximum difference is related to the analysis with $c_w=0.01$, that, nevertheless, presents a value below 5%. It should be recalled that in both differentially heated duct and uniformly heated duct cases, the comparison was addressed on the numerical solutions of the codes.

3.4. Quasi-Two-Dimensional MHD Turbulent Flow

This case regards a quasi-two-dimensional MHD turbulent flow as proposed in [6]. Burr et al. [13] developed an experimental setup consisting of a rectangular stainless steel channel of side length $0.04$ m and wall thickness $6$ mm where the eutectic sodium–potassium alloy is circulated under the presence of a magnetic field. NaK, with density $865$ kg m⁻³ and kinetic viscosity $9.5\times {10}^{-7}$ m² s⁻¹, flows in the $x$-direction, and $\overrightarrow{B}$ is oriented in $z$ and can be varied from $0.25$ T to $2.5$ T. The electric conductivity of the wall is $1.39\times {10}^6$ S m⁻¹, whereas the one of the NaK is $2.8\times {10}^6$ S m⁻¹, from which it results in wall conductance ratios of the side and Hartmann walls of $c_{w,S}=0.0714$ and $c_{w,H}=0.0119$, respectively. The Hartmann numbers investigated are 600, 1200, 2400 and 4800 for Reynolds numbers between $3.3\times {10}^3$ and $1.0\times {10}^5$.

The problem is solved numerically using a RANS $k-\epsilon$ turbulence model that includes two transport equations for the turbulent kinetic energy $k$ (m² s⁻²) and for the dissipation of turbulent kinetic energy $\epsilon$ (m² s⁻³). For $k-\epsilon$, Equation (2) becomes:

$$\rho \frac{\partial \overrightarrow{u}}{\partial t}+\rho \left(\overrightarrow{u}·\nabla \right)\overrightarrow{u}=-\nabla p+\left(\mu +{\mu }_T\right){\nabla }^2\overrightarrow{u}+\overrightarrow{J}\times {\overrightarrow{B}}_0$$

(15)

where ${\mu }_T$ [Pa s] is the eddy viscosity. The two closure equations of the model are:

$$\left(\overrightarrow{u}·\nabla \right)k=\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]+P_k-\rho \epsilon +S_k^L$$

(16)

$$\rho \left(\overrightarrow{u}·\nabla \right)\epsilon =\nabla ·\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}+S_{\epsilon }^L$$

(17)

Here, $P_k$ (W m⁻³) is a source term, ${\sigma }_k$ (−), ${\sigma }_{\epsilon }$ (−), $C_{\epsilon 1}$ (−) and $C_{\epsilon 2}$ (−) are turbulent model parameters. $S_k^L$ (W m⁻³) and $S_{\epsilon }^L$ (W m⁻³) are source terms that include the damping of the turbulent kinetic energy and the dissipation of the turbulent kinetic energy due to the Lorentz force and were modeled by different authors [30,31,32]. The relations selected are the ones proposed by Meng et al. [14], expressed by the following equations:

$$S_k^L=-\sigma B^{2}k{e}^{-C_1^M\sqrt{\frac{\sigma }{\rho }{B}^2\frac{\nu }{k}}}$$

(18)

$$S_{\epsilon }^L=-\sigma B^2\epsilon e^{-C_1^M\sqrt{\frac{\sigma }{\rho }{B}^2\frac{\nu }{k}}}$$

(19)

where $C_1^M$ (−) is a constant with value $30$, and $\nu$ (m² s⁻¹) is the kinematic viscosity. In these relations, $\sqrt{\sigma /\rho B^2\nu /k}$ is the characteristic turbulence damping time, and $\exp \left(-C_1^M\sqrt{\sigma /\rho B^2\nu /k}\right)$ is the decay rate of the turbulent kinetic energy.

The comparison parameters between Burr and the current study are the mean velocity profiles and turbulent kinetic energy of two-dimensional turbulence $k_{2D}$ (m² s⁻²), defined as

$$k_{2D}=\frac{1}{2}\left(\overline{{u^{\prime }}^2}+\overline{{w^{\prime }}^2}\right)$$

(20)

where $u^{\prime }$ (m s⁻¹) and $w^{\prime }$ (m s⁻¹) are the fluctuating term of the velocity for the $x$-th and $z$-th components.

The mesh refinement was carried out until an appropriate wall lift-off in viscous units ${\delta }_w^{+}$ was obtained. In particular, ${\delta }_w^{+}=11.06$ for every case analyzed, that is the lower limit for COMSOL $k-\epsilon$ turbulence model [33]. This value corresponds to the dimensionless wall distance $y^{+}$, where the viscous sublayer meets the logarithmic layer. The total number of elements is $1.28\times {10}^6$. The velocity boundary conditions adopted in the study are non-slip at the duct walls and imposed average velocity at the inlet. The electrical boundary condition is a thin wall condition on the walls.

The mean velocity in $x$-direction $u$, calculated for $Ha=4800$ and various Reynolds numbers, is now compared with Burr results. In Figure 7, the COMSOL solutions and Burr experimental results are displayed.

The main characteristics of the flow are well expressed, and the influence of the Reynolds number on the flow is evident. This is a characteristic of turbulent MHD flows, while the velocity distribution of laminar MHD flows is governed only by $Ha$. Turbulence smoothens out velocity peaks in the side walls that are reduced for increasing Reynolds numbers, and the width of the side layer increases with $Re$ due to turbulent transfer of momentum.

In

Table 3. Local error for the quasi-two-dimensional turbulent flow for $Ha=4800$ and different $Re$. Results in terms of $u (-)$.

$\mathit{R}\mathit{e}$(−)	$\mathbf{Experimental}\mathit{u}$(−)	$\mathbf{COMSOL}\mathit{u}$(−)	${\mathit{\epsilon }}_{\mathit{r}\mathit{e}\mathit{l}}$(%)
$3.3·{10}^3$	$1.638$	$1.644$	$0.368\%$
$3·{10}^4$	$1.547$	$1.492$	$3.52\%$
$6·{10}^4$	$1.442$	$1.482$	$3.31\%$
$1·{10}^5$	$1.405$	$1.447$	$2.98\%$

, the local relative errors between COMSOL and the experimental results, calculated with the following equation, are presented.

$${\epsilon }_{rel}=\left|1-u/u_{exp}\right|$$

(21)

Here, $u$ (−) is the $x$-direction velocity magnitude calculated with COMSOL code and $u_{exp}$ (−) by the experiment, both evaluated at $z=0.45$, that is, the closest point to the side layer, obtained in the experiment. As it can be observed, values agree very well, with a maximum error of about $3.52\%$. Comparing the velocity at $z=0$, there is a relative difference of about 12% between the numerical and experimental values. The code overestimated the bulk velocity, and the same behavior is reported in [14], which uses the same strategy to model the dissipation of the turbulent kinetic energy due to Lorentz forces.

In Figure 8, the comparison between the distributions of the turbulent kinetic energy of two-dimensional turbulence $k_{2D}$ reported by Burr [13] and calculated with COMSOL are presented for $Re=1.0\times {10}^5$ and Hartmann numbers between $600$ and $4800$. The values are captured along the $z$-axis at the midplane $y=0$. The increase in the turbulent kinetic energy as $Ha$ increases can be appreciated, proving that turbulence is promoted by the magnetic field. In both [14] and COMSOL results, turbulence is restrained to the side layers, decreasing fast moving towards the core region, where, in the experimental results, the flow, although weakly, remains turbulent. The anisotropicity of the turbulent flow is particularly evident for the high Hartmann number cases, where $k_{2D}\approx k_{3D}$.

As shown by the results, the code can tackle quasi-two-dimensional MHD flow problems, giving reliable results, particularly in the side layer region. Further improvements are needed to better compute the bulk turbulence that is underestimated by the code.

3.5. Three-Dimensional Turbulent MHD Flow

The benchmark problem on 3D turbulent MHD flow addressed is the one proposed by Smolentsev et al. [6]. The eutectic GaInSn, with density $6360$ kg m⁻³, the electrical conductivity of $3.46\times {10}^{-6}$ S² m⁻¹ and kinetic viscosity of $3.4\times {10}^{-7}$ m² s⁻¹, flows with a maximum flow rate of $2\times {10}^{-3}$ m⁻³ s⁻¹ in plexiglass (insulating) rectangular channel of length $0.5$ m and $100$ mm × $20$ mm cross section, and starting from pure hydrodynamics conditions, is subjected to a nonuniform magnetic field generated by a magnetic obstacle. This is an experimental problem addressed by Andreev et al. [25]. The flow direction is $x$-oriented, and the magnetic field is in the $z$-direction. Starting from a zero value, the magnetic field monotonically increases until it reaches the maximum value of $B_0=0.504 \mathrm{T}$ at the center of the duct, corresponding to $Ha=400$, then it decreases to zero. The magnetic field is slightly nonuniform also in the $y$-direction, but this feature is neglected in the COMSOL model. Further information on the B profile can be found in [25]. The Reynolds number selected is $Re=4000$; therefore, the interaction parameter is equal to $N=40$.

The problem was solved using the Large Eddy Simulation (LES) model [34,35], as suggested by Smolentsev et al. [6]. In particular, the Residual Based Multiscale Variational (RBMV) method [16,34,35] was implemented. In this model, the velocity and pressure fields are decomposed into resolved and unresolved scales:

$$\overrightarrow{U}=\overrightarrow{u}+\overrightarrow{u^{\prime }}$$

(22)

$$P=p+p^{\prime }$$

(23)

where $\overrightarrow{u}$ (m s⁻¹) and $\overrightarrow{u^{\prime }}$ (m s⁻¹) are the resolved scale and the unresolved scale velocities, respectively, and $p$ (Pa) and $p^{\prime }$ (Pa) are the resolved scale and the unresolved scale pressures. In the RBMV method, the unresolved velocity and pressure scales are modeled in terms of the equation residuals for the resolved scales. Further information on RBMV LES modeling can be found [33].

To ensure adequate space discretization, the resolution of wall layers was checked by $\frac{u_{\tau }{h}_w}{\nu }<1$, where the left term is the dimensionless wall distance evaluated at the first mesh cell next to the wall. Here, $u_{\tau }$ is the friction velocity, $h_w$ is the thickness of the first mesh cell next to the wall,$\nu$ is the kinematic viscosity. Time discretization was checked using the relation $C=\frac{U\mathsf{\Delta }t}{h_U}<0.5$, where $C$ is the Courant number, with $U$ flow velocity magnitude, $\mathsf{\Delta }t$ time step and $h_U$ mesh size in the streamline direction. The total number of elements of the selected mesh is $1.8\times {10}^6$. The velocity boundary conditions adopted in the study are non-slip at the duct walls and imposed average velocity at the inlet. The electrical boundary condition is a thin wall condition on the walls.

The selected comparison parameter is the mean velocity profile and the mean electric potential, evaluated at different distances along the channel. In particular, the selected locations correspond to the main flow regions, as described by Andreev. In Figure 9, a comparison between COMSOL and Andreev’s velocity profile at $\frac{x}{H}=-5.3$, with $H$ channel height in $z$-direction, is presented.

The profile is referred to as the first region indicated by the author, characterized by the increasing magnetic field that influences the flow and damps its perturbations, called “turbulence suppression region”. As evident, the agreement between the curves is good, and the integral error index is $0.947\%$. In Figure 10, the electric potential profile, placed at $\frac{x}{H}=0$, is shown.

This region, around the center of the duct where $B$ is maximum, is called “vortical region”. The magnetic field suppresses the fluctuations in the direction parallel to the magnetic field, and the flow becomes quasi-two-dimensional. The integral error index is $4.21\%$. In Figure 11 and Figure 12, the results for the last region, named “wall jet region”, are reported.

This region is located on the remaining part of the channel, where the magnetic field decreases. As observable, the velocity in the middle plane increases greatly, going from $x/H=3$ to $x/H=6$ thanks to the drop of the intensity of the magnetic field, and the quasi-two-dimensional profile stretches to the core of the flow. The agreement between the velocity fields is quite good, and the relative errors in percentage are $6.72\%$ at $x/H=3$ and $7.81\%$ at $x/H=6$.

As presented, the code is capable of representing the characteristic regions of the experimental problem, and the errors obtained are, for every case, well below $10\%$. The results are summarized in

Table 4. Comparison between COMSOL code and [25] experiment. Results in terms of integral error index.

$\mathit{x}/\mathit{H}$(−)	$-5.3$	$0$	$3$	6
${\epsilon }_{rel}$ (%)	$0.947\%$	$4.21\%$	$6.72\%$	$7.81\%$

, and they provide confidence in the capabilities of COMSOL to simulate fully 3D turbulent flows.

4. Discussion

In this paper, a verification and validation procedure was followed, as proposed by Smolentsev et al. [6], and different liquid metal MHD problems were solved, with the aim of verifying the developed models using the commercial software COMSOL Multiphysics. For the magnetoconvection case, the considered benchmarks were the ones tackled by Di Piazza and Bühler [18], presenting the typical conditions that are expected in liquid metal breeding blankets.

The compared parameters showed great agreement for laminar flow problems, both isothermal and non-isothermal. As far as the turbulent cases are concerned, Q2D and 3D, a modified version of the RANS $k-\epsilon$ and the LES RBMV models were adopted, respectively, and the deriving results, in terms of relative errors and the capability of describing the flow features, are very promising.