Magnetohydrodynamics in free surface liquid metal flow relevant to plasma-facing components

While flowing Liquid Metal (LM) Plasma-Facing Components (PFCs) represent a potentially transformative technology to enable long-pulse operation with high-power exhaust for fusion reactors, Magnetohydrodynamic (MHD) drag in the conducting LM will reduce the flow speed. Experiments have been completed in the linear open-channel LMX-U device [Hvasta et al 2018 Nucl. Fusion 58 01602] for validation of MHD drag calculations with either insulating or conducting walls, with codes similar to those used to design flowing LM PFCs for a Fusion Nuclear Science Facility [Kessel et al 2019 Fusion Sci. Technol . 75 886]. We observe that the average channel flow speed decreased with the use of conducting walls and the strength of the applied transverse magnetic field. The MHD drag from the retarding Lorentz force resulted in an increase of the LM depth in the channel that ‘piled up’ near the inlet, but not the outlet. As reproduced by OpenFOAM and ANSYS CFX calculations, the magnitude and characteristics of the pileup in the flow direction increased with the applied traverse magnetic field by up to

. We observe that the average channel flow speed decreased with the use of conducting walls and the strength of the applied transverse magnetic field. The MHD drag from the retarding Lorentz force resulted in an increase of the LM depth in the channel that 'piled up' near the inlet, but not the outlet. As reproduced by OpenFOAM and ANSYS CFX calculations, the magnitude and characteristics of the pileup in the flow direction increased with the applied traverse magnetic field by up to 120%, as compared to the case without an applied magnetic field, corresponding to an average velocity reduction of ∼45%. Particle tracking measurements confirmed a predicted shear in the flow speed, with the surface velocity increasing by 300%, despite the 45% drop in the average bulk speed. The MHD effect makes the bulk flow laminarized but keeps surface waves aligned along the magnetic field lines due to the anisotropy of MHD drag. The 3D fringe field and high surface velocity generate ripples around the outlet region. It was also confirmed that the MHD drag strongly depends on the conductivity of the channel walls, magnetic field, and volumetric flow rate, in agreement with the simulations and a developed analytical model. These validated models are now available to begin to determine the conditions under which the ideal LM channel design of a constant flow speed and fluid depth could be attained. * Authors 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.

