An extensive numerical benchmark of the various magnetohydrodynamic flows

There is a continuous need for an updated series of numerical benchmarks dealing with various aspects of the magnetohydrodynamics (MHD) phenomena (i.e. interactions of the flow of an electrically conducting fluid and an externally imposed magnetic field). The focus of the present study is numerical magnetohydrodynamics (MHD) where we have performed an extensive series of simulations for generic configurations, including: (i) a laminar conjugate MHD flow in a duct with varied electrical conductivity of the walls, (ii) a back-step flow, (iii) a multiphase cavity flow, (iv) a rising bubble in liquid metal and (v) a turbulent conjugate MHD flow in a duct with varied electrical conductivity of surrounding walls. All considered benchmark situations are for the one-way coupled MHD approach, where the induced magnetic field is negligible. The governing equations describing the one-way coupled MHD phenomena are numerically implemented in the open-source code OpenFOAM. The novel elements of the numerical algorithm include fully-conservative forms of the discretized Lorentz force in the momentum equation and divergence-free current density, the conjugate MHD (coupling of the wall/fluid do- mains), the multi-phase MHD, and, finally, the MHD turbulence. The multi-phase phenomena are simulated with the Volume of Fluid (VOF) approach, whereas the MHD turbulence is simulated with the dynamic Large-Eddy Simulation (LES) method. For all considered benchmark cases, a very good agreement is obtained with available analytical solutions and other numerical results in the literature. The presented extensive numerical benchmarks are expected to be potentially useful for developers of the numerical codes used to simulate various types of the complex MHD phenomena.


Introduction
One of the pre-requisites to be able to deal with advanced physical transport phenomena involving the magnetohydrodynamics (MHD) interactions is to have a well-validated and numerically efficient computer code. This still poses a quite challenging task due to a lack of advanced experimental studies that can provide detailed insights into the flow and electromagnetic parameters that can be used to validate computer codes. The essence of the MHD phenomena is usually associated with a flow of highly electrically conducting liquid metals, which are, due to their non-transparency, notoriously difficult to study with standard laser-based optics diagnostics tools.
To validate MHD numerical models, we have to rely on analytical solutions that are based on significant simplifications. In the present manuscript, we are revisiting and proposing an extensive list of possible benchmark cases available in the open literature dealing with various aspects of the MHD phenomena. One of the simplest numerical MHD benchmarks is a fully developed laminar channel, duct, or pipe flow subjected to a uniform magnetic field of different orientations, for which an exact analytical solution exists, Hartmann and Lazarus (1937), Shercliff (1953). The effects of the non-uniform longitudinal magnetic field on a laminar flow of electrically conducting fluid in a pipe were recently numerically simulated in Feng et al. (2015). The open-source computer code OpenFOAM was used and good agreement was obtained between simulations and experiments. The MHD flow in a duct with very thin electrically conducting walls was presented in Tao and Ni (2013). Instead of fully resolving the wall region, a special type of boundary conditions was applied at the wall/fluid interface that takes into account a finite wall conductivity, as proposed in Walker (1981). It should be noted that this approach can be applied only for a very thin wall thickness and small conductance ratios.
Fusion engineering and technology-related research include numerous topics dealing with the MHD phenomena. Smolentsev et al. (2015) provided an extensive review of MHD codes for fusion applications and selected benchmark problems of importance for fusion applications. The proposed benchmarks covered a series of 2D and 3D steady and developing MHD flows in both laminar and turbulent regimes, and the final case also included effects of thermal buoyancy. Gajbhiye et al. (2018) validated their general-purpose solver by analyzing the free convection in a cubical enclosure under a uniform magnetic field and the electro-magnetically driven flow in a toroidal duct. The commercial ANSYS-CFX finite-volume based code was used to simulate a watercooled lithium lead (WCLL) breeding blanket module subjected to a strong uniform magnetic field, Tassone et al. (2017). The commercial multi-physics finite-element code COMSOL was successfully applied to simulate transient natural convection phenomena under influence of the imposed uniform magnetic field, Sahu and Bhattacharyay (2018). Validation of the multi-phase MHD flows is a challenging topic. The number of validation studies dealing with multi-phase MHD phenomena is significantly smaller compared to single-phase MHD phenomena. The analytical solutions for the multi-phase MHD situations are very scarce. One of the recently proposed analytical solutions for a 2D multi-phase MHD flow is presented in Righolt et al. (2016), where the elevation of the liquid-metal/air interface due to the presence of an imposed magnetic field is analytically solved. Numerical simulations of a rising bubble in the liquid metal subjected to an external homogeneous magnetic field of different strengths were studied in Shibasaki et al. (2010). The finite-difference code was used and the terminal bubble velocity dependency on the strength of the imposed magnetic field was analyzed.
Finally, the turbulent MHD phenomena require a special solving strategy due to the necessity to properly capture boththe flow and electromagnetic instabilities. The presence of the fluctuating Lorentz force requires a proper adaptation of the RANS-type of turbulence models (Kenjereš and Hanjalić, 2000;Kenjereš et al., 2004) or applications of the eddy-resolving simulation techniques such as Direct Numerical Simulations (DNS) or Large Eddy Simulations (LES), Kenjereš (2018). Krasnov et al. (2008) compared different sub-grid scale models for the MHD LES channel flow and demonstrated ability of the dynamic Smagorinsky model to properly predict the influence of the imposed magnetic field. Chaudhary et al. (2010) used DNS and analyzed how the increasing strength of a transverse magnetic field could influence the turbulence in the square duct flow. Mao et al. (2017) simulated the MHD flow in the insulated squared duct with different sub-grid scale models and compared data with the DNS results from the previous research of Chaudhary et al. (2010). Additionally, Mao et al. (2017) varied the Hartman number, showing how the turbulence is being suppressed by the imposed magnetic field.
The main goal of the present study is to obtain and validate results from our newly developed OpenFOAM solver over a range of various magnetohydrodynamic flows, and based on these findings, to propose an extensive numerical MHD benchmark, which can be potentially useful for developers of the computer codes for simulations of the MHD phenomena. We are primarily focusing on the influence of the finite electric conductivity of surrounding walls and the multiphase aspects of the MHD phenomena. We have analyzed the following situations: (i) a laminar duct flow with finite conductivity of surrounding walls, (ii) a laminar back-step flow, (iii) a shallow 2D multi-phase cavity, (iv) a rising bubble in the liquid metal, and, finally, (v) a turbulent duct flow with conducting walls. For all mentioned cases we performed a detailed comparative assessment against available analytical solutions or/and numerical results presented in the literature.

