A Formal Verification and Validation of a Low Magnetic Reynolds Number MHD Code for Fusion Applications

IEEE Proof IEEE TRANSACTIONS ON PLASMA SCIENCE 1 A Formal Verification and Validation of a Low Magnetic Reynolds Number MHD Code for Fusion Applications D. Suarez , A. Khodak , Member, IEEE, E. Mas de les Valls, and L. Batet Abstract— As the nuclear fusion research advances, researAQ:1 1 chers and engineers focus more on the design of the required2 systems that complement the nuclear fusion reaction in the3 plasma of a Tokamak. Some proposals for breeding blankets as4 well as plasma-facing components’ protection systems are based5 in liquid metal flows under the Tokamak intense magnetic fields.6 This creates the situation where induced magnetic field can be7 neglected and the low magnetic Reynolds number (Re) electric8 potential formulation can be used to close the magnetohydro-9 dynamic (MHD) system of equations. In the last few years,10 many different laboratories have developed their own MHD11 codes to study the liquid metal flow. A formal verification and12 validation of such codes is necessary to enhance the reliability13 of the numerical results and to make sure that design decisions14 are based on safe grounds. The fusion community has made15 the effort of proposing standardized benchmark cases through16 which any MHD code should demonstrate its reliability. This17 work contains the formal validation and verification activities of18 the MHD code developed some years ago in UPC and currentlyAQ:2 19 candidate to contribute to the research done at PPPL. The code20 is implemented over OpenFOAM which makes it easily modifi-21 able. Among these benchmark cases, there are high Hartmann22 number (Ha), 3-D flows, and magneto-convective interaction23 cases.24 Index Terms— Computational fluid dynamics (CFD), liquid25 metal flow, magnetohydrodynamics (MHDs), verification and26 validation.27 I. INTRODUCTION28 LIQUID metal flows under the intense magnetic fields29 generated in the Tokamak can be encountered in some30 design proposals such as the divertor (or in general the plasma31 facing components) film flow protection system [1] or inside32 the breeding blankets, such as the dual-coolant lead-lithium33 (DCLL), the water-cooled lead-lithium (WCLL) or the helium-34 cooled lead-lithium (HCLL) blanket [2]. In all of them, the35 AQ:3 Manuscript received 4 February 2022; revised 30 June 2022; accepted 23 August 2022. This work was supported in part by the European Union through the Euratom Research and Training Programme within the framework of the EUROfusion Consortium under Grant 101052200—EUROfusion and in part by the Associació/Col·legi d’Enginyers Industrials de Catalunya with Fundació Caixa d’Enginyers. The review of this article was arranged by Senior Editor G. H. Neilson. (Corresponding author: D. Suarez.) AQ:4 D. Suarez, E. Mas de les Valls, and L. Batet are with the AQ:5 Universitat Politecnica de Catalunya, 08034 Barcelona, Spain (e-mail: [email protected]). A. Khodak is with the Princeton Plasma Physics Laboratory, Princeton, NJ 08540 USA. Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPS.2022.3203801. Digital Object Identifier 10.1109/TPS.2022.3203801 liquid behavior is that of a magnetohydrodynamic (MHD) 36 flow, since the fluid is electrically conducting and electric 37 current will be induced in it. The electric current under the 38 magnetic field will generate the so-called Lorentz forces, that 39 will strongly affect the flow structure. 40 The highly intense external magnetic field is much higher 41 than the induced magnetic field intensity predicted in the 42 above-mentioned applications for the liquid metal flows. 43 Such condition is known as the low magnetic Reynolds 44 number (Re) [3]. 45 Even though liquid metal MHD flows can be found in tran46 sient cases, the characteristic frequency of these problems is 47 small enough to neglect the transient terms in electric potential 48 equation; i.e., the so-called quasi-static approximation can be 49 used for liquid metal MHD flows. 50 Computational simulations provide an important contribu51 tion to understand the flow characteristics of the liquid metal 52 under such circumstances. The equations that govern the MHD 53 phenomenon must be implemented in a computational fluid 54 dynamics (CFD) code. The code must later be verificated and 55 validated (V&V) as promoted by AIAA in 2002 [4] for general 56 purpose CFD codes, to provide reliability to the results that 57 may outcome after its use. 58 The fusion community, in 2015, agreed to perform such 59 process for liquid metal MHD codes over the selected bench60 mark cases [5]. In such document, five cases were proposed to 61 cover a wide range of MHD flows, from laminar to turbulent 62 flows, which are of interest for fusion applications. 63 In 2020, a code-to-code comparison of different 64 magneto-convective MHD codes was performed in [6]. 65 This exercise is foreseen as a step forward in the capability 66 to predict this type of flows, specially of interest for breeding 67 blanket applications. 68 Several research groups have gone through the V&V 69 process of their codes, following the steps proposed by 70 Smolentsev et al. [5]. Some works that present parts of this 71 effort are mentioned in the following lines. 72 In 2017, Khodak [7], in PPPL, described the implementation 73 of the MHD governing equations into the general purpose CFD 74 code CFX, part of the ANSYS Workbench. Two-dimensional 75 Shercliff and Hunt cases were validated as well as the 3-D 76 fringing magnetic field case. 77 In 2018, Patel et al. [8], from UCLA, validated the codes 78 HIMAG and FLUIDYN for Shercliff and Hunt cases, and later 79 applied them to a 3-D case with some 90◦bends. 80 0093-3813 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. IEEE Proof 2 IEEE TRANSACTIONS ON PLASMA SCIENCE In 2018, Gajbhiye et al. [9], in IIT Hyderabad, validated81 their code with the Shercliff, Hunt, magneto-convective flow,82 and 3-D fringing magnetic field. They studied also a curved83 pipe MHD flow, with very good accuracy.84 In 2018, Sahu and Bhattacharyay [10], from IPR, validated85 the COMSOL code for laminar 2-D cases (Shercliff and Hunt),86 3-D fringing magnetic field, buoyant MHD flow and eventually87 studied a conduit with several 90◦bends.88 In 2021, Blischchik et al. [11], presented a large set of89 benchmark cases of MHD flows to validate the solver imple-90 mented over OpenFOAM in TUDelft. They present the cases91 of Shercliff, Hunt, backward facing step, multiphase cavity92 flow, rising bubble in liquid metal, and turbulent conjugate93 MHD flow.94 In 2021, Alberghi et al. [12], from the Polytecnico di Turin,95 have done a step further in the V&V of the COMSOL code,96 working out the application of turbulence models for the97 solution of the MHD turbulent and quasi-2-D-turbulent flows.98 In some cases, either in the works listed here or else-99 where, the code has been validated with benchmark cases100 that contain an interesting specific feature of the phenomena101 under investigation. Commonly, these types of add-hoc vali-102 dations have not followed the formal procedure proposed by103 Smolentsev et al. [5].104 The work presented here is focused on the formal verifica-105 tion and validation of the MHD code developed at UPC by106 Mas de les Valls in her Ph.D. thesis [13], and used in previous107 works such as [14] for multiple regions, or [15] for magneto-108 convective cases, among others.109 II. CODE DESCRIPTION110 The set of equations implemented in the model are the111 Navier-Stokes equations coupled to the relevant set of reduced112 quasi-static Maxwell equations. The flow is considered incom-113 pressible. The main electromagnetic variable is the electric114 potential, which is solved using the Poisson’s equation derived115 from the application of the divergence of the Ohm’s law.116 The electric potential formulation suits perfectly with the117 low magnetic Re of the liquid metal flows considered in the118 fusion reactors. The electric current density is treated using119 the conservative scheme proposed by Ni et al. [16]. The set120 of equations solved are121 ∂ρ ∂t+∇·(ρv)=0 (1) 122 ∂ρv ∂t+∇·(ρvv)=−∇ p+∇·(τ)+j×B(2) 123 ∇·(σm∇φ)=∇·σm(v×B)(3) 124 j=σm(−∇φ+v×B)(4) 125 where vis fluid velocity vector, ttime, ppressure, ρdensity126 of the fluid, τthe stress tensor, jelectric current density vector,127 φelectric potential, σmelectric conductivity of the fluid, and128 Bexternal magnetic field vector. The dimensionless numbers129 that describe the ratio of forces are: the Hartmann number130 (Ha =BL(σm/μ)1/2), the Re (Re =UL/ν), and the Stuart131 number or the interaction parameter (N=Ha2/Re). It is noted132 that Band Ucorrespond to the magnetic field magnitude133 and the mean velocity magnitude. The parameter Lis the134 characteristic length of the problem, that depends on the case 135 and ν=μ/ρis the kinematic viscosity. For those cases in 136 which the surrounding wall of the liquid metal is electrically 137 conducting, (3) and (4) are solved in the solid region. At the 138 boundaries, the conservation of jis guaranteed. An important 139 dimensionless number that defines the effect of the wall in the 140 flow behavior is the wall conductance ratio (cw=σwtw/σL). 141 The solver uses the PISO approach to resolve the pressure142 velocity coupling, which is the usual procedure for incom143 pressible flows [17]. The equations have been implemented 144 and solved on the OpenFOAM suite [18], that is based on 145 the finite-volume method approach, which guarantees the 146 volumetric conservation of the calculated variables. Unless 147 otherwise specified, the interpolation method used in this 148 work between values at cell centers to cell faces and vice 149 versa is linear interpolation. Similarly, second-order upwind 150 discretization scheme has been used for advective terms. 151 III. BENCHMARK CASES 152 The selected benchmark cases for the validation process 153 have been: problems A1 and A2 (fully developed laminar 154 steady MHD flow), problem B2 (3-D laminar steady MHD 155 flow), and problem E (MHD flow with heat transfer) as defined 156 in [5]; and the code-to-code comparison of a 3-D laminar 157 unsteady MHD flow with heat transfer [6]. 158 A. Problem A (Fully Developed Laminar Steady MHD Flow) 159 The proposed benchmark for the 2-D fully developed lam160 inar MHD flow is the classical problem of a conducting fluid 161 flowing in a rectangular channel with a transvers uniform 162 magnetic field. To proof the capability of the code to deal 163 with insulating walls and with electrically conducting walls, 164 two subcases of problem A were proposed: classical Shercliff’s 165 case [19] (problem A1, for insulated walls) and Hunt’s 166 case [20] (problem A2, for conducting Hartmann walls). 167 In both subcases the characteristic length of the problem L168 is half width of the channel in the direction of the magnetic 169 field. 170 1) A1 (Shercliff Case): In the classical Shercliff case, the 171 formation of two types of boundary layers can be clearly 172 observed. The electric currents (j) induced in the fluid will be 173 aligned perpendicular to the magnetic field (B), and will there174 fore produce the Lorentz force, which will oppose to the flow 175 of the fluid in the bulk central region of the channel. At the 176 edges of the rectangular channel, near the walls parallel to the 177 magnetic field, the electric current will recirculate and align 178 with the magnetic field. In this region, the Lorentz force will be 179 negligible, since jand Bwill be mostly parallel. This boundary 180 layer is called side layer, and its thickness is approximately 181 ∼L/(Ha)1/2. When the electric current reaches the corner, the 182 jdirection will change again, to be perpendicular to B, but in 183 this case in the opposite direction than in the bulk region. This 184 happens in the wall perpendicular to the magnetic field, and 185 creates the so-called Hartmann boundary layer, of thickness 186 ∼L/Ha, where the Lorentz force pushes the fluid. 187 2) A2 (Hunt Case): In the original work by Hunt, he derived 188 mathematically the dimensionless flow rate ˜ Qand the velocity 189 profile for an arbitrary conducting thin wall next to the 190 IEEE Proof SUAREZ et al.: FORMAL VERIFICATION AND VALIDATION OF A LOW MAGNETIC Re MHD CODE 3 TABLE I A. FULLY DEVELOPED LAMINAR STEADY MHD FLOW Hartmann boundary layer, while the wall next to the side191 layer is electrically insulated. In this case, the electric current192 density in the side layers is, like in the Shercliff case, aligned193 with the magnetic field, and the Lorentz force there is also194 negligible. The electric current in the Hartmann layer, however195 is completely different. The electric current density lines196 encounter a new path out of the fluid domain, through the197 electrically conducting thin wall. This prevents the returning j198 from proving the thrust that is found in the Hartmann boundary199 layer of the Shercliff case. The flow distribution in such cases200 is the known as M-shape velocity profile, where the velocity201 jets appear in the side boundary layers, where the Lorentz202 force is negligible.203 The results of the Shercliff and Hunt cases (shown in204 Table I) have been compared with the analytical solutions205 provided in [19] and [20] for the dimensionless flow rate206 ˜ Q=�1 −1�1 −1˜ Ud ˜yd ˜z, where ˜z=z/aand ˜y=y/L. In this 207 case, the direction of Bis aligned with y. The dimensionless208 velocity is ˜ U=U/[L2ν−1ρ−1(−dp/dx)], where −dp/dx is209 the uniform pressure gradient in the flow direction. Selected210 discretization schemes were central differencing, of secondAQ:6 211 order and the steady state for A cases was considered reached212 when (˜ Ui−˜ Ui−1)/( ˜ Ui−1) < 10−10. A grid sensitivity study213 was carried out for those cases with Ha =500. It is noted214 that for Shercliff case it is clear that better resolution results in215 better agreement with the analytical data, but the improvement216 in the Hunt case is below the decimal digits on the table.217 Profiles of Hartmann and side boundary layers for the most218 challenging cases (Ha =15 000) are shown in Fig. 1 (noncon-219 ducting walls) and Fig. 2 (conducting Hartmann walls). The220 agreement of the obtained result with the analytical profile221 confirms the accuracy of the calculated ˜ Q. In the figure,222 one can see that the boundary layers extend further than the223 estimated distance ∼L/(Ha)1/2(for side layers) and ∼L/Ha224 (for Hartmann layers). In Fig. 1, it can be seen that in A1225 cases the flow is carried by the central flat region. In Fig. 2,226 it can be appreciated that clearly, in the A2 case, most of the227 flow is carried by the side boundary layers jets. At the end of228 the jet, next of the bulk region, velocity shows negative values,229 as predicted by the analytical solution.230 Fig. 1. Hartmann boundary layer (top) and side boundary layer (bottom) of nonconducting walls case with Ha =15 000. Analytical data from [20] (problem A1). For A1 and A2 cases, six cells have been placed both 231 inside the Hartmann boundary layer distance (∼L/Ha) and 232 in the side boundary layer distance (∼L/(Ha)1/2). Although AQ:7 233 side boundary layer is larger, the code accuracy is sensible to 234 the cells located in the Hartmann boundary layer. Table I also 235 includes the mesh details for each case. The number of cells 236 shown correspond to a fully squared channel; however thanks 237 to the symmetry of the case, the actual simulation considered 238 half number of cells in the direction of the magnetic field. 239 B. Problem B (3-D Laminar Steady MHD Flow) 240 The selected benchmark case for 3-D laminar steady MHD 241 flow was the nonuniform (fringing) transverse magnetic field at 242 the exit of a magnet. Two geometries for the straight conduits 243 were proposed. The first one (B1 problem) is a circular pipe 244 flow and the second one (B2 problem) is a rectangular chan245 nel flow. Experimental results for both cases were reported 246 together in the work by Reed et al. [21]. In these problems, 247 the flow of the liquid metal through a fringing transverse 248 magnetic field induces electric currents in the direction of the 249 IEEE Proof 4 IEEE TRANSACTIONS ON PLASMA SCIENCE Fig. 2. Hartmann boundary layer (top) and side boundary layer (bottom) of conducting Hartmann walls case with Ha =15 000. Analytical data from [20] (problem A2). TABLE II B. DIMENSIONLESS PARAMETERS FOR THE BENCHMARK PROBLEM B2 TABLE III B. COEFFICIENTS TO FIT (5) fluid flow, which represents the main difference from the A250 problems. AQ:8 251 The present work has focused on the validation of252 problem B2, where the channel is rectangular and perpendic-253 ular to the magnetic field. The fluid flows from left to right254 through the rectangular channel in the MHD region as fully255 developed (−15 <x/L<−5), it encounters the variation256 of the transverse magnetic field around −5<x/L<5,257 and eventually becomes hydrodynamic (non-MHD) flow (5 <258 x/L). The flow characteristics (dimensionless variables Ha,259 cw, and N) are shown in Table II. The dimensionless magnetic260 field ˜ B=B/B0has been fit to the reference work field261 through a function (5) [7], where the coefficients are indicated262 in Table III263 ˜ B(x/L)=1 2·�1−tanh�(x/L−w)·a+(x/L−w)2·b264 +(x/L−w)3·c��.(5)265 No turbulence models have been applied although the high266 Re number for the hydrodynamic region (Re ∼15 000). The267 mesh used for this case was 460 ×182 ×142 in the flow268 direction x, the direction of the magnetic field (y), and the269 z-direction, respectively. The mesh allocated five cells in the270 Hartmann boundary layer and 15 in the side layer. To solve this271 steady state case, a relaxation factor was applied to velocity272 field and upwind discretization was applied for the advective273 term to smooth the pressure difference results.274 Fig. 3. Cross-sectional pressure difference and magnetic field distribution along the rectangular channel. Benchmark data from [21] (problem B2). Two results achieved in this work are shown in Fig. 3, 275 together with the benchmark data. The first one is the result 276 of applying the magnetic field described in (5), without con277 sidering any correction regarding the curl-free nature of any 278 magnetic field. The field is therefore 1-D, and the rotational 279 ∇×B=0; this case is named “inconsistent.” The second 280 result shown corresponds to the simulation where the external 281 magnetic field has been corrected to guarantee the curl-free 282 (∇×B=0) and the solenoidal (∇·B=0) nature of the 283 applied fringing magnetic field. The corrected field is then 284 2-D, and is named “consistent.” 285 The followed procedure to achieve a consistent magnetic 286 field has been obtained from the very detailed study carried out 287 by Albets-Chico et al. [22]. This study analyzes the influence 288 of the consistency of the fringing magnetic field on the numeri289 cal simulation results. More specifically, we have followed two 290 correction levels (iterations), as considered enough by authors. 291 The numerical results shown in Fig. 3 emphasize the impor292 tance of providing a consistent magnetic field to the model, 293 since the consistent field reaches much better agreement with 294 the experimental data than the inconsistent field, which yields a 295 pressure difference peak value 15% smaller than the consistent 296 field result. The code can smoothly reproduce the laminar 297 steady 3-D phenomena, observed by the magnitude of the 298 pressure difference between Hartmann and side boundary 299 layers. 300 Albets-Chico et al. [22] also state that among the contribu301 tors to the disagreements observed by the community between 302 experimental and numerical results are the fitting method for 303 the magnetic field. During the development of this work, 304 in a preliminary simulation, the fitting function described 305 in [23] was used. The results were much different than the 306 expected, showing a much smaller and wider peak. Using 307 fitting method (5) provided much better results. 308 A different but also important issue related to this simulation 309 is the width of the electrically conducting wall. Although pre310 vious validation of other code [7] used 3-D formulation of the 311 magnetic field with ∇×B=0, some problems on the accuracy 312 of the calculated pressure difference were reported. The cited 313 work [7] suggested that the information regarding the wall 314 thickness (that was not available in the original experimental 315 IEEE Proof SUAREZ et al.: FORMAL VERIFICATION AND VALIDATION OF A LOW MAGNETIC Re MHD CODE 5 work [21]) had a strong influence on the magnitude of the316 peak. Thinner walls provided a larger pressure difference peak.317 The best agreement with the experiment was obtained when318 the ratio between liquid conductivity and wall conductivity319 was 10. This ratio was used in our simulation.320 A different interesting topic might require attention to the321 results of this simulation. The code validation performed by322 Khodak [7] showed a relevant nonzero pressure difference323 between boundary layers in the MHD region that our sim-324 ulation did not capture. Our results show, indeed, pressure325 differences in the fully developed MHD region but they are326 small enough to be negligible in Fig. 3.327 C. Problem E (MHD Flow With Heat Transfer)328 The selected benchmark for the MHD buoyancy-driven con-329 vection problem is based on the experimental and theoretical330 study described in [24].331 The liquid metal flows in natural convection inside a vertical332 enclosure of square cross section (40 ×40 ×300 mm3)333 where two vertical walls have different temperatures and the334 others are adiabatic. The magnetic field is perpendicular to the335 temperature gradient.336 The momentum equation has been adapted to include the337 buoyant term modeled using the Boussinesq hypothesis338 ∂ρv ∂t+∇·(ρvv)=−∇ p+∇·(τ)+j×B+ρβθg(6) 339 where βis the thermal expansion coefficient, gis gravity340 vector, and θis the difference between the cell temperature341 and the mean temperature of the enclosure. The dimensionless342 parameters that describe the flow behavior are the Ha and the343 Grashof number (Gr =gβT(2a)3/ν2). ais half width of the344 channel in the heat flux radial direction. It is noted that for345 this problem, the characteristic length for Gr number and for346 Ha number is 2a, so that Ha =B(2a)(σm/μ)1/2. Two series347 of simulations have been run, one for Gr number equal to348 4·107to compare to the numerical results and another one349 for Gr =3·107to compare with the experimental results,350 provided in [24].351 The enclosure walls are electrically insulating so only one352 region has been modeled.353 For all cases, the mesh used has 250 ×84 ×60 cells,354 in the vertical direction, magnetic field direction, and thermal355 gradient direction, respectively.356 The main goal of reproducing problem E is to obtain the357 Nusselt number (Nu) for different Ha numbers. Experimental358 and numerical results from the original work [24] show a359 maximum of the Nu around Ha =200, although the values360 were not the same for experimental and numerical analysis.361 The definition of the Nu number is362 Nu =�Sqds kLM(Th−Tc)2a(7) 363 where Sis the surface of the hot wall, q is the heat364 transferred from the wall to the liquid metal, kLM is the thermal365 conductivity, and Thand Tcare the temperatures of hot and366 cold walls.367 Fig. 4 shows the instantaneous velocity magnitude dis-368 tribution in the insulated vertical enclosure of problem E369 Fig. 4. Instantaneous velocity magnitude distribution in the mid-plane of the vertical enclosure for Ha =0, 100, 200, and 900 for Gr =4×107 (problem E). Fig. 5. Transient evolution of the Nu for different Ha of 3-D calculation of MHD flow with heat transfer (Gr =4×107) in a finite enclosure (problem E). (Gr =4×107). The evolution of the Nu for different Ha 370 numbers can be seen in Fig. 5. These results show very good 371 agreement with the original work [24] numerical simulations 372 (Gr =4×107) and capture well the magneto-convective 373 phenomena. The mean and maximum/minimum values shown 374 in Fig. 6 have been computed for the second half of the 375 time range shown in Fig. 5 (0.025 <˜ t<0.05), where 376 ˜ t=t/(L2×ν). They agree very well with those presented 377 in the original work. A very good agreement is also observed 378 IEEE Proof 6 IEEE TRANSACTIONS ON PLASMA SCIENCE Fig. 6. Comparison of Nu for different Ha of MHD flow with heat transfer in a finite enclosure for different Gr. Benchmark data from [24] (problem E). Fig. 7. Case results. Time-averaged velocity and temperature distributions and time-averaged velocity profiles along the duct. Velocity profiles are slightly rotated (code-to-code comparison). in the Ha where the maximum Nu can been found, Ha ∼200.379 The experimental results of the original work show a lower Nu380 number for all Ha numbers (Gr =3×107), not only compared381 to our calculation but to their own calculations.382 D. 3-D Laminar Unsteady MHD Flow With Heat Transfer383 (Code to Code Comparison)384 Finally, the code-to-code comparison work among different385 liquid metal MHD codes [6] has been chosen as the latest386 validation case.387 The selected case is especially interesting since it contains388 several features of the liquid metal flow in blankets. The389 characteristic dimensionless numbers for this benchmark case390 are Ha =220, Re =3040, and Gr =2.88 ×107. The selected391 mesh for the bulk region of this case was 631 ×76 ×102,392 containing six cells in the Hartmann boundary layer and393 18 cells in the side boundary layer.394 Fig. 8. Case results. Mean velocity and mean temperature profiles at different axial locations. HIMAG data from [6] (code-to-code comparison). Fig. 9. Case results. Time-averaged temperature at the duct axis. HIMAG data from [6] (code-to-code comparison). The case considers a liquid metal flowing in the direction 395 of gravity through a rectangular channel with electrically 396 conducting walls. The inlet velocity is flat and the magnetic 397 field at the entrance is zero. In the first region, the flow 398 develops the boundary layer and the increasing magnetic field 399 forces the generation of the characteristic side layer jets. After 400 this first development region, the flow finds a wall heated 401 atafixed heat flux, that triggers buoyant effect in the flow 402 and some recirculation appears. After the heated region, the 403 magnetic field decreases and the flow becomes hydrodynamic 404 again, carrying still the effects of the buoyancy. 405 Among the data to report after the execution of the bench406 mark case, there is the time-averaged velocity and temperature 407 profiles at the midplane for different axial locations. Time 408 averaging has been calculated from 300 to 600 s. These results 409 can be observed in Figs. 7 and 8. The time-averaged tem410 perature evolution along the axial direction can be observed 411 IEEE Proof SUAREZ et al.: FORMAL VERIFICATION AND VALIDATION OF A LOW MAGNETIC Re MHD CODE 7 Fig. 10. Case results. Velocity and temperature fluctuations at (0, 0, 0). HIMAG probe was located at (0, −0.00575, 0.00575). HIMAG data from [6] (code-to-code comparison). in Fig. 9. The velocity and temperature fluctuations on the412 central point of the domain can be observed in Fig. 10. The413 results of our calculations were compared with those obtained414 by HIMAG code. The reader should note that the probe was415 located in different place in HIMAG than in our calculation.416 IV. CONCLUSION417 In this work, the low magnetic Reynolds MHD code for418 liquid metals developed at UPC [13] has been validated in419 different conditions.420 Most of the validation and verification standard cases rec-421 ommended in [5] have been reproduced including: the fully422 developed 2-D MHD flow (problem A), 3-D laminar steady423 MHD flow (problem B2), and MHD flow with heat transfer424 (problem E).425 To complement the standard validation cases, the code-426 to-code comparison exercise [6] based on 3-D laminar427 unsteady MHD flow with heat transfer has been successfully428 reproduced.429 The code has proved to provide excellent agreement with430 the benchmark cases and suits well in MHD flows in channels431 and containers.432 Some more research must be done in improving the solver433 performance in round pipes (B1 case, associated with the434 correct discretization procedures of non-orthogonal meshes).435 ACKNOWLEDGMENT436 The authors would like to thank the contribution of Dr.437 Smolentsev providing the data for code comparison. They438 would also like to thank EUROfusion for the allocation of439 HPC capacity on Marconi-Fusion and Marconi100.440 Views and opinions expressed are however those of441 the author(s) only and do not necessarily reflect those of the442 European Union or the European Commission. Neither the443 European Union nor the European Commission can be held444 responsible for them.445 REFERENCES446 [1] A. Khodak and R. Maingi, “Modeling of liquid lithium flow in porous447 plasma facing material,” Nucl. Mater. Energy, vol. 26, Mar. 2021,448 Art. no. 100935.449 [2] G. Federici, L. Boccaccini, F. Cismondi, M. Gasparotto, Y. Poitevin, and 450 I. Ricapito, “An overview of the EU breeding blanket design strategy 451 as an integral part of the DEMO design effort,” Fusion Eng. Design,452 vol. 141, pp. 30–42, Apr. 2019. 453 [3] P. Davidson, Introduction to Magnetohydrodynamics. Cambridge Texts 454 in Applied Mathematics, 2001. AQ:9455 [4] Guide Guide for the Verification and Validation of Computational Fluid 456 Dynamics Simulations, AIAA G-077–1998, AIAA, 2002. AQ:10 457 [5] S. Smolentsev et al., “An approach to verification and validation of MHD 458 codes for fusion applications,” Fusion Eng. Des., vol. 100, pp. 65–72, 459 Nov. 2015. 460 [6] S. Smolentsev et al., “Code-to-code comparison for a PbLi mixed461 convection MHD flow,” Fusion Sci. Technol., vol. 76, no. 5, 462 pp. 653–669, Jul. 2020. 463 [7] A. Khodak, “Numerical analysis of 2-D and 3-D MHD flows relevant 464 to fusion applications,” IEEE Trans. Plasma Sci., vol. 45, no. 9, 465 pp. 2561–2565, Sep. 2017. 466 [8] A. Patel, G. Pulugundla, S. Smolentsev, M. Abdou, and 467 R. Bhattacharyay, “Validation of numerical solvers for liquid metal 468 flow in a complex geometry in the presence of a strong magnetic field,” 469 Theor. Comput. Fluid Dyn., vol. 32, no. 2, pp. 165–178, Apr. 2018. 470 [9] N. L. Gajbhiye, P. Throvagunta, and V. Eswaran, “Validation 471 and verification of a robust 3-D MHD code,” Fusion Eng. 472 Des., vol. 128, pp. 7–22, Mar. 2018. [Online]. Available: 473 https://www.sciencedirect.com/science/article/pii/S0920379618300358 474 [10] S. Sahu and R. Bhattacharyay, “Validation of COMSOL code 475 for analyzing liquid metal magnetohydrodynamic flow,” Fusion 476 Eng. Des., vol. 127, pp. 151–159, Feb. 2018. [Online]. Available: 477 https://www.sciencedirect.com/science/article/pii/S0920379618300115 478 [11] A. Blishchik, M. van der Lans, and S. Kenjereš, “An extensive numerical 479 benchmark of the various magnetohydrodynamic flows,” Int. J. Heat 480 Fluid Flow, vol. 90, Aug. 2021, Art. no. 108800. 481 [12] C. Alberghi, L. Candido, R. Testoni, M. Utili, and M. Zucchetti, “Further 482 verification and validation of comsol magnetohydrodynamic models for 483 liquid metal breeding blankets technologies,” NOT PUBLISHED YET, 484 Tech. Rep., 2021. AQ:11485 [13] E. M. de les Valls, “Development of a simulation tool for MHD flows 486 under nuclear fusion conditions,” Ph.D. dissertation, Dept. Phys. Nuclear 487 Eng., Universitat Politècnica de Catalunya, Barcelona, Spain, 2011. 488 [14] D. Di Giulio, D. Suarez, L. Batet, E. M. de les Valls, and L. Savoldi, 489 “Analysis of flow channel insert deformations influence on the liquid 490 metal flow in DCLL blanket channels,” Fusion Eng. Des., vol. 157, 491 Aug. 2020, Art. no. 111639. 492 [15] D. Suarez, E. Iraola, C. Lampón, E. M. de les Valls, and L. Batet, “Liquid 493 metal MHD flow influence on heat transfer phenomena in fusion reactor 494 blankets,” Fusion Eng. Des., vol. 170, Sep. 2021, Art. no. 112503. 495 [16] M.-J. Ni, R. Munipalli, N. B. Morley, P. Huang, and M. A. Abdou, 496 “A current density conservative scheme for incompressible MHD flows 497 at a low magnetic Reynolds number. Part I: On a rectangular collo498 cated grid system,” J. Comput. Phys., vol. 227, no. 1, pp. 174–204, 499 Nov. 2007. 500 [17] R. Issa, “Solution of the implicitly discretised fluid flow equations by 501 operator-splitting,” J. Comput. Phys., vol. 61, pp. 39–65, Jan. 1985. 502 [18] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby, “A tensorial approach to 503 computational continuum mechanics using object-oriented techniques,” 504 Comput. Phys., vol. 12, no. 6, p. 620, 1998, doi: 10.1063/1.168744. 505 [19] B. Y. J. A. Shercliff, “Mathematical proceedings of the Cambridge 506 Cambridge philosophical society: Steady motion of conducting fluids 507 in pipes under transverse magnetic fields,” Math. Proc. Cambridge Phil. 508 Soc., vol. 49, pp. 136–144, Jan. 1953. 509 [20] J. C. R. Hunt, “Magnetohydrodynamic flow in rectangular ducts,” 510 J. Fluid Mech., vol. 21, pp. 577–590, Apr. 1965. 511 [21] C. B. Reed, B. F. Piclolglou, T. Q. Hua, and J. S. Walker, “Alex results— 512 A comparison of measurements from a round and a rectangular duct with 513 3-D code predictions,” in Proc. IEEE 12 Symp. Fusion Eng., Monterey, 514 CA, USA, Jun. 1987, pp. 1267–1270. 515 [22] X. Albets-Chico, E. V. Votyakov, H. Radhakrishnan, and S. Kassinos, 516 “Effects of the consistency of the fringing magnetic field on direct 517 numerical simulations of liquid–metal flow,” Fusion Eng. Des., vol. 86, 518 no. 1, pp. 5–14, Jan. 2011, doi: 10.1016/j.fusengdes.2010.07.014. 519 [23] I. Kirillov, “Present understanding of MHD and heat transfer phenom520 ena for liquid metal blankets,” Fusion Eng. Des., vol. 27, nos. 1–2, 521 pp. 553–569, Mar. 1995. 522 [24] G. Authié, T. Tagawa, and R. Moreau, “Buoyant flow in long vertical 523 enclosures in the presence of a strong horizontal magnetic field. Part 2. 524 Finite enclosures,” Eur. J. Mech. B, Fluids, vol. 22, no. 3, pp. 203–220, 525 May/Jun. 2003. 526