Introduction
The use of flowing Liquid Metal (LM) Plasma-Facing Components (PFC) is a promising option for the protection of the solid divertors of magnetic fusion devices [1][2][3]. Fusion science experiments were carried out in tokamaks, such as CDX-U [4], T-11M [5], NSTX [6], FTU [7], HT-7 [8], and EAST [9], as well as other devices including IFMIF [10], TELS [11], and Magnum-PSI [12,13]. Most of these studies concluded that LM PFCs are beneficial for exhausting high heat flux; under the proper conditions, energy confinement is also improved. Using a large area of liquid lithium as the limiter, CDX-U observed a significant increase in the Ohmic plasma confinement [4]. Employing a Capillary Porous System (CPS) liquid lithium limiter as an evaporator, FTU extended the effective range of lower hybrid current drive towards the domain relevant for fusion reactors with a low recycling lithium wall [7]. Meanwhile, Morgan et al demonstrated that a CPS with liquid Sn or lithium can operate under an 18 MW m −2 heat load from a Magnum-PSI linear device without apparent damage [12,13]. A tray of static liquid lithium with a CPS structure was employed as the divertor plate in NSTX: lithium became saturated with implanted hydrogenic species, with impurities accumulating on the surface [6]. The experience reinforces the need for forced convection to be incorporated into the flowing LM divertor designs under development. A slow-Flowing Liquid Lithium limiter (FLiLi) with ∼100 µm film was tested on HT-7 and EAST [9,14]. The limiter was shown to be compatible with high-performance plasmas under different operating conditions in the EAST, and even improved plasma performance in many cases. Due to these advantages, several under-construction and future fusion devices are planning to test flowing LM PFCs. The Italian Divertor Test Tokamak (DTT) project is considering using tin and/or lithium as the LM material [15]. Recently, Khodak and Maingi used the customized ANSYS CFX solver to simulate the free-surface flow of LM on top of a solid divertor assembly and included the effect of heat transfer [16]. Smolentsev et al carried out simulations of a fast-flowing liquid lithium divertor to assess the viability of such PFC designs for the Fusion Nuclear Science Facility (FNSF), identifying regimes of (1) accelerating flow, (2) decelerating flow, (3) non-monotonic gradients with hydraulic jumps, and (4) the desirable constant flow speed and LM channel depth [17].
Several designs of LM-PFCs have been proposed in the past four decades. These designs broadly fall under two categories: fast-flowing and slow/static LM surfaces. Slow/static LM-PFCs, being self-healing, primarily target improving the plasma performance [18] and protecting the underlying substrate [19]. Heat is transferred from the liquid to the solid through thermal conduction and dissipated by an active cooling mechanism. A layer of vapor may form on the surface of the liquid, shielding it from the heat flux. In designs that experience vapor shielding, excessive evaporation remains a concern due to the limited capacity to control the liquid temperature [3]. Moreover, impurity compounds that form on the slowly-moving LM surface or CPS surface could discount the benefits of the improved heat transfer and LM free-surface refreshment [3,6,20]. Another implementation of LM-PFCs is to use a fast-flowing liquid film that traverses the length of the divertor target at high velocities [1,3,[21][22][23][24]. Yet another candidate is the use of vertically oriented LM jets [25]. The considerable Magnetohydrodynamic (MHD) drag that conducting liquids experience in magnetic fields is of significant interest [25][26][27][28][29]. While these effects have been evaluated for fusion blankets, the applications considered for divertor PFCs, and possible distortions to the flow channels, represent a complex design challenge for free surface flow concepts.
Past analytical studies of LM free-surface MHD flow commonly focused on the steady, fully developed film flow in a rectangular channel of arbitrary electrical conductivity with a constant applied coplanar or inclined magnetic field [24,25,[30][31][32]. In such flows, the governing equations are reduced to two spatial dimensions in the cross-sectional plane of the flow, greatly facilitating their solution. Perhaps the most significant theoretical result for LM film flows has been obtained by Shishko [33], who employed the Galerkin-Kantorovich method to derive an analytical solution using two essential functions in the Z-direction, the parabola, and the Hartmann profile. However, the analytical model has not yet been validated with experimental data.
Results of experiments on LM flow under the influence of magnetic fields are presented in past studies [26,31,32,[34][35][36]. The experimentally achieved Hartmann numbers (formally defined in section 2.1) <∼100 were relatively low, and the relatively high Reynolds number led to turbulent flow. Based on the assumed achievement of fully developed flow, the LM height along the flow (streamwise) direction calculated by 1D or 2D theoretical methods qualitatively matched experimental results [31,32,36]. It is unknown if the steady, fully developed film flow regime can be realized in a practical LM freesurface divertor application. Ying, Gao, Morley et al conducted extensive experimental and numerical studies on jet and film flows using a Gallium alloy in the Magneto-ThermO fluid Research (MTOR) facility, which utilized tokamak-like 3D magnetic fields. Their results indicated that the combination of spanwise and surface-normal magnetic fields caused various adverse effects, such as gradual or abrupt increases in the film thickness and detachment of the LM flow channel from the side walls [23,25,35]. They found that these 3D effects cannot be produced by a 2D model. Using the 3D simulation code HIMAG, however, they could only qualitatively reproduce the free surface LM behavior observed in the experiment [37]. To alleviate the power requirements of fast-flowing solutions, injecting an electric current through the LM for propulsion [16,22] and also the use of the Thermoelectric-MHD effect [38] were proposed. External J × B force was used to offset the negative effects from MHD drag, as well as oppose any J × B forces created by plasma currents. Finally, as an alternative to having a single flow path, the divertorlets concept introduces many flow paths by juxtaposing vertically oriented inlets and outlets [39,40].
Although efforts have been made to study the MHD effects on LM film flows from both engineering and laboratory perspectives, fully-validated predictive modeling capability is still lacking. Some issues remain to be solved, e.g. (i) validating 3D code with experimental data in the developing LM flow regime and (ii) characterizing the surface state in the MHDdominated free-surface flow. Additional solutions to mitigate the MHD drag are needed, including mitigation of drag from conducting walls.
To address this gap in the LM PFC knowledge, we conducted an experimental study of MHD effects in a free-surface LM flow with different types of conductive and insulating walls in the external magnetic field generally parallel to the free surface and perpendicular to the flow direction. Two 3D Computational Fluid Dynamics (CFDs) codes were compared with experimental data, with one of the codes also having been used to simulate the FNSF [41]. The remainder of this paper is organized as follows. Section 2 explains the LM experimental system and simulation setup. Experimental data on the MHD effect on the free-surface LM flow and numerical code validation are presented in section 3. Section 4 discusses free-surface flow physics, surface state, contact angle, 3D magnetic filed effect, and proposes a potential solution to alleviate the MHD drag. The last section summarizes our findings and forecasts the future work.

Experiment setup
The experimental results presented in this paper were obtained using the LM eXperiment-Upgrade (LMX-U) test facility [42][43][44], which was designed and built to investigate freesurface, liquid-metal flows, and MHD effects relevant to LM-PFC development. A depiction of Liquid Metal eXperiment Upgrade (LMX-U) can be found in figure 1. LMX-U uses a GaInSn eutectic alloy, Galinstan, as the working LM because it is a liquid at room temperature and has low reactivity and non-toxic properties (see table 1). A gear pump is used to circulate Galinstan throughout the closed-loop system continuously. A flow meter monitors the flow rate and shows that a constant flow rate is maintained for a given rotary pump speed despite variations in LM hydrostatic pressure. Argon is continuously injected into the channel to minimize LM oxidation, and its pressure is kept slightly above atmospheric.
A sliding laser sheet and camera measured LM depth at different points along the channel. The laser-sheet depth measurements were calibrated by comparing the video data to corresponding height measurements taken using solid blocks with known thicknesses [43]. Measurements indicate that the variation of flow height in the line parallel with the magnetic field is negligible. Therefore, LM height in the spanwise direction is an averaged value. The averaged velocity can be calculated by u = Q/ (zw), where 'Q' is volume flow rate, 'z' is the liquid height, and 'w' is the inner chute width. To measure the surface velocity of the LM, small spherical insulating particles with a diameter of ∼5 mm were introduced into the flow from the inlet. Particle is made of plastic or wood, with a density much less than the Galinstan (around 500-100 kg m −3 ). These particles floated on the surface and were carried along with the ambient LM flow. An overhead camera recording data at 120 frames per second was used to record particle motion and calculate the particle velocity, which is expected to closely approach the top surface velocity. During LMX-U operation, Galinstan was pumped through a rectangular duct with a width = 109 mm and a length of 1200 mm. The chute was held in the horizontal position for this work, parallel to the floor. The substrate and side walls (Hartmann walls) are made from insulating acrylic. To investigate the impact of electrical conductivity of the boundary conditions on the flow, a single piece of metal sheet was bent into a '⊔' shaped liner, which fit into the acrylic channel tightly. Four liners were embedded into the insulating base to test different materials. Table 2 lists the physical properties of the different liners used in experiments. The wall conductance ratio is defined as C W = 2t w σ w / (wσ LM ), where t w , σ w , σ LM , and 'w' are side wall thickness, wall conductivity, LM conductivity, and the width of the channel within the liner, respectively. It is noted that, except where specifically noted, nozzles were not used to control the initial inlet LM height, to avoid the onset of a hydraulic jump [44].
The LMX duct was installed horizontally within the air gap of a C-shaped, water-cooled electromagnet with a length of 74 cm in the flow direction. The magnet provided a magnetic field up to 0.33 T with nearly uniform (∼4% field strength variation) across the width of the test section. Along the streamwise direction, the electromagnet's output, measured using a LakeShore model 410 gauss meter is shown in figure 2. A smooth fit was used to calculate the magnetic field gradient as input to the MHD drag simulations. Magnetic field density measurements were taken at the middle of the channel bottom where it contacts the floor (Y = 0, Z = 0); see figure 1 for the coordinate system). X = 0 was set as a reference point at the edge of the electromagnet near the inlet. The magnetic field varied strongly near both ends of the electromagnet: two symmetric fringe regions with positive and negative gradients were observed, comparable to the toroidal  magnetic field gradients in a midsize tokamak [23].
Quantifying key dimensionless quantities is essential to determine the regime of a particular flow. The Reynolds number, defined by Re = uR h /ν, quantifies the ratio of inertial forces to viscous forces, where R h = zw/(w + 2z) is a characteristic length. The Hartmann number, Ha = 0.5Bw √ σ LM / (νρ), relates the magnitude of electromagnetic forces to the viscous forces. The magnetic Reynolds number, R em = uwσµ 0 , which is the ratio of magnetic advection to magnetic diffusion. The Stuart number, N = Ha 2 /Re, quantifies the magnitude of electromagnetic forces to inertial forces.

Governing equations
The Galinstan flow is described as an incompressible, laminar flow of electrically conductive fluid with homogeneous physical properties. The flow in a magnetic field is governed by the MHD equations, namely the conservation equations for momentum (equation (1)), mass (equation (2)), and current (equation (3)). The MHD drag is reflected by the Lorentz force term ⃗ J × ⃗ B in the Navier-Stokes equation, which also includes surface tension force, given by ∇z represents the curvature of the interface. Ohm's law (equation (4)) resolves the electric current, and the electric potential is determined by the Poisson equation realized by introducing (equation (4)) into (equation (3)): where ⃗ v is the local flow velocity, p is the hydrodynamic pressure, µ is the dynamic viscosity, ⃗ J is the current density, ⃗ B is the applied magnetic field, σ is the electrical conductivity, and φ is the electric potential.

Computational model
Analytic solutions to these equations are available for simple problems. However, one can solve these equations numerically by implementing the MHD equations into an existing CFD code, which can serve as a design tool after validation. The analysis presented in this paper allows validation of a customized version of the commercial multi-purpose CFD code ANSYS CFX [41,46] and OpenFOAM, an open-source CFD toolbox with a built-in electromagnetic module [47]. For both codes, free-surface MHD flow is modeled in the framework of a fully 3D Volume of Fluid (VOF) method with the free surface shape determined by the solution. Nonuniform, structured, hexahedral meshes are used in all simulations shown in this paper. In VOF methods, the auxiliary variable, α, is the volume fraction. It ranges between zero and one and represents the ratio of primary to secondary fluid volume in a discrete computational cell. Therefore, the interface between two immiscible fluids is represented by a sharp jump in α. If the volume fraction is advected accurately and the interface thickness is kept constant in space and time, VOF schemes have been shown to demonstrate a high level of accuracy and mass conservation. Electric current and potential distributions are properly computed, accounting for the different electrical conductivities of the two materials. For OpenFOAM, the Multi-Dimensional Limiter for Explicit Solution method is an algebraic VOF method that modifies the advection of the volume fraction by adding an interface compression velocity term [48]. This greatly reduces the smearing of the interface and attempts to control its thickness throughout the duration of the simulation. The partitioned approach is used to achieve coupling between two phases. Using this approach, the equations are solved separately in their respective domains, and information is exchanged across domains iteratively by imposing appropriate boundary conditions [49]. The surface tension force acting on the interface is modeled in this work using the continuum surface force method [50]. In the simulation analysis, the surface item is included in the calculation. In practice, the effect of the surface tension is neglectable in the LMX channel flow because (1) the channel is relatively large compared to the capillary length; (2) the surface tension at the interface is reduced when LM is moving; and (3) the dynamic contact angle, θ, is set as 90 • with respect to gravity direction for all walls without specific notation, for the sake of the simplicity. However, the contact lines, which are boundaries between solid walls and the LM, must be specified. We observed convex meniscuses for acrylic and stainless-steel walls, but not for copper and brass walls, due to a better wettability between Galinstan and copper than acrylic and stainless steel. A slip boundary condition at the top of the channel was used, implying zero viscous stress and no penetration, while a no-slip boundary condition was used for all other inner walls of the channel. The continuity of the normal component of electric current between solid and liquid was employed to model conductive walls, while the zero electric flux boundary condition was used for insulated walls. Furthermore, the symmetry of the geometry in the Y-direction was assumed allowing the simulation of a half channel with a symmetry boundary condition. It is noted that the OpenFOAM simulation is a predictive analysis; with magnetic field, channel geometry and material, and flow rate serving as inputs. The ANSYS CFX setup used in the present analysis assumes a flow domain range from −350 mm to 882 mm (figure 1) without an assumption of symmetry and with the inlet and outlet pressure and a flow rate defined by the experimental data.
To adequately resolve the thin Hartmann and side layers, meshes were refined near the walls. From a numerical point of view, we have ensured that the grid resolution used in the simulations is sufficient, and the obtained results remain nearly constant with increased resolution. E.g. OpenFOAM, two meshes composed of hexahedra are used, with the parameters listed in table 3. The results with two meshes are discussed in section 4.
Due to the relatively small magnetic Reynolds number (∼0.1) in the LMX condition, the induced magnetic field is negligible relative to the applied field. This allows us to reduce the MHD equations to a single current conservation equation, equation (4). The main advantage of using this approach for evaluating the Lorentz force is numerical simplicity and efficiency since only a single scalar Poisson equation needs to be solved. The MHD forces exerted on the flow are computed using the non-induction approximation that utilizes the electric potential as a solution variable in OpenFOAM. This assumption is examined by comparing it to results obtained from the magnetic vector potential solver. Both the electrical potential and magnetic vector potential solvers are able to produce results that agree well with the analytical solution with comparable accuracy in the Schercliff flow and Hunt flow [49]. A code, titled 'FreeMHD', was used to solve the electric potential Poisson equation and calculate the Lorentz force. The FreeMHD code has been developed as a valuable tool for the analysis of free-surface MHD flows, using the OpenFOAM framework as its foundation. FreeMHD can be found on Github: https://github.com/PlasmaControl/FreeMHD 4 .

Analytical model
To understand the LM accumulation in the chute effectively in terms of Lorentz force, which acts on the LM flow, a theoretical analysis is conducted. Conservations of mass and momentum describe steady, free-surface flows. For a control volume with upstream and downstream cross-sections, a schematic diagram of this analysis is shown figure 3.
The mass conservation is given by: where 'z' is the LM surface height; 'du' and 'dz' are LM flow velocity and height variation in this control volume, respectively.
Assuming the flow along the Y direction (parallel with the B direction) is negligible, and the flow is steady and incompressible in a prismatic rectangular channel, the momentum conservation of a control volume can be written as: x wzdx is the Lorentz force in the control volume (the LM flowing in the perpendicular magnetic field experiences a net electromagnetic body force); τ is average shear stress acting on the channel bottom and sides. β = where 'A' is area element for integration and 'v' is the local velocity, is a momentum coefficient. To simplify, β is unity if the velocity is assumed to be uniform.
For obtaining a solution, a fully developed channel flow is assumed. For the open channel flow, based on the Darcy-Weisbach equation [51], the momentum balance is expressed as: where λ is the Darcy friction factor. Combining equations (6) and (7), we obtain ( For the flow with the pressure gradient induced by the MHD drag, equation (8) can be converted as: A complete expression for the MHD flow is complex. The MHD drag is approximated as F MHD ≈ C M σuB 2 , including B 2 , u or Q, and C M , which is the dimensionless MHD drag force coefficient. Note that the coefficient C M is intrinsic to the specific wall conductivity and geometry. As stated by Molokov and Reed [29], a simplified and insightful theoretical result for LM flows has been obtained by Shishko [33], who provided a solution for the fully developed flow on the inclined plate. The LMX flow is developing since on the horizontal plate fully developed flow is not possible. This could cause a mismatch between analytical calculations and experiments. He used a Galerkin method with two basic functions in the Y-direction, the parabola, and the Hartmann profile, to approximate the velocity profile along the magnetic field lines. Results of Shishko model were confirmed by the fully developed CFD analysis using CFX [52]. In the analytical model we used a drag formula from Shishko paper and applied it at each cross-section for the developing flow on the horizontal plate as shown below: where, E = 4Ha (1+CwHa) (1+Cw) , D = − 1 . Substituting equation (10) into equation (9), we obtain: Starting from the downstream end of the channel, the MHD retarding force in each cross-section is integrated toward the inlet, causing the LM to accumulate, i.e. pileup, and increase in height upstream. With the integration of equation (11), the relationship between LM height and X location can be expressed as equation (12) x = x out − ρwQ µE where the x out and z out are the X-location and LM height of the outlet, respectively. Identifying that the LM flow is over the horizontal datum, and the surface height is larger than critical depth, z c = 3 √ Q 2 / (gw 2 ), the LM flow approaches free overfall without a hydraulic jump [44]. The theory calculation starts from the electromagnet edge near the outlet, x out = 814 mm and z out = 8.99 mm. At the location of X = 814 mm, the magnetic field in the fringe decays to half of unfirm value. The z out = 8.99 mm is the calculated LM height at X = 814 mm, assuming the flow from the magnetic field fringe to the overfall edge is Darcy-Weisbach laminar flow, which is not impacted by the magnetic field. For the laminar flow, the Darcy friction factor is expressed as the equation (13 Combining equation (13) with equation (8), we obtain Starting from the overfall edge, the X-location (x L ) and corresponding flow height (z L ) can be calculated by solving the equation (14): where x c = 882 mm and z c = 8.85 mm. The liquid height at the overfall edge is set as the critical depth, z c .

MHD drag with a conducting wall
Significant MHD drag on the Galinstan flow was observed when inserted copper liners were used as walls. Figure 4 displays the evolution of LM height and averaged velocity near the inlet (X = −250 mm) as a 0.3 T magnetic field was turned on and off. Before the electromagnet was turned on, the average velocity was about 0.16 m s −1 , and the LM height was ∼16 mm, with Re ∼ 6700. After the electromagnet was turned on, the LM height grew to ∼30 mm in less than 5 s and reached a stable elevation of ∼ 35 mm in the next 15-20 s. The average velocity was reduced proportionately by ∼45%. Once the electromagnet was turned off, the LM height dropped quickly (in <5 s) to the previous height, and the depth returned to the initial value without a magnetic field. The pump drives the LM from the bottom up into the chute through the rectangular opening, similar to an underground fountain. Note that Galinstan was pumped into the copper chute at X = −350 mm vertically with constant pressure. Due to the vertical velocity component in the flow near the inlet, ripples formed on the free surface but were quickly damped in the streamwise direction. At X = −250 mm, the location of the data measurement, the LM height approached a relatively steady state with fluctuations <0.5 mm. The phenomenon of LM height increase is analogous to the pressure drop due to MHD drag in an internal duct or pipe flow. Figure 5 shows the LM surface height profiles along the flow direction before and during the ∼0.3 T magnetic field applied after the liquid height reached a steady state. The blue curves (no magnetic field) show a modestly decreasing LM height from ∼16 mm to ∼10 mm from upstream to downstream. In comparison, the orange curves (B = ∼0.3 T) show a strongly decreasing LM height from 35 mm to 8 mm along the channel length. Without an applied magnetic field, viscous drag reduces the flow height from the inlet to the outlet; this phenomenon is also termed 'pileup'. With the 0.3 T magnetic field (Ha = ∼635, N = ∼60) applied, the LM height near the inlet rose by 120%, while the outlet height decreased slightly. The LM volume in the channel increased by ∼65% from B = 0 to B = 0.3 T (more LM volume stays in the channels). While observing the same LM height in the spanwise direction, it is believed that the ratio of the integrated area under the curve of the LM height is the ratio of the volume. To calculate the area, trapezoids were used to approximate the curve, and the sum of the areas of these trapezoids was calculated. Note that this experimental flow is developing because the chute length of LMX-U is generally insufficient to achieve a fully developed flow. Similarly, a future tokamak reactor cannot obtain a fully developed flow since divertor length is insufficient. It is noted that although two magnetic fringes with large magnetic gradients (figure 2) exist at the upstream and downstream ends of the flow, the slope of the surface height, ∇h = dz/dx , became moderate in these two regions. Thus, it is evident that the magnetic field intensity plays a stronger role in the LM pileup than the magnetic field gradients that occur near the electromagnet. The height of the LM surface was measured along the entire chute length via a laser profilometer, except where three metal fasteners blocked the view (indicated by the three data gaps in figure 5). The liquid height was effectively an average value across the whole width because measurements showed little LM height variation in the transverse direction.

Analytical model validation
To validate the analytical model, comparisons between experiments and calculations under the same conditions are conducted. The solid red line of figure 5 shows that LM surface height, calculated by the theory, varies gradually with the X location. While the calculation follows the trend of the observed LM height, the computed slope is ∼30% less than the data. This mismatch probably comes from several assumptions in section 2.4.
For closed channel flow, the flow state of the LM is significantly affected by the magnetic field, flow rate, and wall electrical conductivity [28]. To investigate the effect of these parameters on MHD drag in open channel flow, we performed numerous calculations and experimental investigations. Figure 6 shows the LM height measured at different X locations in steady state changes with Ha, scanned from ∼160 to ∼670. The LM height near the beginning of the electromagnet (X = −240 mm, 0 mm, and 150 mm) increased linearly as a function of the Ha. The computed LM surface height at X = 150 mm evolution as Ha is superposed in figure 6. The theory confirms the linear relationship, with a computed slope within ∼10% of the measured slope. The MHD drag effect can be observed for highly conductive copper walls at a small magnetic field, 0.08 T, Ha ∼ 160.
Similarly, both experimental and theoretical results indicate that LM accumulation is linearly related to the flow rate ( figure 7). The linear fit slopes of theoretical calculations at X = 0 mm and X = 400 mm are in agreement with the experimental data within an approximately 10% deviation. The theoretical LM height does not entirely correspond with  the experimental data; however, by comparing X = 400 mm with X = 0 mm, the difference between the two gradually diminishes towards the outlet, also as illustrated in figure 5. To account for the influence of flow rate on the LM height measurement, the LM height at X = 720 mm is subtracted. The theoretical slopes, generated by scanning the magnetic field and flow rate, exhibit excellent consistency with experimental data. Consequently, these calculations prove useful in comprehending the evolution of LM height and forecasting the trend of the LM surface height.

CFD simulation validation
A second purpose of the current experiments is to validate the numerical codes. To our best knowledge, no simulation code has been fully validated against experimental results for free-surface MHD channel flows. In this study, ANSYS CFX and OpenFOAM were used to simulate the experimental results. For both codes, the simulation setup is the same as the experiment parameters. Laminar flow was set for simulation cases with a 0.3 T magnetic field. In practice, flow laminarization should occur because Ha/Re exceeds the critical Ha/Re = 0.0075 in a transverse magnetic field [28]. For future reactors, e.g. FNSF, a strong toroidal magnetic field, ∼4-6 T, leads to a laminar flow even with a high velocity of ∼10 m s −1 [17]. A comparison between numerical and experimental results is shown in figure 8. Both ANSYS CFX and OpenFOAM agree with LM height profiles measured in experiments. OpenFOAM results also match experimental LM surface height near the inlet and outlet, with discrepancies <5% for both magnetic fields. ANSYS CFX properly interpolates LM surface height shape along the streamwise direction between the inlet and outlet. In addition, for the fluid height profile without a magnetic field, ANSYS CFX obtains an agreement with experimental data by including a turbulence model, shown as the comparison between the black line and the black dot in figure 8.
Since the LM-PFC in fusion applications is exposed to plasma heat and particle fluxes, the free-surface velocity is crucial for deciding the heat transfer and recycling characteristics. Figure 9 displays surface velocity distribution along the streamwise direction measured by the particle tracking method     and calculated by OpenFOAM and ANSYS CFX. Both simulation codes produce a surface velocity trend consistent with experimental results, with surface velocity ramping up gradually along the streamwise direction. The surface velocity is increased by 400% above the average velocity from two simulations and experimental data. There is a slight discrepancy in the velocity value between simulations and experiments, possibly because air friction reduces the velocity of the introduced particles below the fluid surface velocity. Figure 10 compares experimental surface velocities and simulation velocity profiles in X-Z (Y = 0) plane between 0.2 T and 0.3 T magnetic fields. The increment of surface velocity with increasing B is consistent between the experimental data and the simulations, indicating that the MHD effect with increasing magnetic field accelerates the surface flow. Figures 10(b) and (c) show that the liquid enters the uniform magnetic field region, and MHD drag diverts most (∼60%-90%) of the flow toward the free surface. For the higher magnetic field case, the velocity in bulk tends to be constant, and a 'slow regime' occurs, characterized by velocity distribution where most of the flow is carried by one jet, which forms near the free surface with virtually no movement of the LM in bulk. For small Ha numbers, the bulk flow becomes faster and is characterized by monotonic velocity profiles along the flow depth. An increase in Ha corresponds with the progressive definition of a surface region at the LM flow, clearly shown in figure 10(c).

Surface state and 3D characters
In most practical applications, 3D effects are seldom avoidable, and the applied magnetic field often varies in magnitude and direction in the domain. In LMX, although the flow is impacted by the two fringes of the magnetic field, no side wall detachment or inclined surface (X-Y plane) along the magnetic field direction (Y direction) was observed due to lacking an appreciable magnetic field of the Y component. Figure 11 shows the surface height in the spanwise direction is pretty constant in the three X locations, which is simulated by the 3D codes. Figure 11(c) displays the velocity profile and induced electric current streamlines in the Y-Z cross section with X = 500 mm for the 0.3 T case. Since electric currents cannot cross the top surface, all induced currents flow through the top side layer. Since the electric conductivity of the copper wall is near ∼20 times that of Galinstan, the induced currents preferentially flow through the lateral walls and close through the bottom wall. The Lorentz force is minimized in the top side layer near the free surface because all currents flow parallel to the magnetic field. In contrast, the bulk flow region, where induced currents interact with the magnetic field, contains most of the Lorentz force as the retardant force. Therefore, at a constant flow rate, a considerable part of the flow is through the side layer near the free surface, which may compensate for the reduction of fluid flow through the bulk flow region. The absence of the no-slip condition on the free surface also contributes to the enhanced surface velocity. For a closed duct flow, the generated MHD effect modifies velocity profiles, and two no-slip conditions exit, forming an 'M' shape flow [28]. In open-channel flow, the bottom floor has a nonslip boundary, but the top surface is free to slip. By analogy with the Hunt flow, the high-velocity layer thickness is proportional to Ha −0.5 [28]. Figure 10(c) depicts an almost constant thickness of the high-velocity layer along the streamwise direction in the uniform magnetic field region. In the X-direction, the downwards-sloping free surface accelerates the liquid further. The ANSYS CFX calculation also obtains similar current streamline and velocity distribution. As the LM flows into the outlet region, inertial force keeps the LM flow at a high velocity, resulting in a lower LM height than the LM flow without a magnetic field. In the region near the inlet where the magnetic field is weak and averaged velocity is low, mixing due to vortices generated from the inlet renders the velocity profile relatively uniform in the Z-direction ( figure 11(a)). As the liquid enters the uniform magnetic field region, MHD drag in the core region retards flow and the fast surface region appears ( figure 11(b)). Suppression of turbulence due to MHD drag is also clearly seen in figures 11(b) and (c), as the flow becomes laminar as it enters the uniform magnetic field region.
However, surface fluctuations are changed by the co-planar and the 3D fringe field with a strong ∇B, shown in figure 12, which compares the LM surface state captured by the overhead camera with and without the magnetic field. 3D fluctuation due to the high Re in a turbulent state clearly appears in figure 12(a) without the magnetic field. Appling 0.3 T magnetic field led the surface more reflective and waves aligned along the magnetic field emerged as the red arrows indicated in figure 12(b). This is consistent with theory prediction [53] and other experiments [26]. Due to the anisotropy of MHD drag, which acts to hinder flow across magnetic field lines, vorticity perpendicular to the magnetic field is greatly subdued, and turbulent eddies are elongated along magnetic field lines. Figure 12(b) also shows ripples were generated in the downstream region. As approaches the outlet, the surface velocity is gradually increased, which could overcome surface restoring force, such as hydrostatic pressure and MHD drag. OpenFOAM simulations show, due to the presence of a strong ∇B, electric current lines near the free surface curl and flow in the same direction as the LM flow, creating a Lorentz force that points upwards. These two effects possibly combine to create small surface ripples that were also found in OpenFOAM simulation with the fine mesh, shown as the oscillation of the surface velocity in the region of X > ∼600 mm (orange line in figure 9). The non-uniformity of the magnetic field induces 3D electric currents, which not only flow within the plane of the applied magnetic field but also have a component that flows in the direction of the overall flow. This in turn causes an additional 3D MHD drag at the fringing region.
For a realistic fusion reactor design, including the variation of the toroidal and poloidal magnetic fields, fully 3D simulations of the flow in an LM divertor are carried out [49]. A flat substrate does not conform to the toroidal field lines, resulting in a significant surface-normal component that grows towards the chute's extremities in the toroidal direction as the field lines intersect the film at an increasing angle. In combination with the poloidal field component, this misalignment of the toroidal field creates an asymmetric surface-normal field that strongly impedes the flow on one side of the chute and accelerates it on the other. The flow velocity contains a degree of asymmetry due to the transverse gradient of the surface-normal magnetic field. Furthermore, stream-wise electric currents that arise due to the magnetic field gradient in the direction of the flow cause the liquid to separate from the dividing wall on the high surface-normal field side. These results are not consistent with previous studies that utilize reduced-order models that do not account for the influence of magnetic field variation [17]. This indicates that fully 3D simulations that do not omit any field or geometry variations are important to produce realistic predictions of the flow behavior under fusion machine conditions.

Influence of the wall conductance
The electric conductance of the channel walls influences the current distribution in the fluid and determines the flow pattern. For a closed rectangular flow, configurations of wall conductivities have been studied extensively, including Shercliff flow, Hunt flow, and Ufland flow [28]. The currents close through the relatively thin Hartmann and side layers for insulating walls. Because these layers are very thin, their electric resistance is high, so the current magnitude is small. For highly conducting walls, a significant fraction of currents may close through the walls in addition to that in the viscous layers. This increases the total magnitude of currents compared to insulating conditions. Consequently, we expect stronger Lorentz forces and a higher MHD drag with increasing wall conductance. Figure 13 shows the LM height difference between the inlet and outlet as a function of the wall conductance ratio, as table 2 shows. As expected, the LM height difference increased as the wall conductance ratio increased both in experimental data and simulations. Good agreement can also be observed between experimental data and numerical calculations for conducting walls and insulating walls. Note that the LM thickness difference roughly reflected the integrated MHD drag on the LM flow. As the wall conductance ratio increased from 0 (insulating wall) to 0.13 (a moderately conductive wall), the LM height difference increased by ∼300%, while the LM height difference increased only by ∼30% as the conductance ratio increased from 0.13 to 0.75. This result indicates that the MHD drag is more sensitive in the region with a small conductance ratio. To minimize the MHD drag in a high magnetic field environment, the C W should fall into the region with a small ratio by choosing a material with low conductance, e.g. stainless steel or tungsten. From the definition of wall conductance ratio, C W = 2t w σ w / (wσ LM ), is linearly relevant to the wall thickness (t w ). Reducing the wall thickness is a considerable method. This finding suggests that wall thickness can become a key factor in the magnitude of the MHD drag.

Discussion
One of the most distinguishing aspects of this work in comparison to other MHD duct flow experiments is the presence of a free surface. For the free surface flow without the confinement from the top wall, the pressure drop contributed by the retarding force is reflected in the varied surface height. From the downstream to the upstream side, the integrated retarding force on every Y-Z cross-section leads the LM to accumulate in the channel, and so the LM height increases gradually from the outlet to the inlet. The LM height profile is similar to the closed flow with high pressure in the inlet and low pressure in the outlet, but the free-surface flow leads to more LM inventory in the chute (∼65% more for the LMX channel). The 1D analytical model and simulation confirm this kind of accumulation and indicate the flow rate, magnetic field, and gravity are the key factor to the 'pileup'. In a fusion machine, where the divertor plate is typically inclined or vertical to the floor, gravity could be an important factor in driving the flow. As a result, free surface flow in the inclined plate is driven by gravity and inertia. Viscous and Lorentz forces oppose the flow. Depending on the force balance, the surface height can be (1) uniform, (2) accelerating, (3) decelerating, and (4) nonmonotonic, e.g. as in a hydraulic jump. Ideally, the flow thickness in the divertor should be as uniform as possible so that the plasma interacts with a 'flat surface' with minimal leading edges to concentrate the heat flux.
The surface velocity of the free surface is greatly enhanced by the strong MHD effect. It could be argued that a high freesurface film velocity is beneficial for fusion applications, both in divertors and flowing first walls since it allows for the fast refreshment of the liquid armor. In addition, Smolenstev et al reported that this kind of MHD flow for FNSF allows for about 20 K lower temperature at the strike point as compared to the slug flow [17]. Hirooka et al reported that hydrogen and helium recycling is noticeably reduced from a dynamic liquid Li surface [54]. At the interface between the liquid and the atmosphere, waves can be created by displacing fluid elements near the surface. The inertia of the LM with high surface velocity tends to generate surface instability. The anisotropy of MHD drag, which acts to hinder flow across magnetic field lines, vorticity perpendicular to the magnetic field is greatly subdued, and turbulent eddies are elongated along magnetic field lines. Gravity and surface tension act to restore surface equilibrium. This oscillation was observed in LMX flow and OpenFOAM simulation with a highly refined mesh (figures 14 (bottom) and 9). Once inertial force begins to dominate over surface tension, droplet formation becomes a very significant concern. Since the LM has a relatively high surface tension coefficient, the surface tension acts to restore equilibrium when perturbations result in significant stretching of the fluid surface.
Taking the open-channel MHD flow with stainless steel liners as an example, we present two velocity profiles in the X−Z plane calculated with two meshes, whose parameters are listed in table 3. The relatively coarse mesh is capable of resolving the details of the flow while reducing the computational cost. Using the fine mesh, more details of the vortices forming as the liquid exits the uniform magnetic field region are resolved. However, the difference of the liquid surface heights along the streamwise is negligible, as comparing the top and bottom panels of figure 14. This comparison verifies the sufficiency of the mesh resolution in the current study.
As mentioned in section 2.2, the contact angle is set as 90 • in the simulation. In experiments, we observed convex meniscuses for acrylic and stainless-steel walls, but not for copper and brass walls, due to a better wettability between Galinstan and copper than acrylic and stainless steel. To investigate the surface tension effect on the contact lines, a direct comparison of θ = 90 • and 160 • , calculated by OpenFOAM,  is shown in figure 15. Changing the contact angle leads to a slight difference in LM heights, with a maximum value of ∼1.5 mm. With a magnetic field, the LM height is observed to be decreased in the condition of the convex meniscus (180 • >θ>90 • ). This is the inverse of what happens in non-MHD flows.
The use of electrically insulating materials in contact with the LM significantly reduces the MHD drag compared to conducting walls by forcing the currents generated in the viscous layers. Insulating material coating (e.g. Al 2 O 3 [55], Y 2 O 3 [56]) on structure material is a possible alternative method to mitigate the MHD drag. However, the compatibility of LM with various insulating materials for long pulse lengths is challenging and requires significantly more development. It was shown that the MHD drag is sensitive to the conductance ratio, suggesting that reducing the thickness of walls is one way to mitigate the MHD drag. Multiple layers of material, e.g. a thin RAFM steel contacting with LM, an insulating layer, and a RAFM steel base, could be compressed together for a new structural material design. The thin RAFM layer is expected to reduce the total induced current and induced MHD drag significantly.
It was shown that the analytical model can qualitatively represent the trend of the LM free surface flow observed in the experiment and thus can be used in the initial design analysis of the free surface LM PFC. Detailed design requires CFD simulations, which quantitatively reproduce measured flows within the LMX-U condition. In a tokamak, where the divertor plate is typically inclined or vertical to the floor, gravity could be utilized to alleviate pileup issues. Nevertheless, a surfacenormal magnetic field also has a detrimental effect on the flow. The optimized curvature of the chute geometry in both toroidal and poloidal directions is a tentative design approach to minimize or eliminate the normalized component magnetic field throughout the film path [49]. In addition, partitioning the flowing liquid torus into multiple chutes using vertical walls may be helpful in limiting the induced toroidal currents [57]. However, further investigation of chute geometry, wall conductance, and other factors is necessary to design an optimal flow path that can achieve a stable flow and effectively meet the operational requirements.

Conclusion and outlook
To mimic the toroidal field in tokamaks, the MHD drag effect on an open-channel LM flow has been investigated with insulating and conducting walls in LMX-U. MHD drag causes the LM to ramp up from the outlet to the inlet, deforming the LM surface shape. At the region near the channel inlet, the LM depth significantly increased by up to ∼120% more than the case without a magnetic field, corresponding to an average velocity reduction of ∼45%. With an applied transverse magnetic field, the simulation demonstrates that the LM flow velocity profile was modified. With a magnetic field, the average velocity was reduced, but surface velocity was increased up to 300%. Simulations with OpenFOAM and ANSYS CFX reproduced experimental observations, showing that (1) induced currents interact with the magnetic field, causing a retarding Lorentz force; (2) near the free surface, the induced current flows parallel to the magnetic field and MHD drag is eliminated, leading to a high surface velocity. The MHD effect makes the bulk flow laminarized but keeps surface waves aligned along the magnetic field lines due to the anisotropy of MHD drag. The surface fluctuations are mitigated in the uniform magnetic field region except for waves parallel with the magnetic field line. The 3D fringe field and high surface velocity generate ripples around the outlet region. No side wall detachment or inclined surface along the magnetic field direction was observed due to lacking a normal magnetic field component. It was also observed that the MHD drag strongly depends on the conductivity of channel walls, magnetic field, and flow rate. The developed 1D analytical model semi-quantitively reproduces the experimental results. The use of a thin wall with good compatibility with the LM or electrically insulating materials significantly reduces the MHD drag compared to conducting walls by forcing the currents generated in the viscous layers.
While this is a good first step, we note that the LMX-U results differ in dimensionless parameters as compared to an FNSF (Re ∼ 1 × 10 5 , u in ∼ 10 m s −1 , Ha ∼ 1.5E5). Spanning the gap from LMX-U to FNSF can be done by further upgrades to LMX with a higher B ∼ 2 T (Ha ∼ 5E3) and higher initial velocity ∼3 m s −1 (Re ∼ 1E5). Also, experiments and simulations need to be carried out in a complex geometry magnetic field, e.g. MTOR-like toroidal magnetic field with a surface normal component. Once fully validated, numerical simulations such as by ANSYS CFX and OpenFOAM will provide LM-PFC designers with a useful tool to predictively compute free-surface liquid-metal behavior for future fusion devices. Nonetheless, the present tools can be used to determine the conditions under which a desirable flow speed and flow depth could be realized.