Governing equations for a single-phase MHD
We consider an incompressible electrically conductive fluid with liquid metal properties. The fluid is affected by the imposed external (constant) magnetic field through the Lorentz force. Conservation of mass and momentum are used to describe the MHD flow (under the assumption that the imposed magnetic field is known), and are written as: where U is velocity, p is pressure, ν is the kinematic viscosity, ρ is density, J is the current density and B is the imposed magnetic field. In the momentum equation, the MHD interactions are accounted for through the Lorentz force term. In addition to the velocity and pressure, also the current density (J) needs to be calculated. For the one-way coupled MHD phenomena, i.e. when the following conditions are valid where Re m is the magnetic Reynolds number, Pr m is the magnetic Prandtl number, L is the characteristic length and λ is the magnetic diffusion, ν is the kinematic viscosity, the Ohm's law for a moving conducting fluid is used where σ is the electrical conductivity of the fluid. By imposing the divergence-free current density condition in the Ohm's law, i.e.
the final Poisson's equation for the electric potential (ϕ) is obtained and can be written as In addition to Re m and Pr m (given in Eqn. (3)), the hydrodynamic Reynolds and Hartmann number are used as typical MHD nondimensional parameters:

Governing equations for a multi-phase MHD: volume of fluid method
In the current study, the Volume of Fluid (VOF) method is applied to the multi-phase MHD flow simulations. In addition to the Lorentz force, also the surface-tension and gravitational forces need to be included into the momentum equation:

∂U ∂t
where f g is the gravity force term, γ is the surface tension, k is the curvature of the interface (calculated as k = ∇⋅ ∇α |∇α| ), ν av is the phase averaged viscosity (calculated as ν av = α⋅ν 1 + (1 − α)⋅ν 2 , where '1' and '2' are phase indicators), ρ av is the phase averaged density (calculated as ρ av = α⋅ρ 1 + (1 − α)⋅ρ 2 ) and the volume fraction α is described by the following transport equation: where U r is the artificial compression velocity used for the interface sharpening, which is calculated as: where n f is the normal vector of the cell surface, ψ is the mass flux through the face, S f is the cell surface area, and C α is a coefficient that is used to control the interface thickness. There is no artificial interface compression when C α = 0. In order to control the spurious velocities which appear near the interface due to the sharp change of α, the volume fraction function is smoothed by the following Laplacian filter (Hoang et al., 2013;Mukherjee et al., 2018): where α is the resulting smooth volume fraction function, while subscripts c and f indicate the cell center and cell face, respectively. Using the smooth function α in Eqn. (9), instead of the original function α will suppress these parasitic velocities. In the current study, the filter (11) is applied twice for each time step.

The eddy-resolving MHD turbulence: large Eddy simulation
Turbulence modeling is performed by the Large Eddy Simulation (LES) method which is a good alternative to a more computationally expensive Direct Numerical Simulations (DNS). We apply the spatial filtering operation to Eqn. (2), which finally can be written as:

∂U ∂t
where ('') indicates the spatially filtered value and (τ sgs ) is the sub-grid scale stress tensor, and p * = p + 1 3 τ ′ I. The eddy viscosity concept is applied for the closure of the subgrid stress tensor as: where S is the modulus of the strain rate tensor, ν sgs is the subgrid scale turbulent viscosity, and C s is Smagorinsky coefficient. In the present work, we have adopted the dynamic approach to locally estimate values of the Smagorinsky coefficient, Lilly (1992), as follows: where Δ is the first filter (calculated as Δ = (Δ x Δ y Δ z ) 1/3 ), Δ is the second filter (calculated as Δ = 2Δ) and '〈…〉' means the local spatial averaging over the cell faces.

Numerical details 2.4.1. The conservative form of the Lorentz force
The additional Lorentz force in the momentum equation is traditionally treated in a non-conservative way (i.e. by applying the volume integration of the source term). This can potentially lead to significant numerical errors, especially for flow regimes with high Hartmann numbers. Similarly, the total electric current density must be conserved too. Both of these requirements are achieved through the application of the Four Steps Projection Method (FSPM) proposed by Ni et al. (2007), which can be summarized through the following four steps: 1. Calculate the magnetic flux at cell faces: where the cell -face electric conductivity (σ f ) is calculated by applying the harmonic average between different phases, and (S f ) is the cell surface area vector. 2. Use Eq. (17) to solve the discretized electrical potential equation and find electric potential (ϕ) at the cell centers: where 'm' indicates the number of cell faces. 3. Calculate the current density flux at cell faces using the surfacenormal gradient of electric potential (ϕ): where (J n ) is the cell face normal component of the current density. 4. Finally, use the current density flux from Eq. (19) and calculate the fully conservative form of the Lorentz force as: where (Ω c ) is the volume of cell, (r c ) is the cell center distance vector and (r f ) is the face center distance vector.

Conjugate MHD: taking into account electric conductivity and thickness of surrounding walls
The finite electric conductivity and finite thickness of surrounding walls have a significant impact on the fluid flow. This is due to the effects of the current density transfer between a liquid layer and solid walls, which is directly influencing the intensity and direction of the local Lorentz force in the near-wall region. To include the fluid/wall interface effects, we have developed an approach similar to traditional conjugate heat transfer, but now instead of the heat flux transfer, we focus on the distribution of the electric potential and current density in both domains. Transport equations of the electric potential in liquid (L) and solid (S) wall domains can be written as: Note that the source term (the RHS of Eqn.(21) is absent for the solid wall domain. Along the fluid/wall interface ( Fig. 1), the conservation and continuity of the electric current density (J) needs to be kept. This is achieved by imposing following set of the boundary conditions at the interface: The electric current density in the computational cell center is calculated in the same manner for both liquid and solid part of where the harmonic average is used to interpolate the electric conductivity at the interface, needed for calculation of the current density flux at the cell faces (J n,f ).

The computer code
The integrated MHD solver, which includes all above-listed transport equations, for both single-and multi-phase MHD phenomena is based on the finite-volume open-source computer code OpenFOAM-extend 4.0, Weller et al. (1998). Coupling between pressure and velocity field is performed with the PISO algorithm, Issa et al. (1986).

Laminar duct flow with conjugate MHD
In the first test case, we address a laminar pressure-driven flow of an electrically conducting fluid in the rectangular duct subjected to a transverse magnetic field, Fig. 2. The duct has the square cross-section (where L-is the half-width), length of 20L and d s is the thickness of side-walls (Hartmann walls). The Reynolds number is kept constant at Re = 10 and Hartmann number is varied in the 0⩽Ha⩽10 4 range. At the inlet, a uniform velocity profile is imposed. At all walls, the no-slip velocity boundary conditions are applied. At the outlet, a zero-pressure boundary condition is imposed. The uniform transverse magnetic field is imposed. To deal with the finite-thickness surrounding walls, we introduce characteristic wall conductance parameter, defined as: Three types of electric boundary conditions for the walls perpendicular to the magnetic field (Hartmann walls) are considered: (i) arbitrary conductive walls with varied wall conductance parameter (0.005⩽C d ⩽40), (ii) fully electrically insulated walls (∂ϕ/∂n = 0 and d S = 0), and finally, (iii) fully conductive walls (ϕ = 0 and d S = 0). The walls parallel to the magnetic field (Shercliff walls) are considered as electrically insulated for all cases.
Although the final steady-state results are validated against analytical solutions, the solution procedure is performed in a time-dependent mode. This time-dependent approach is not numerically efficient, but our final goal is to have a well-validated solver able to simulate MHD phenomena in transient and turbulent flow regimes, so we adopted a time-dependent solution approach for all benchmark cases presented here. The second-order central difference scheme (CDS) is applied for both convective and diffusive terms of discretised momentum equation, whereas the second-order backward scheme is used for time integration.
For all simulations the same hexahedral non-uniform orthogonal mesh is used with (Nx× Ny × Nz = 80 × 100 × 100) fluid control volumes for the fluid domain and (Nx × Ny × Nz = 80 × 10 × 100) solid for the solid domain, respectively. In making the spatial distribution of the nonuniform mesh, special attention is devoted that characteristic Hartmann and Shercliff boundary layers (with a typical thickness of δ Ha = L/Ha and δ Sh = L/Ha 1/2 ) are properly resolved. This is achieved by placing between 5 and 10 control volumes with a typical grid expansion ratio of 1.14 in the region bounded by the wall and the edge of the boundary layer (at δ Ha ).
Contours of the calculated streamwise velocity and electric potentialafter reaching steady state in the center of the duct (x = 10L) -are shown in Fig. 3. For the MHD neutral case (Ha = 0) the velocity exhibits a typical symmetric parabolic-like distribution, Fig. 3(a). By imposing the transverse magnetic field (Ha = 100) and by keeping all duct-walls  electrically insulated, a flattening of the velocity distribution occurs in the central part of the duct, whereas thin Hartmann boundary layers are generated along opposite vertical walls, Fig. 3(b). Next, by keeping the same strength and direction of the imposed magnetic field, and by changing electric properties of the vertical walls from fully insulated to walls with a finite thickness and conductivity (i.e. C d = 0.1), we observe a dramatic reorganization of the velocity with peaks in the proximity of the Shercliff walls, Fig. 3(c). Finally, by making Hartmann walls fully electrically conducting (C d →∞), the velocity distribution with two peaks is still present, Fig. 3(d). The electric potential contours exhibit close to a linear distribution in the vertical direction for electrically insulated and finite-conductivity Hartmann walls, Fig. 3(f) and (g). In contrast to this, the perfectly electrically conducting Hartmann walls impose almost a uniform distribution in the central part of duct, Fig. 3 (h).
The numerical solutions are compared next against the following analytical solutions: (1) Shercliff's solution for the electrically insulated   walls, Shercliff (1953), (2) Hunt's solution for the electrically fully conductive walls, Hunt (1965), and (3) Sloan's solution for the walls with the arbitrary electrical conductivity and thickness, Sloan and Smith (1966). For all simulated cases, an excellent agreement between present numerical simulations and analytical solutions is obtained, confirming an adequate implementation and validation of the conjugate MHD solver, Fig. 4.
To illustrate the sensitivity of the numerical solution, we perform a mesh dependency study with three mesh levels: (i) the coarse mesh

The 2D MHD laminar back-step flow
Next, we consider the two-dimensional backward-facing step flow in a laminar flow regime subjected to a uniform vertical magnetic field, Fig. 7. In contrast to the previous case, this configuration is expected to produce a more complex flow pattern with a well-defined recirculation region in the lower part of the domain. The channel height is L and its length is 15L. The lower and upper boundaries of the channel are no-slip walls. The upper half of the left boundary is the inlet, while the lower half is the solid wall. The inlet velocity is defined as: For the right boundary, a simple zero-gradient condition is imposed. All walls are treated as perfectly electrically insulated. The simulation domain and all boundary conditions are selected such that they match exactly the numerical study of Mramor et al. (2014), who applied a MHD extension of the Local Radial Basis Function Collocation Method (LRBFCM). The entire simulation domain is represented by an orthogonal numerical mesh with (Nx × Ny = 600 × 50) = (3 × 10 4 ) total control volumes. Two values of the Reynolds number are simulated (Re = 300 and 800, where Re = u x L/ν) over a range of Hartmann numbers (0⩽Ha⩽50). The second-order linear upwind differential scheme is used for convective terms, the second-order central differencing scheme (CDS) is used for diffusion terms, and the second-order backward scheme for the time-integration. The contours of the non-dimensional streamwise velocity (u/u 0 where u 0 = (u x )| x=0 , i.e. the inlet integrated velocity profile), at Re = 800 and different strengths of the imposed magnetic field (Ha = 0, 5, 10 and 50) are shown in Fig. 8. It can be seen that with a magnetic field increase, the recirculation length reduces, and flow becomes much more uniform. At Ha = 0, two large recirculation regions along the upper and lower walls are generated. With Ha increase, the recirculation region along the upper wall disappears, while the recirculation long the lower wall is still present, but its length is significantly reduced. This reduction of the recirculation region is further illustrated in zoom-in plots, where we superimposed contours of the streamwise velocity and streamlines, as shown in Fig. 9. At the highest value of Ha = 50, the recirculation can be observed only in a very small region attached to the lower part of the inlet plane. A   Fig. 7. The sketch of the simulation domain of the 2D laminar MHD back-step test case.  comparison of obtained profiles of horizontal (u/u 0 ) and vertical (v/u 0 ) velocity components at the exit plane with values presented in the literature (Mramor et al., 2014), are shown in Fig. 10. It can be seen that the horizontal velocity profiles become more flat with the magnetic field increase for both Reynolds numbers. The vertical velocity component almost completely disappears at higher values of Ha. A very good agreement between the present profiles and results from the literature (Mramor et al., 2014) is obtained for all presented cases. To demonstrate that the obtained results at present mesh of (Nx × Ny = 600 × 50) (M2) (3 × 10 4 CVs) are grid independent, one coarser (Nx × Ny = 300 × 25) (M1) (0.75 × 10 4 CVs) and one finner (Nx × Ny = 1200 × 100) (M3) (1.2 × 10 5 CVs) numerical mesh are generated, and results are compared in Fig. 11. A good agreement between different mesh levels is obtained, with a slight overprediction of the local maxima of the non-dimensional vertical velocity (v/u 0 ) at y/L = 0.7 for the coarse mesh.

The multi-phase two-dimensional shallow cavity flow with MHD
The first example of the MHD multi-phase test case is a shallow cavity subjected to combined effects of the imposed non-uniform mag-netic field and electric potential difference. The two-dimensional cavity with characteristic length L and partially filled with the electrically conductive liquid (where d is the liquid layer height and d≪L) is shown in Fig. 12. The upper part of the cavity if filled with air (σ air = O (10 − 15 ) S/m, i.e. negligible electric conductivity). The external magnetic field is aligned with the negative z-direction (perpendicular to the cavity) and its linear distribution is defined as: where α b = 0.1 defines a distribution parameter. The no-slip velocity boundary condition is imposed at all walls (i.e. bottom and side-walls). The gravity force is aligned with the negative y-coordinate direction. The side-walls are kept at constant (but different) electric potential (ϕ 1 = − 1 2 Δϕ, ϕ 2 = 1 2 Δϕ, where Δϕ is the imposed electric potential difference). The bottom wall is perfectly electrically insulated (∂ϕ/∂n = 0 and C d = 0). Because of the imposed magnetic field and electric potential difference, the generated Lorentz force within the fluid will drive the flow. This fluid motion will be opposed by a joint combination of the viscous, gravitational, and surface tension forces. To account for additional free-surface related physical mechanisms, the following set of non-dimensional parameters is introduced, Righolt et al. (2016): In addition to the redefined Reynolds (Re * ) and Hartmann (Ha * ), also the Bond (Bo * ) and capillary (Ca * ) numbers are introduced. The characteristic non-dimensional velocity is calculated as: Because of the large number of possibilities based on the various combinations of characteristic non-dimensional numbers, in the present work we kept constant Re * = A and Ha * = 1, while we change Bo * and Ca * . We also kept the identical aspect ratio of the domain, A = d/L = 0.1. The two-dimensional orthogonal, non-uniform mesh (Nx × Ny = 50 × 200) with rectangular control volumes is used. The central differencing scheme (CDS) is used for the diffusive and convective terms of transport equations. The time-integration is performed with the secondorder backward scheme. For this particular case, the different values of the interface compression coefficient (0⩽C α ⩽1) did not have any significant impact on the obtained solutions due to a smooth free-surface deformation. The local variation of the resulting Lorentz force generates the flow of electrically conducting fluid (initially at rest) in the lower part of the cavity with characteristic elevation of the free-surface, as shown in Fig. 13.
This non-dimensional vertical elevation (h/d) of the free-surface, as a function of Ca * and Bo * numbers, is shown in Fig. 14. It can be seen that an excellent match between the present numerical results (CFD) and analytical solutions is obtained for all calculated cases. Note that a vertical elevation of the free-surface increases with an increase in both Ca * and Bo * . The horizontal profiles of the non-dimensional horizontal (u/u * ) and vertical (v/u * ) velocity in the proximity of the left-wall are shown in Fig. 15(a) and (b), respectively. The vertical profile of the nondimensional horizontal velocity at the central vertical line is shown in Fig. 15(c). Again, an excellent agreement between the present simulation (CFD) and analytical solution from the literature (Righolt et al., 2016) is obtained, proving the capability of the MHD multi-phase solver. To confirm that the presented solutions are grid independent, we analyzed the non-dimensional free-surface elevation (h/d) for three mesh sizes: (i) the coarse mesh (M1) (N x × N y = 25 × 100), (ii) the intermediate (previously presented results) mesh (M2) (N x × N y = 50 × 200), and (iii) the fine mesh (M3) with (N x × N y = 100 × 400). A good agreement between results at different mesh resolutions confirms the full mesh convergence of the presented results, Fig. 16.

The 3D rising gas bubble in liquid metal subjected to a longitudinal magnetic field
A rising gas bubble (with an initial diameter d b = L/2) is submerged into the liquid metal confined in the 3D rectangular box (with height 3L, width and depth L) is analyzed next, Fig. 17. This test case is based on the study of Shibasaki et al. (2010). All boundary surfaces are electrically insulated walls (∂ϕ/∂n = 0,C d = 0) with imposed no-slip boundary conditions. The external magnetic field is aligned with the y-coordinate and the gravity is oriented in the opposite direction. The problem is fully defined with the following set of non-dimensional parameters: where G is the Galilei number, Γ is the Tension number, and subscripts (G) and (L) indicate the gas and liquid phase, respectively. The nondimensional velocity, pressure, and time are defined as: We kept constant G = 4⋅10 4 , Γ = 2⋅10 6 and varied 0⩽Ha⩽200 to study the influence of the magnetic field strength on the rising bubble behavior. The electrical conductivity ratio is σ G /σ L = 2.49⋅10 − 7 . The orthogonal mesh is created with the mesh size (N x × N y × N z = 60 × 180 × 60), identical to the mesh used in Shibasaki et al. (2010). The second-order linear-upwind scheme is used for the convective terms in both momentum and volume fraction equations, whereas the backward scheme is used for time integration. Because of a sharp jump of the electrical properties at the phase interface, we have applied the harmonic interpolation scheme for the electric conductivity. For this case, the interface compression coefficient (C α ) had stronger effect on the final shape of the rising bubble. The selected value of C α = 0.1 proved to be a good choice for both multi-phase benchmarks presented here. The obtained characteristic bubble shape, current density streamlines, contours of the vertical velocity and pressure in the central vertical plane at an arbitrary time instant t/t * = 0.02 and for Ha = 50, are shown in Fig. 18. The current density streamlines form close loops around the bubble with opposite directions above and below the bubble. The velocity contours portray an updraft region in the center of the domainabove and below the bubble, whereas the down-drafts are generated along the side walls. Contours of the pressure exhibit almost linear distribution in the vertical direction, with small deviations in the proximity of the bubble surface. It can be seen that the resulting shape of the bubble strongly depends on the imposed magnetic field strength, Fig. 19. The higher Ha leads to the bubble stretching in the direction of the imposed magnetic field (y-direction) and to a reduction of its rising velocity. We compare our results with a numerical study of Shibasaki et al. (2010) who applied the finitedifference (FDM) multi-phase MHD code. Comparison of the computed terminal velocity for different values of Ha is shown in Fig. 20. After an

initial slight increase in the terminal velocity for intermediate values of
Ha < 50, a gradual decrease is obtained with a further increase of the imposed magnetic field. The agreement between the current simulations and data presented in Shibasaki et al. (2010) is good up to Ha = 50. After reaching this peak value, larger differences are observed, but qualitatively similar trends are observed. Differences for larger values of Ha number can be partially explained by the use of different discretization approaches (the present finite-volume vs. finite-difference of Shibasaki et al. (2010)), the application of different convective schemes (the present second-order linear-upwind vs. the third-order UTOPIA scheme of Shibasaki et al. (2010)), as well as due to the absence of the mesh-dependency study of Shibasaki et al. (2010)). We also performed additional simulations with a second-order quadratic-upwind scheme for convective terms in momentum equations, but this resulted in marginal differences of rising velocity (less than 1%) compared to the linearupwind scheme (see Table 1). Finally, we complete a mesh-dependency study for two different Hartmann numbers Ha =50 and 200, and three meshes: (i) the coarse mesh (M1) (Nx × Ny × Nz = 30 × 90 × 30) = (0.081 × 10 6 ) total CVs, (ii) the present mesh (M2) (Nx × Ny × Nz = 60 × 180 × 60) = (0.64 × 10 6 ) total CVs and (iii) the fine mesh (M3) (Nx × Ny × Nz = 120 × 360 × 120) = (5.1 × 10 6 ) total CVs. Results in Table 2 demonstrate that the finest mesh (M3) provides the best agreement with the reference data. However, the difference in terminal velocity values between intermediate (M2) and fine (M3) mesh is only 1%, while the total number of CVs is four times larger. Based on this small difference, we conclude that results are grid independent already at the mesh (M2).

A conjugate MHD duct flow in a fully developed turbulent regime
The final test case is a conjugate MHD square duct flow in a fully developed turbulent regime. The duct height and width are L, and its length is 16L. The imposed magnetic field is aligned with the y-axis and perpendicular to the flow direction, Fig. 21. The periodic boundary conditions are imposed in the streamwise (x-coordinate) direction. All other surfaces are walls with imposed no-slip velocity boundary conditions. The lower and upper walls (Shercliff walls) are fully electrically insulated (∂ϕ/∂n = 0 and C d = 0). The front and back walls (perpendicular to the imposed magnetic field -Hartmann walls) are considered to have three different types of electric boundary conditions: (i) the finite conducting walls with the wall conductance parameter C d = (σ S d S )/(σ L L) = 0.05, (ii) fully electrically insulated walls, and (iii) fully conductive walls (ϕ = 0 and C d →∞). We apply the wall-resolving dynamic large-eddy simulation (LES) approach. The numerical mesh contains (N x × N y × N z = 240 × 120 × 120 ) fluid = (3.456 × 10 6 ) CVs and (N x × N y × N z = 240 × 12 × 120 ) wall = (0.3456 × 10 6 ) CVs in the fluid and wall regions, respectively. The non-dimensional mesh parameters are Δy + wall = Δz + wall ≈ 0.6, Δy + core = Δz + core ≈ 6 and Δx + ≈ 25. The central differencing scheme (CDS) is used for spatial discretization and the second-order backward scheme for temporal discretization. The flow is defined with the following set of the non-dimensional parameters: Re = 5602 and Ha = 21.2. Furthermore, a simulation with Ha = 0 is performed in order to provide a comparison with the non-MHD neutral case. The selected value of the Reynolds number assures that a fully developed turbulence is generated and maintained. All simulations are statistically averaged over at least 100 flow-through times, and the spatial averaging procedure is applied to accelerate the convergence of the flow statistics. The instantaneous coherent structures colored by streamwise velocity for various Ha, and Hartmann walls conductivities (expressed through the wall-conductivity parameter, C d ), are shown in Fig. 22. Under the action of the imposed transverse magnetic field, by changing the electric properties of the walls, the coherent structures start to be suppressed in the proximity of Hartmann walls, as seen from the side-views of the duct shown in Fig. 22(a)-(c).
Contours of the long-term time-averaged streamwise velocity, turbulent kinetic energy, and electric potential, for various wall conductivities, are shown in Fig. 23. Starting from a symmetrical distribution, contours of the mean streamwise velocity start to be suppressed in the direction of the imposed magnetic field (y-direction). This behavior is caused by the reorganization of the electric current density streamlines. In contrast to the fully closed current loops in the fluid region (for electrically insulated walls), a finite electric conductivity of walls makes that current density loops also enter these regions, causing significant changes in resulting Lorentz force components in the y-and z-directions, respectively.
The contours of the turbulent kinetic energy portray the reorganization from fully symmetrical distributions for the neutral case with characteristic peaks in the proximity of duct walls, Fig. 23(e), to the nonsymmetrical distributions for non-insulated walls, Fig. 23(f)-(h). It can be seen that the levels of turbulent kinetic energy are suppressed in the proximity of Hartmann walls for the fully insulated and walls with finite conductivity, Fig. 23(f)-(g). At the same time, distributions of the   Note that scaling factor of 2× is applied in the xdirection to provide a more compact view. turbulent kinetic energy along Shercliff walls are just slightly affected by changing Hartmann walls conductivity, Fig. 23(f)-(g). Interestingly, for the fully conducting walls, the turbulence kinetic energy along Shercliff walls is reduced in comparison to values along Hartmann walls, Fig. 23 (h). The contours of the electric potential also illustrate significant changes between the fully electrically insulated Fig. 23 and fully conducting Hartmann walls Fig. 23(l), with the latter exhibiting significantly more pronounced non-uniform distribution in the vertical direction. Next, we move to a more detailed comparison of the characteristic long-term time-averaged first-and second-order statistics profiles, Figs. 24 and 25. The present results obtained with the dynamic LES approach are validated against two reference studies: Gavrilakis (1992) who simulated an MHD neutral turbulent duct flow, and Chaudhary et al. (2010) who simulated MHD turbulent duct flow at Ha = 21.2 with fully electrically insulated wallsboth using the fullyresolving Direct Numerical Simulations (DNS) approach. The timeaveraged mean streamwise velocity profiles in the proximity of Hartmann and Shercliff walls are shown in Fig. 24(a) and (b), respectively. It can be seen that a very good agreement between the present and DNS results from the literature is obtained at both locations. The profiles of the non-dimensional rms values of the fluctuating streamwise velocity   24. The long-term time-averaged non-dimensional streamwise velocity (the semi-log plots of U + vs. y + and z + ) profiles (where U + = U/U τ ,y + = yU τ /ν, and) in the proximity of Hartmann (a) and Shercliff wall (b), respectively. along the identical locations reveal an interesting behavior, Fig. 25. In the proximity of Shercliff walls, an increase of the wall conductivity produced a gradual decrease of the rms values, Fig. 25(b). In contrast to this behavior, the distributions in the proximity of Hartmann walls indicate an initial suppression for the fully electrically insulated walls, followed by an increase for the fully conducting walls, Fig. 25(a). Again, a good agreement with available DNS references (for the non-MHD situation and the MHD case with fully insulated walls) is obtained confirming suitability of here used dynamic LES approach.

Summary and conclusion
We have presented a comprehensive numerical benchmark study addressing a range of single-and multi-phase one-way coupled MHD flows. The single-phase cases included the conjugate MHD flows in ducts with varied electric conductivity of the wallin both laminar and turbulent flow regimes, and the laminar back-step flow subjected to a transverse magnetic field. The multi-phase cases covered a twodimensional MHD cavity and a rising bubble in a liquid metal flowsboth simulated with the volume of fluid (VOF) approach. We have implemented an extended set of MHD transport equations in the opensource code OpenFOAM. Our particular focus was to extend the existing set of MHD benchmarks and to provide a detailed comparison with similar studies in the literature. We also proposed a novel methodology and benchmark for a conjugate MHD in a turbulent duct flow with an arbitrary wall conductivity (expressed in terms of the wall conductance parameter). For the multi-phase flows, we have introduced a recently proposed analytical solution of a two-dimensional partially-filled cavity flow subjected to an external magnetic field. An excellent agreement was obtained for all cases for which analytical solutions are available. For considered test cases without analytical solutions, a very good agreement was obtained with available numerical studies from the literature. It is concluded that here developed and validated version of the computer code can be used for advanced fundamental and industrial/technological studies involving various aspects of the MHD phenomena.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.