# On Boundary Conditions for Lattice Kinetic Schemes for Magnetohydrodynamics
Bounded by Walls with Finite Electrical Conductivities

Eman Yahia<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Colorado Denver, 1200 Larimer street, Denver, CO 80204, U.S.A.
Kannan N. Premnath<EMAIL_ADDRESS>Department of Mechanical
Engineering, University of Colorado Denver, 1200 Larimer street, Denver, CO
80204, U.S.A.

(August 28, 2024)

###### Abstract

Magnetohydrodynamic (MHD) flow of liquid metals through conduits play an
important role in the proposed configurations for harnessing fusion energy,
and various other engineering and scientific problems. The interplay between
the magnetic fields and the fluid motion gives rise to complex flow physics.
Moreover, the electrical conductivity of the bounding solid walls relative to
that of the fluid, in addition to the Hartmann number, has a major influence
on the structure of the flow field and the induced magnetic fields in such
cases. An effective approach to represent such effects is via the Shercliff
boundary condition for thin conducting walls relating the induced magnetic
field and its wall normal gradient at the boundary via a parameter referred to
as the wall conductance ratio [1]. Within the framework of the lattice
Boltzmann (LB) method for MHD, a lattice kinetic scheme (LKS) involving a
vector distribution function for the magnetic fields was proposed by Dellar
[2]. However, the prior LB studies only considered boundary implementations
for the limiting special cases involving either insulating or perfectly
conducting walls without accounting for the finite conductivity effects of the
container walls. In this paper, we present two new boundary schemes that
enforce the Shercliff boundary condition in the LB schemes for MHD. It allows
for the specification of the wall conductance ratio to any desired value based
on the actual conductivities and length scales of the container and the wall
thicknesses and can be made to vary freely between the special cases of zero
(for insulating walls) and infinity (for perfectly conducting walls). One
approach is constructed using a link-based formulation involving a weighted
combination of the bounce back and anti-bounce back of the distribution
function for the magnetic field and the other approach involves an on-node
moment-based implementation. Both the boundary schemes are parameterized by
the relaxation time of the distribution function of the magnetic field and the
wall conductance ratio, and prescribe the closure relations for the incoming
distribution functions from the walls. Moreover, their extensions to
representing moving walls are presented. The boundary schemes are validated
for canonical body force driven or shear driven MHD flows for a wide range of
the values of the wall conductance ratio. The accuracy and second order grid
convergence of these new and more general boundary schemes for LKS for MHD are
demonstrated.

Lattice Boltzmann method, Magnetohydrodynamics, Boundary conditions, Wall
electrical Conductivities, Link-based approach, Moment-based approach

###### pacs:

47.11.Qr,05.20.Dd,47.27.-i

††preprint: PREPRINT

##  I Introduction

The modern study of magnetohydrodynamics (MHD) can be traced back to the
pioneering works of Hartmann [3, 4] and Alfvén [5]. The former performed an
analytical and experimental investigation of liquid metals in channels, while
the latter predicted the existence of the MHD waves, which was then confirmed
experimentally [6]. Much of the early investigations centered on MHD problems
related to astrophysics and the geodynamo [7, 8]. On the other hand, there
have been increasing number of technological problems involving MHD in power
generation and materials processing applications [9], and in systems such as
tokamaks for the demonstration of the nuclear fusion process [10].

Since a successful harnessing of fusion energy would represent a major
breakthrough advance in meeting the future energy needs with essentially no
carbon footprint, there have many ongoing research endeavors in this
direction, including the experimental ITER project, aimed to overcome the
associated challenges. From a fluid mechanics perspective, fusion energy
systems are of great current interest as they involve a variety of
configurations involving the flow of liquid metals under high magnetic fields
and subjected to different conditions thereby representing rich avenues for
studying and understanding the multiphysics associated with wall bounded MHD
flows (see e.g., [11, 12, 13, 14, 15]). For example, the transport of the
lead-lithium alloy and confined by solid wall structures is being considered
both as a breeder and coolant as part of the liquid blanket module in the
proposed fusion reactor configurations.

The fluid motions of the liquid metals under imposed magnetic fields are
represented by the Navier-Stokes (NS) equations augmented with a locally
varying Lorentz force that arises from the interaction between the
electrically conducting fluid in motion and the magnetic field. The evolution
of the components of the magnetic field is based on a combined representation
of the different laws given by the Maxwell’s equations of electromagnetism and
referred to as the magnetic induction equation, which is, in turn, influenced
by the fluid motions [16, 17, 18]. Taken together, these macroscopic equations
for MHD, which are briefly summarized in the next section (see Sec. II), are,
in general, strongly coupled system for the spatio-temporal variations of the
velocity and magnetic fields.

An important element associated with solving the MHD equations is the
representation of the boundary conditions. While the velocity field is
subjected to the well-known no-slip conditions at the boundary, the boundary
conditions for the magnetic field requires some further consideration. Since
the metallic fluid in motion is confined by a solid wall, whose material, in
general, can have a finite electrically conductivity, the effect of the latter
needs to be represented on the former consistently. Indeed, the electrical
properties of the solid wall can significantly affect the structure of both
the magnetic and velocity fields for the coupled MHD system (see e.g., [10,
19, 20]). One approach to accommodate this would be to solve the magnetic
induction equations within both the liquid and the solid and then impose the
appropriate conjugate continuity conditions at their boundary locations.
However, due to the possibly large disparities between the electrical
conductivities of the liquid and solid, such an approach becomes stiff and
computationally challenging to solve.

An alternative and effective strategy would be to impose the effect of the
solid with its finite conductivity only as a boundary condition under a thin
conducting wall assumption first proposed by Shercliff [1], which avoids the
need to directly solve for the magnetic field variations within the solid.
Essentially, this approach, whose derivation is given for completeness in Sec.
III, involves relating the induced magnetic field and its wall normal gradient
at the boundary via a parameter referred to as the wall conductance ratio
$c_{w}$ [10]. Some early analytical studies that demonstrate the influence of
this parameter on certain wall-bounded canonical MHD problems are presented in
Refs. [21, 22, 23]. Moreover, over the years, the Shercliff boundary condition
accommodating the electrical properties of the wall have been widely used in
various numerical investigations that simulate MHD within confining
geometries, including those related to the fusion engineering applications –
for example, see Refs. [24, 25, 26, 27, 28, 29, 15]. Indeed, the main
characteristic electro-magnetic parameters in such MHD problems are the
Hartmann number, which is related to the relative effect of the imposed
magnetic field strength and defined later, and the wall conductance ratio.
Thus, a key feature of any effective numerical scheme for realistic liquid
metal MHD in wall bounded conduits requires a strategy that accommodates the
finite conductivity of walls via the Shercliff boundary condition.

The lattice Boltzmann (LB) method [30, 31, 32, 33] is a mesoscopic
computational technique for fluid dynamics and other transport problems, and
has seen remarkable success during the last three decades [34, 35]. Emerging
from the lattice gas automata with further contributions from certain key
elements of the kinetic theory, it has evolved into a highly parallelizable
numerical algorithm that computes the distribution functions varying as a
result of a local collision step, which is followed by a lock-step advection
or streaming along certain discrete particle velocity directions. Such
discrete velocities, which are referred to as the lattice, respect certain
symmetric properties in order to recover the required invariance properties
associated with the macroscopic behavior of the fluid motions, such as the
isotropy of the viscous stress tensor. It may be noted that the standard
convention to represent a lattice set is $DdQq$, where $d$ represents the
number of spatial dimensions and $q$ denotes the number of discrete particle
directions. The effect of collision is often modeled via a relaxation process
that represents the approach of the distribution functions or their moments to
their equilibria at certain characteristic relaxation times. The boundary
conditions on the velocity field in simulating the NS equations, for example,
need to be specified in terms of the incoming distribution functions from the
boundary nodes. Since the boundary conditions have a major influence on the
numerical accuracy, stability and computational efficiency, the development of
the required closure relations has attracted significant focus and much effort
(see e.g., Ref. [36]). For representing the no-slip velocity boundary
condition, one of the most popular approaches is the link-based bounce-back
type boundary scheme and its variants (e.g., Refs. [37, 38]). As an
alternative in this regard, one could enforce a subset of the hydrodynamic
moments of the distribution functions directly on the lattice nodes (e.g.,
Refs. [39, 40, 41, 42, 43]).

The LB method has the following advantages compared to the more conventional
numerical techniques. First, its streaming step is linear, while its collision
step models all nonlinearity locally and are thus characterized by its simpler
implementations; on the other hand, given that the convective term of the
Navier-Stokes equations is nonlinear and nonlocal, their direct solution
techniques based on classical approaches require more cumbersome
implementation procedures. Moreover, the pressure field is computed locally
via an equation of state in the LB method, which avoids the need to solve the
time-consuming Poisson-type elliptic equation required by the conventional
methods. Second, the local algorithmic procedure of the LB approach lends
itself for a natural implementation on modern parallel computing systems for
efficient simulations of large size flow problems. Third, the boundary
conditions associated with complex geometries can be accommodated with ease in
the LB method via relatively simple rules on the particle distribution
functions on Cartesian grids, which circumvents the often time-consuming grid
generation process required with traditional approaches. And, fourth, the
coupling of the complex multiphysics phenomena with fluid motions can be
represented readily in the LB formulations via utilizing additional
distribution functions and specialized collision models, with their
implementations following the same general framework of the collide-and-stream
algorithmic steps.

The LB method has also been applied for simulating MHD. Some of the earlier
approaches [44, 45] in such cases were inspired from the developments in the
lattice gas automata [46]. These, along with other formulations [47, 48], led
to complicated algorithms for simultaneously computing the velocity and
magnetic fields, which compromised their implementation and efficiency.
Recognizing the features unique to the vectorial magnetic induction equation,
Dellar [2] proposed a separate evolution equation for the vector distribution
function, referred to as the lattice kinetic scheme (LKS) for representing the
dynamics of the magnetic field. On the other hand, the usual scalar
distribution function is solved for the velocity field. This strategy resulted
a significant simplification when compared to the earlier approaches LB
approaches for MHD. It was further extended by Pattison _et al_ [49], in
particular, in using grid refinement approaches for resolving the Hartmann
layers at high Hartmann number MHD flows relevant to fusion MHD problems and a
more general collision model for the velocity field solver for enhanced
stability. Furthermore, a moment system representation for the evolution of
the vector distribution function was proposed and analyzed by Dellar [50], who
also extended it for handling a class of MHD problems with current-dependent
resistivity [51]. Also, recently, De Rosis _et al_ [52] presented a one-stage,
simplified LBM for efficient simulation of MHD. Here, we mention that the
prior investigations involving LB schemes for MHD have provided relatively
limited focus on the boundary conditions for the vector distribution functions
for the magnetic field. Most of the studies treat the wall to be insulated – a
limiting special case situation, which was implemented using an anti-bounce
back scheme in Ref. [2], an extrapolation method in Ref. [49] and a moment
method in Ref. [43]. As such, none of the prior investigations seem to have
explicitly considered the finite electrical properties of the wall, which, as
mentioned above, can have significant influence on both the magnetic and
velocity field distributions within the fluid in motion, and hence are
essential in faithfully representing realistic wall-bounded MHD flows,
especially involving liquid metals.

This paper addresses the above mentioned limitation and constructs and
analyzes two new boundary schemes for the vector distribution function for MHD
by incorporating the actual electrical properties of the wall via the
Shercliff condition for thin conducting walls [1]. The first approach is
constructed using a link-based formulation, which, as we shall see, would be
based on a weighted combination of the bounce back and anti-bounce back of the
distribution function for the induced magnetic field. The second boundary
scheme developed for the conducing walls will be an on-node moment based
approach. The resulting closure relations for the incoming distribution
functions from the walls for both the boundary schemes will be shown to be
parameterized by the relaxation time of the distribution function of the
magnetic field and the wall conductance ratio. Both the schemes will be
further generalized for representing MHD involving moving walls. The resulting
boundary condition specification approaches can use the wall conductance ratio
based on the true electrical conductivity of the wall relative to that of the
fluid and also involving the actual length scale of the container relative to
the wall thickness. We will then test and validate our boundary formulations
for well-defined class of body force driven as well as shear driven MHD flows
for a wide range of the wall conductance ratio, and also investigate the order
of accuracy for grid convergence.

The outline of this paper is as follows. Section II presents the governing
macroscopic equations for MHD, which is then followed by a derivation of the
Shercliff thin wall boundary condition in Sec. III. The LB algorithms for
evolving a vector distribution function for the magnetic field and a scalar
distribution function for the velocity field are discussed briefly in Secs. IV
and V, respectively. Then, Secs. VI and VII present the derivations of the new
link-based and on-node moment-based boundary schemes, respectively, for the
magnetic induction equation in wall-bounded geometries. The results for
validation and grid convergence of the two boundary schemes for the magnetic
distribution function for the canonical body-force driven and shear-driven MHD
flows are discussed in Sec. VIII. Finally, a summary of the developments and
conclusions arising from this work as well as its future outlook are provided
in Sec. IX. Moreover, various additional details supporting the developments
of this work are presented in the appendices. For completeness, A discusses
the moment-based formulation for imposing the no-slip velocity condition, and
the two boundary schemes constructed in the main sections are extended for
moving walls in B. Finally, a detailed derivation of the analytical solution
of the shear-driven MHD flow between two parallel plates with contrasts in the
wall conductance ratios, which is utilized for validation of the proposed
boundary schemes, is presented in C.

##  II Macroscopic Equations of MHD

The governing equations for MHD is a coupled system of equations involving the
evolution of the density and momentum fields given by the standard Navier-
Stokes equations including the Lorentz force, and a magnetic induction
equation representing the variation of the magnetic fields. It may be noted
that the magnetic induction equation is obtained from the Maxwell’s equations
of electromagnetism by combining the Faraday’s law, Ampère-Maxwell’s law and
Ohm’s law under the constraints of the solenoidal magnetic field, i.e., no
magnetic monopoles [16, 17, 18, 10]. Thus, the MHD equations can be written as

$\displaystyle\partial_{t}\rho+{\bm{\nabla}}\cdot(\rho\bm{u})$
$\displaystyle=$ $\displaystyle 0,$ (1a)
$\displaystyle\partial_{t}(\rho\bm{u})+{\bm{\nabla}}\cdot(\rho\bm{u}\bm{u})$
$\displaystyle=$
$\displaystyle-{\bm{\nabla}}p+{\bm{\nabla}}\cdot\bm{\mathsf{T}}+\bm{F}_{Lorentz}+\bm{F}_{ext},$
(1b)
$\displaystyle\partial_{t}\bm{B}+{\bm{\nabla}}\cdot(\bm{u}\bm{B}-\bm{B}\bm{u})$
$\displaystyle=$
$\displaystyle\eta\nabla^{2}\bm{B},\quad\quad\eta=\frac{1}{\sigma\mu_{m}},$
(1c)

where $\rho$, $\bm{u}$ and $\bm{B}$ are the density field, velocity field, and
magnetic field, respectively, $p$ is the pressure, $\bm{F}_{Lorentz}$ is the
Lorentz force, $\bm{F}_{ext}$ is an external body force. Here,
$\bm{\mathsf{T}}$ is the viscous stress tensor given by
$\bm{\mathsf{T}}=\rho\nu\left[\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{\dagger}\right]$,
where $\nu$ is the kinematic viscosity of the fluid, and $\eta$ is the
magnetic resistivity, which is related to the electrical conductivity of the
fluid $\sigma$ and the magnetic permeability $\mu_{m}$. The Lorentz force
$\bm{F}_{Lorentz}$ appearing in Eq. (1b) reads as

$\bm{F}_{Lorentz}=\bm{J}\times\bm{B},$ (2)

where

$\bm{J}=\frac{1}{\mu_{m}}\bm{\nabla}\times\bm{B}$ (3)

is the local current density in the fluid. The above system of equations are
to be supplemented with appropriate boundary conditions at the interface
between the fluid and the solid walls. The velocity field satisfies the
standard no-slip condition. A treatment of the boundary condition for the
magnetic field that accounts for the finite wall conductivity is discussed
next.

##  III Shercliff Boundary Condition on the Induced Magnetic Field for Thin
Conducting Walls

We will now present a derivation of the of the boundary condition on the
induced magnetic field for thin conducting walls proposed by Shercliff [1].
Figure 1 illustrates a cross section of a fluid with an electrical
conductivity $\sigma$ bounded by a solid wall whose electrical conductivity is
$\sigma_{s}$, where $H$ is the characteristic length scale of the confinement
of the fluid and $\delta$ is the wall thickness. The boundary between the
fluid and solid is represented by $\Gamma$ with an outward unit normal
$\hat{n}$ and a unit tangent $\hat{t}$.

Figure 1: Illustration of a cross section of an electrically conducting fluid
bounded by a solid wall with a finite electrical conductivity for the
derivation of the Shercliff thin wall boundary condition.

Let $B$ denotes the cross-section normal component of the induced magnetic
field in the fluid, while $B_{s}$ represents a similar component within the
solid wall. On the interface $\Gamma$, the no-slip velocity condition, and the
continuity of the tangential component of the electrical field $E$, and of the
cross-section normal component of the induced magnetic field $B$ are
satisfied, which read as

$\bm{u}=\bm{u}_{s}=0,\quad\quad E_{t}=E_{t,s},\quad\quad
B=B_{s},\quad\quad\mbox{on}\;\;\Gamma.$ (4)

Using the Ohm’s law for the current density
$\bm{J}=\sigma(\bm{E}+\bm{u}\times\bm{B})$ in the tangential direction along
$\Gamma$, and by exploiting the conditions appearing in the no-slip condition,
i.e., $J_{t}|_{\Gamma}=\sigma E_{t}|_{\Gamma}$ or
$E_{t}|_{\Gamma}=J_{t}|_{\Gamma}/\sigma$, the continuity condition of the
electric field in the above (Eq. (4)) becomes

$\left(\frac{J_{t}}{\sigma}\right)_{\Gamma}=\left(\frac{J_{t,s}}{\sigma_{s}}\right)_{\Gamma}.$
(5)

Then, writing the tangential current density $J_{t}$ as a conjugation of the
normal derivative of magnetic field component $B$, i.e.,
$J_{t}=(1/\mu_{m})\partial B/\partial\hat{n}$, Eq. (5) can be rewritten as

$\left(\frac{1}{\sigma}\frac{\partial
B}{\partial\hat{n}}\right)_{\Gamma}=\left(\frac{1}{\sigma_{s}}\frac{\partial
B_{s}}{\partial\hat{n}}\right)_{\Gamma},$ (6)

It may be noted that the induced magnetic field $B$ in the fluid satisfies a
component of the magnetic induction equation given in Eq. (1c), while the
component $B_{s}$ in the solid obeys

$\bm{\nabla}\cdot(\eta_{s}\bm{\nabla}B_{s})=0,\quad\quad\eta_{s}=1/(\sigma_{s}\mu_{m}).$
(7)

To proceed further, following Shercliff [1], we now assume that the solid
conducting walls are relatively thin, i.e., $\delta\ll H$. This means that the
solution of the Laplace’s equation (Eq. (7)) for $B_{s}$ can be well
approximated by a function _linear_ in the spatial coordinates, analogous to
the solution of the Laplace’s equation within a membrane. Let $\xi$ be the
local coordinate normal to the wall (see Fig. 1). Since $\hat{n}$ is the
outward unit normal by definition, it follows that $\xi=-\hat{n}H$. Then, we
can prescribe an ansatz $B_{s}=c_{1}\xi+c_{2}$ and subject to
$B_{s}(\xi=0)=B_{s}|_{\Gamma}$ and $B_{s}(\xi=-\delta)=0$, whose solution
immediately yields $B_{s}(\xi)=B_{s}|_{\Gamma}(1+\xi/\delta)$ or
$B_{s}(\hat{n})=B_{s}|_{\Gamma}(1-\hat{n}H/\delta)=B|_{\Gamma}(1-\hat{n}H/\delta)$,
where the last expression arises from the continuity of the magnetic field at
the interface via Eq. (4). Then, the wall normal gradient of the induced
magnetic field within the solid wall evaluated at the interface becomes the
following:

$\left(\frac{\partial
B_{s}}{\partial\hat{n}}\right)_{\Gamma}=-B|_{\Gamma}\frac{H}{\delta}.$ (8)

Using the above expression (Eq. (8)), which essentially is only related to the
magnetic field of interest from the fluid side evaluated at the interface and
the inherent length scales of the problem, in Eq. (6), we obtain

$\left(\frac{1}{\sigma}\frac{\partial
B}{\partial\hat{n}}\right)_{\Gamma}=-\frac{1}{\sigma_{s}}\frac{H}{\delta}B|_{\Gamma}.$
(9)

Finally, defining

$c_{w}=\frac{\sigma_{s}\delta}{\sigma H}$ (10)

as the _wall conductance ratio_ parameter, which is a dimensionless quantity,
and using it in Eq. (9), we obtain the Shercliff boundary condition for the
induced magnetic field involving a thin conducting wall as [10]

$\frac{\partial
B}{\partial\hat{n}}+\frac{B}{c_{w}}=0\quad\quad\mbox{on}\;\;\Gamma$ (11)

given in terms of the unit normal $\hat{n}$ gradient. This can also be
rewritten in terms of the local spatial coordinate, via $\xi=-\hat{n}H$, as

$-\frac{\partial
B}{\partial\xi}+\frac{B}{c_{w}H}=0\quad\quad\mbox{on}\;\;\Gamma$ (12)

It may be noted that this general conducting wall boundary condition
degenerates to the following special cases:

$\displaystyle c_{w}=0:\quad B|_{\Gamma}=0\quad\quad\quad\quad\mbox{(Perfectly
insulating wall)}.$ $\displaystyle c_{w}\rightarrow\infty:\left(\frac{\partial
B}{\partial\xi}\right)_{\Gamma}=0\quad\quad\mbox{(Perfectly conducting
wall)}.$

The boundary condition given in Eq. (11) or (12) is of mixed or Robin type
relating the induced magnetic field with its wall-normal gradient at the
boundary via the wall conductance ratio. Here, we emphasize that the Shercliff
boundary condition avoids the need for computing conjugate MHD solutions
requiring the calculation of the magnetic field within the solid walls $B_{s}$
explicitly. While this can be achieved in principle by computing the magnetic
fields in both the fluid and solid media with the necessary matching of the
conditions at the interface, it often leads to numerically stiff situations
when there is large disparities in the conductivities leading to computational
challenges and inefficiencies with time consuming calculations. Indeed, the
thin wall boundary condition is widely used in studying wall-bounded MHD of
liquid metals in various applications (see e.g., [21, 22, 23, 24, 25, 26, 27,
28, 29, 15]). Hence, in this work, we develop new LB boundary schemes that
implements this boundary condition that effectively accounts for the finite
electrical conductivity of the wall. Also, it would be interesting to develop
LB formulations for solving conjugate MHD problems which may be useful in
certain cases and will form a subject for a future investigation.

##  IV Lattice Kinetic Scheme (LKS) for Magnetic Induction Equation

Before constructing the boundary schemes, we will first briefly discuss the
solution procedure for the MHD equations given in Eq. (1) using the LB
methods. In this regard, we will now look at the LKS introduced by Dellar [2]
for solving the magnetic induction equation (Eq. (1c)) to obtain the magnetic
field. For maintaining simplicity of our presentation without losing
generality, we will only consider the LB formulations in two-dimensions (2D)
and involving the simplest collision model in the following. Recognizing the
inability of scalar distribution functions to represent the antisymmetric
tensor $\bm{u}\bm{B}-\bm{B}\bm{u}$ appearing as a flux term in the left side
of the magnetic induction equation given in Eq. (1c), Dellar [2] proposed the
use of a vector distribution function, whose LB evolution will now be
discussed. The D2Q5 lattice used for this purpose is shown in Fig. 2(a).

Figure 2: (a) D2Q5 lattice for solution of the magnetic induction equation
using a vector distribution function and (b) D2Q9 lattice for the solution of
the Navier-Stokes equations using a scalar distribution function.

The 2D LKS in the index notation reads as

$g_{\alpha i}(\bm{x}+\bm{e}_{\alpha}\delta_{t},t+\delta_{t})-g_{\alpha
i}(\bm{x},t)=-\frac{1}{\tau_{m}}(g_{\alpha i}-g_{\alpha i}^{eq}),$ (13)

where $i\in\\{x,y\\}$ represent the Cartesian spatial directions and
$\alpha=\\{0,1,2,3,4\\}$ denote the discrete particle directions. The right
side of Eq. (13) represents the relaxation of the $i$-th component of the
vector distribution function to its equilibrium $g_{\alpha i}^{eq}$ with a
relaxation time $\tau_{m}$ during a time step $\delta_{t}$. The left side of
Eq. (13) shows the perfect-shift streaming of $g_{\alpha i}$ along the
characteristic directions. Here, the equilibrium distribution function
$g_{\alpha i}^{eq}$ can be written as

$g_{\alpha i}^{eq}=W_{\alpha}\left\\{B_{i}+\frac{e_{\alpha
j}}{c_{sm}^{2}}(u_{j}B_{i}-B_{j}u_{i})\right\\},$ (14)

which is related to the magnetic field vector component $B_{i}$ and the rank-2
antisymmetric tensor $u_{j}B_{i}-B_{j}u_{i}$, and the associated weights
$W_{\alpha}$ are given by

$W_{\alpha}=\begin{cases}1/3&\quad\alpha=0\\\
1/6&\quad\alpha=1,2,3,4\end{cases}$ (15)

Here, typically $c_{sm}^{2}=1/3$. Then, defining the first velocity moment of
the vector distribution function $g_{\alpha i}$ as the following rank-2 tensor
$\Lambda_{ji}$ in the component form as

$\Lambda_{ji}=\sum_{\alpha=0}^{4}e_{\alpha j}g_{\alpha i},$ (16)

whose equilibrium parts are represented by

$\Lambda_{ji}^{(0)}=\sum_{\alpha=0}^{4}e_{\alpha j}g_{\alpha
i}^{eq}=u_{j}B_{i}-B_{j}u_{i},$ (17)

and the (leading) non-equilibrium parts based on the Chapman-Enskog analysis
[53] (discussed in detail later) are given by

$\Lambda_{ji}^{(1)}=\sum_{\alpha=0}^{4}e_{\alpha j}g_{\alpha
i}^{(1)}=\sum_{\alpha=0}^{4}e_{\alpha j}(g_{\alpha i}-g_{\alpha
i}^{eq})=-\tau_{m}c_{sm}^{2}\partial_{j}B_{i}.$ (18)

The solution of the LKS (Eq. (13)), by rewriting it in the form of collide-
and-stream steps, yields the distribution function $g_{\alpha i}$, whose
zeroth and first velocity moments yield the magnetic field $B_{i}$ and the
current density $J_{i}$, respectively, which can be written as

$B_{i}=\sum_{\alpha=0}^{4}g_{\alpha i},$ (19)

and

$J_{i}=\frac{1}{\mu_{m}}\epsilon_{ijk}\partial_{j}B_{k}=-\frac{1}{\mu_{m}}\frac{1}{\tau_{m}c_{sm}^{2}}\epsilon_{ijk}\Lambda_{jk}^{(1)}=-\frac{1}{\mu_{m}}\frac{1}{\tau_{m}c_{sm}^{2}}\epsilon_{ijk}(\Lambda_{jk}-\Lambda_{jk}^{(0)}),$
(20)

where $\epsilon_{ijk}$ is the rank-3 permutation (Levi-Civita) tensor. The
above lattice kinetic algorithm simulates the magnetic induction equation (Eq.
(1c)), where the magnetic resistivity $\eta$ is related to the relaxation time
$\tau_{m}$ via

$\eta=\frac{1}{\sigma\mu_{m}}=c_{sm}^{2}\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}.$
(21)

##  V Lattice Boltzmann Method (LBM) for Fluid Flow

We will now present a LB scheme for solving the flow field part of the MHD
equations (Eq. (1a) and (1b)). Again, for maintaining brevity, we now present
the simplest LBM for simulating the fluid motion based on the single
relaxation time (SRT) model [54], which can readily replaced with any of the
more general or sophisticated collision model when required. The 2D SRT-LBM
for computing the hydrodynamics is based on the scalar distribution function
$f_{\alpha}$ that evolves on the D2Q9 lattice (see Fig. 2(b)) having a greater
degree of isotropy than that used in the previous section to recover the
viscous stress tensor. It can be written as

$f_{\alpha}(\bm{x}+\bm{e}_{\alpha}\delta_{t},t+\delta_{t})-f_{\alpha}(\bm{x},t)=-\frac{1}{\tau}(f_{\alpha}-f_{\alpha}^{eq})+\left(1-\frac{1}{2\tau}\right)S_{\alpha}\delta_{t},$
(22)

where the equilibrium distribution function $f_{\alpha}^{eq}$ and the source
term $S_{\alpha}$, which incorporates the various body forces including the
Lorentz force, can be written, respectively, as

$f_{\alpha}^{eq}=w_{\alpha}\rho\left\\{1+\frac{\bm{e}_{\alpha}\cdot\bm{u}}{c_{s}^{2}}+\frac{(\bm{e}_{\alpha}\cdot\bm{u})^{2}}{2c_{s}^{4}}-\frac{\bm{u}\cdot\bm{u}}{2c_{s}^{2}}\right\\},$
(23)

and

$S_{\alpha}=\frac{(\bm{e}_{\alpha}-\bm{u})}{\rho
c_{s}^{2}}\cdot(\bm{F}_{Lorentz}+\bm{F}_{ext})f_{\alpha}^{eq}.$ (24)

Here, the weights $w_{\alpha}$ are given by

$w_{\alpha}=\begin{cases}4/9&\quad\alpha=0\\\ 1/9&\quad\alpha=1,2,3,4\\\
1/36&\quad\alpha=5,6,7,8\end{cases}$ (25)

In the above, typically, $c_{s}^{2}=1/3$. The solution of the LBM (Eq. (22)),
by rewriting it in the form of collide-and-stream steps (see e.g., [36]),
results in the distribution function $f_{\alpha}$, whose various velocity
moments give the fluid density and the velocity fields, which can be
represented as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha},$ (26) $\displaystyle\rho\bm{u}$
$\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha}\bm{e}_{\alpha}+\frac{1}{2}(\bm{F}_{Lorentz}+\bm{F}_{ext})\delta_{t},$
(27)

and the pressure is given by $p=\rho c_{s}^{2}$. These hydrodynamic fields
satisfy the NS equations (Eq. (1a) and (1b)), where the kinematic viscosity of
the fluid $\nu$ is related to the relaxation time $\tau$ as

$\nu=c_{s}^{2}\left(\tau-\frac{1}{2}\right)\delta_{t}.$ (28)

##  VI Link-based Boundary Scheme for Implementation of Shercliff Boundary
Condition for the Magnetic Field

We will now present a link-based boundary scheme to determine the unknown
incoming distribution function from conducting walls for the LKS given in Sec.
IV so that the boundary condition for the induced magnetic field as presented
in Eq. (12) is enforced. A schematic of the arrangement of the nodes and
lattice links around the boundary is shown in Fig. 3.

Figure 3: Schematic of the arrangement of nodes and lattice links around
conducting wall boundary for the link-based boundary condition scheme.

Here, $\bm{x}_{b}$, $\bm{x}_{f}$, and $\bm{x}_{g}$ refer to the boundary node,
nearest fluid node, and ghost node, respectively, and $\bm{e}_{\alpha}$ and
$\bm{e}_{\bar{\alpha}}$ are the incoming and outgoing particle velocity
directions, where $\bm{e}_{\bar{\alpha}}=-\bm{e}_{\alpha}$. The outward unit
normal and tangent to the wall at the boundary node are represented by
$\hat{n}$ and $\hat{t}$, respectively, and the local wall normal coordinate
pointing into the fluid is denoted by $\xi$. We assume that the wall is
stationary in this section, and that the boundary node $\bm{x}_{b}$, where Eq.
(12) is to be satisfied, is located half-way between the nearest fluid and
ghost nodes, and our goal is to construct the unknown incoming distribution
function $g_{\alpha}(\bm{x}_{f},t+\delta_{t})$ at the nearest fluid node. Let
$B_{i}$ be the induced magnetic field component for which the Shercliff BC is
applicable in the presence of an imposed external magnetic field.

First, using $\partial
B_{i}/\partial\xi=\bm{e}_{\alpha}\cdot\bm{\nabla}B_{i}$, it follows that Eq.
(12) can be rewritten in the following form that is more convenient for
analysis:

$-\bm{e}_{\alpha}\cdot\bm{\nabla}B_{i}+\frac{B_{i}}{c_{w}H}=0\quad\quad\mbox{on}\;\;\bm{x}_{b}.$
(29)

Then, using the Chapman-Enskog (C-E) multiscale expansions [53]

$g_{\alpha i}=g_{\alpha i}^{(0)}+\epsilon g_{\alpha
i}^{(1)}+\cdots,\quad\quad\partial_{t}=\partial_{t_{0}}+\epsilon\partial_{t_{1}}+\cdots,$
(30)

where $\epsilon=\delta_{t}$ is a small perturbation parameter, and applying
the Taylor series expansion to the left side of Eq. (13), we get the
$O(\epsilon)$ evolution equation as

$\partial_{t_{0}}g_{\alpha i}^{(0)}+e_{\alpha j}\partial_{j}g_{\alpha
i}^{(0)}=-\frac{1}{\tau_{m}}g_{\alpha i}^{(1)},$ (31)

from which we obtain the non-equilibrium part of the distribution function
$g_{\alpha i}^{(1)}$, which can be written as

$g_{\alpha i}^{(1)}=-\tau_{m}\partial_{t_{0}}g_{\alpha
i}^{(0)}-\tau_{m}e_{\alpha j}\partial_{j}g_{\alpha i}^{(0)}.$ (32)

Substituting Eq. (32) in Eq. (30), we obtain an estimate of the distribution
function $g_{\alpha i}$ in terms of its equilibrium part $g_{\alpha i}^{(0)}$
and its space and time derivatives as

$g_{\alpha i}=g_{\alpha i}^{(0)}-\tau_{m}\delta_{t}\partial_{t_{0}}g_{\alpha
i}^{(0)}-\tau_{m}\delta_{t}e_{\alpha k}\partial_{k}g_{\alpha i}^{(0)}.$ (33)

Also, using $\bm{e}_{\bar{\alpha}}=-\bm{e}_{\alpha}$, the distribution
function for the opposite direction to that given in Eq. (33) follows as

$g_{\bar{\alpha}i}=g_{\bar{\alpha}i}^{(0)}-\tau_{m}\delta_{t}\partial_{t_{0}}g_{\bar{\alpha}i}^{(0)}+\tau_{m}\delta_{t}e_{\alpha
k}\partial_{k}g_{\bar{\alpha}i}^{(0)}.$ (34)

Now, the equilibrium distribution function given in Eq. (14) in view of Eq.
(17) can be rewritten as a function of the magnetic field $B_{i}$ and the
anti-symmetric tensor $\Lambda_{ji}^{(0)}$ as

$g_{\alpha i}^{(0)}=W_{\alpha}\left\\{B_{i}+\frac{e_{\alpha
j}}{c_{sm}^{2}}\Lambda_{ji}^{(0)}\right\\},$ (35)

and using $W_{\bar{\alpha}}=W_{\alpha}$ and
$\bm{e}_{\bar{\alpha}}=-\bm{e}_{\alpha}$ again, the corresponding expression
in the opposite direction reads as

$g_{\bar{\alpha}i}^{(0)}=W_{\alpha}\left\\{B_{i}-\frac{e_{\alpha
j}}{c_{sm}^{2}}\Lambda_{ji}^{(0)}\right\\}.$ (36)

From the last two equations, we can then write the sum and difference of the
equilibria in the two opposite directions as

$\displaystyle g_{\alpha i}^{(0)}+g_{\bar{\alpha}i}^{(0)}$ $\displaystyle=$
$\displaystyle 2W_{\alpha}B_{i},$ (37a) $\displaystyle g_{\alpha
i}^{(0)}-g_{\bar{\alpha}i}^{(0)}$ $\displaystyle=$ $\displaystyle
2W_{\alpha}\frac{e_{\alpha j}}{c_{sm}^{2}}\Lambda_{ji}^{(0)}.$ (37b)

On the other hand, combining Eqs. (33) and (34), the sum and difference of the
distribution function in the opposite directions follow as

$\displaystyle g_{\alpha i}+g_{\bar{\alpha}i}$ $\displaystyle=$
$\displaystyle(g_{\alpha
i}^{(0)}+g_{\bar{\alpha}i}^{(0)})-\tau_{m}\delta_{t}\left[\partial_{t_{0}}(g_{\alpha
i}^{(0)}+g_{\bar{\alpha}i}^{(0)})+e_{\alpha_{k}}\partial_{k}(g_{\alpha
i}^{(0)}-g_{\bar{\alpha}i}^{(0)})\right],$ (38a) $\displaystyle g_{\alpha
i}-g_{\bar{\alpha}i}$ $\displaystyle=$ $\displaystyle(g_{\alpha
i}^{(0)}-g_{\bar{\alpha}i}^{(0)})-\tau_{m}\delta_{t}\left[\partial_{t_{0}}(g_{\alpha
i}^{(0)}-g_{\bar{\alpha}i}^{(0)})+e_{\alpha_{k}}\partial_{k}(g_{\alpha
i}^{(0)}+g_{\bar{\alpha}i}^{(0)})\right].$ (38b)

Substituting for $(g_{\alpha i}^{(0)}+g_{\bar{\alpha}i}^{(0)})$ and
$(g_{\alpha i}^{(0)}-g_{\bar{\alpha}i}^{(0)})$ from Eq. (37) in Eq. (38), we
finally obtain $g_{\alpha i}\pm g_{\bar{\alpha}i}$ in terms of the macroscopic
field variables in their derivatives as

$\displaystyle g_{\alpha i}+g_{\bar{\alpha}i}$ $\displaystyle=$ $\displaystyle
2W_{\alpha}B_{i}-2W_{\alpha}\tau_{m}\delta_{t}\left[\partial_{t_{0}}B_{i}+\frac{e_{\alpha_{j}}e_{\alpha_{k}}}{c_{sm}^{2}}\partial_{k}\Lambda_{ji}^{(0)}\right],$
(39a) $\displaystyle g_{\alpha i}-g_{\bar{\alpha}i}$ $\displaystyle=$
$\displaystyle
2W_{\alpha}\frac{e_{\alpha_{j}}}{c_{sm}^{2}}\Lambda_{ji}^{(0)}-2W_{\alpha}\tau_{m}\delta_{t}\left[\frac{e_{\alpha_{j}}}{c_{sm}^{2}}\partial_{t_{0}}\Lambda_{ji}^{(0)}+e_{\alpha_{k}}\partial_{k}B_{i}\right].$
(39b)

Next, we evaluate Eq. (39) on the boundary node $\bm{x}_{b}$, where due to the
no-slip BC $u_{i}(\bm{x}_{b})=0$ and from the scaling $O(B_{i})\sim O(u_{i})$,
we neglect terms of $O(u_{i}^{2})$ or smaller. Thus,
$\partial_{t_{0}}B_{i}(\bm{x}_{b})\sim\partial_{k}\Lambda_{ji}^{(0)}\sim
O(u_{i}^{2})$ and $\Lambda_{ji}^{(0)}(\bm{x}_{b})=0$ due to the no-slip BC.
Hence, the expressions in Eq. (39) for $g_{\alpha i}\pm g_{\bar{\alpha}i}$ at
$\bm{x}_{b}$ reduce to

$\displaystyle g_{\alpha i}(\bm{x}_{b})+g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$ $\displaystyle 2W_{\alpha}B_{i}(\bm{x}_{b}),$ (40a)
$\displaystyle g_{\alpha i}(\bm{x}_{b})-g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$
$\displaystyle-2W_{\alpha}\tau_{m}\delta_{t}e_{\alpha_{k}}\partial_{k}B_{i}(\bm{x}_{b}),$
(40b)

which expresses the fact the sum is related to the magnitude of the magnetic
field, while the difference is related to its gradient, with both of them
being evaluated at the boundary node.

Now, the key idea in constructing a link-based boundary condition
implementation for the LKS is to express the boundary node estimates of the
combinations of the distribution functions $g_{\alpha i}(\bm{x}_{b})\pm
g_{\bar{\alpha}i}(\bm{x}_{b})$ given above in terms of the distribution
functions at the nearest fluid node $\bm{x}_{f}$, i.e., $g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})$, which is unknown and the required incoming
quantity from the wall, and $\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$, which
is the outgoing and a known quantity. Here, $\tilde{g}_{\bar{\alpha}i}$ refers
to the post-collision value of $g_{\bar{\alpha}i}$, which follows from Eq.
(13) as

$\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)=g_{\bar{\alpha}i}(\bm{x}_{f},t)-\frac{1}{\tau_{m}}[g_{\bar{\alpha}i}(\bm{x}_{f},t)-g_{\bar{\alpha}i}^{(0)}(\bm{x}_{f},t)].$
(41)

In this regard, we apply the Taylor series in space to the above two
quantities $g_{\alpha i}(\bm{x}_{f},t+\delta_{t})$ and
$\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$ by assuming the boundary node
$\bm{x}_{b}$ is located half-way from the nearest fluid node $\bm{x}_{f}$,
which yield

$\displaystyle g_{\alpha i}(\bm{x}_{f},t+\delta_{t})$ $\displaystyle\approx$
$\displaystyle g_{\alpha
i}(\bm{x}_{b},t+\delta_{t})+\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}g_{\alpha i}(\bm{x}_{b},t+\delta_{t}),$ (42a)
$\displaystyle\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$ $\displaystyle\approx$
$\displaystyle\tilde{g}_{\bar{\alpha}i}(\bm{x}_{b},t)-\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}\tilde{g}_{\bar{\alpha}i}(\bm{x}_{b},t),$ (42b)

where the last equation exploited the fact that
$\bm{e}_{\bar{\alpha}j}=-\bm{e}_{\alpha j}$. Now, on $\bm{x}_{b}$, we take
$g_{\alpha i}(\bm{x}_{b},t+\delta_{t})=g_{\alpha i}(\bm{x}_{b},t)=g_{\alpha
i}(\bm{x}_{b})$ and
$\tilde{g}_{\bar{\alpha}i}(\bm{x}_{b},t)=g_{\bar{\alpha}i}(\bm{x}_{b},t)=g_{\bar{\alpha}i}(\bm{x}_{b})$.
Moreover, in evaluating the spatial gradients $\partial_{j}g_{\alpha i}$ and
$\partial_{j}\tilde{g}_{\bar{\alpha}i}$ on $\bm{x}_{b}$ in the above last two
equations, we consider the leading order equilibrium contributions for the
distribution functions, i.e., $\partial_{j}g_{\alpha
i}(\bm{x}_{b})\approx\partial_{j}g_{\alpha i}^{(0)}(\bm{x}_{b})\approx
W_{\alpha}\partial_{j}B_{i}(\bm{x}_{b})$ and
$\partial_{j}\tilde{g}_{\bar{\alpha}i}(\bm{x}_{b})\approx\partial_{j}g_{\bar{\alpha}i}^{(0)}(\bm{x}_{b})\approx
W_{\alpha}\partial_{j}B_{i}(\bm{x}_{b})$. Thus, we can simplify the Eqs. (42)
to

$\displaystyle g_{\alpha i}(\bm{x}_{f},t+\delta_{t})$ $\displaystyle\approx$
$\displaystyle g_{\alpha
i}(\bm{x}_{b})+W_{\alpha}\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}),$ (43a)
$\displaystyle\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$ $\displaystyle\approx$
$\displaystyle\tilde{g}_{\bar{\alpha}i}(\bm{x}_{b})-W_{\alpha}\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}).$ (43b)

These last two equations in Eq. (43) can be rearranged further as

$\displaystyle g_{\alpha i}(\bm{x}_{b})$ $\displaystyle\approx$ $\displaystyle
g_{\alpha i}(\bm{x}_{f},t+\delta_{t})-W_{\alpha}\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}),$ (44a) $\displaystyle
g_{\bar{\alpha}i}(\bm{x}_{b})$ $\displaystyle\approx$
$\displaystyle\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)+W_{\alpha}\frac{\delta_{t}}{2}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}).$ (44b)

Now, combining these two equations in Eq. (44) as their sum and difference, we
can then express them in the following equivalent forms:

$\displaystyle g_{\alpha i}(\bm{x}_{b})+g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$ $\displaystyle g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t),$ (45a)
$\displaystyle g_{\alpha i}(\bm{x}_{b})-g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$ $\displaystyle g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)-W_{\alpha}\delta_{t}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}).$ (45b)

Then, using Eq. (45) in Eq. (40), we eliminate the intermediate quantities
$g_{\alpha i}(\bm{x}_{b})\pm g_{\bar{\alpha}i}(\bm{x}_{b})$ and relate the
unknown incoming distribution function $g_{\alpha i}(\bm{x}_{f},t+\delta_{t})$
to the known outgoing distribution function
$\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$ and the boundary node values of the
magnetic field $B_{i}(\bm{x}_{b})$ and its gradient
$\partial_{j}B_{i}(\bm{x}_{b})$. As results of these steps, we get

$\displaystyle\frac{1}{2}\left[g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)\right]$
$\displaystyle=$ $\displaystyle W_{\alpha}B_{i}(\bm{x}_{b}),$ (46a)
$\displaystyle-\frac{1}{2\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}}\left[g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)\right]$
$\displaystyle=$ $\displaystyle W_{\alpha}\delta_{t}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b}).$ (46b)

The relations in Eq. (46) together need to satisfy the constraint given by the
Shercliff boundary condition for thin conducting walls, which can be rewritten
from Eq. (29) by multiplying it with the weighting factor $W_{\alpha}$ as

$W_{\alpha}B_{i}(\bm{x}_{b})-c_{w}HW_{\alpha}e_{\alpha
j}\partial_{j}B_{i}(\bm{x}_{b})=0.$ (47)

Substituting the expressions in Eq. (46) into Eq. (47), we get

$\left[1+\frac{c_{w}H}{\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}}\right]g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})+\left[1-\frac{c_{w}H}{\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}}\right]\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)=0.$
(48)

Based on Eq. (48), we then introduce a weighting factor $\theta$ related to
the wall conductance ratio $c_{w}$ and the relaxation time $\tau_{m}$ as

$\theta=\frac{c_{w}H}{\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}}.$ (49)

Then, Eq. (48) can be rewritten as

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=\frac{1}{1+\theta}[-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)]+\frac{\theta}{1+\theta}[+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)],\quad\quad
0\leq\theta\leq\infty.$ (50)

Equation (50) represents the new link-based BC closure which reconstructs the
incoming distribution function for the induced magnetic field component from
thin conducting walls via a weighted combination of the anti-bounce back
(first term) and bounce back (second term) of the outgoing post-collision
distribution function. It follows that the above boundary scheme reduces to
the following limiting cases:

(a) Insulated wall: $c_{w}=\theta=0$, which corresponds to

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)\quad\mbox{(Anti-
bounce back scheme)},$

(b) Perfectly conducting wall: $c_{w}=\theta\rightarrow\infty$, which
corresponds to

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)\quad\mbox{(Bounce
back scheme)},$

On the other hand, the implementation of the link-based no-slip velocity
boundary condition via the standard half-way bounce back scheme [37] implies
that
$f_{\alpha}(\bm{x}_{f},t+\delta_{t})=\tilde{f}_{\bar{\alpha}}(\bm{x}_{f},t)$.
A derivation of the extension of the above boundary formulation given for a
stationary conducting wall to a moving conducting wall is presented in B.1.

###  VI.1 Example implementation: Hartmann flow through a channel bounded by
conducting walls

We will now illustrate the implementation of the link-based boundary scheme
for the body force-driven flow between two conducting plates in the presence
of an external magnetic field. A schematic of this flow set up is shown in
Fig. 4, where an electrically conducting fluid with conductivity $\sigma$
bounded by plates, each of which has a thickness $\delta$ and electrical
conductivity $\sigma_{s}$ and a separated by a distance $2a$. The fluid is in
motion due to a body force $F_{x}$ applied in the $x$ direction in the
presence of an imposed magnetic field $B_{0}$ normal to the bounding walls,
i.e., in the $y$ direction.

Figure 4: Sketch of the Hartmann flow bounded by two parallel plates with
finite electrical conductivities.

The analytical solution for this generalized Hartmann flow problem with finite
wall conductivity effects was given by Chang and Lundgren [21], which will be
summarized later in Sec. VIII when discussing the numerical results. For this
case, as shown in Fig. 4, at each boundary, the Shercliff thin wall boundary
condition is applicable for the induced magnetic in the $x$ direction, i.e.,
$B_{x}$, while the magnetic field in the $y$ direction $B_{y}$ is based on the
imposed magnetic field $B_{0}$. Here, the wall conductance ratio $c_{w}$ and
the weighting factor $\theta$ can be written as

$c_{w}=\frac{\sigma_{s}\delta}{\sigma
a},\quad\quad\theta=\frac{c_{w}a}{(\tau_{m}-\frac{1}{2})\delta_{t}}.$

Referring to Fig. 2(a) for the lattice directions, the link-based boundary
closure relations for the incoming distribution functions corresponding to the
components of the magnetic field $B_{x}$ and $B_{y}$ for the bottom and top
walls follow from the above considerations read as

$\displaystyle y=-a:\quad g_{2x}$ $\displaystyle=$
$\displaystyle\frac{\theta-1}{\theta+1}\tilde{g}_{4x},$ $\displaystyle g_{2y}$
$\displaystyle=$ $\displaystyle-\tilde{g}_{4y}+2W_{2}B_{0},$ $\displaystyle
y=+a:\quad g_{4x}$ $\displaystyle=$
$\displaystyle\frac{\theta-1}{\theta+1}\tilde{g}_{2x},$ $\displaystyle g_{4y}$
$\displaystyle=$ $\displaystyle-\tilde{g}_{4y}+2W_{4}B_{0}.$

That is, the components of $g_{2x}$ and $g_{4x}$ for the bottom and top walls,
respectively, follow from Eq. (50), while that for the wall normal components
$g_{2y}$ and $g_{4y}$ that satisfy the Dirichlet boundary condition based on
$B_{0}$ directly result from Eq. (46a).

##  VII On-Node Moment-Based Boundary Scheme for Implementation of Shercliff
Boundary Condition for the Magnetic Field

We will now present an alternative moment-based approach, which implements the
boundary conditions directly by means of constraints on the moments of the
distribution functions on the nodes rather than on some cut-link locations.
Originally proposed by Noble [39] for the D2Q7 lattice and then applied for
the D2Q9 lattice by Bennett [40] for specifying the no-slip velocity boundary
conditions, it was subsequently extended for the distribution functions for
the magnetic fields for representing the specific case of the insulated wall
MHD by Dellar [43]. Here, we will consider a more general situation involving
MHD with finite wall conductivities, which require the specification of the
mixed or Robin-type boundary condition.

Consider again the Hartmann flow case bounded by thin conducting walls
discussed at the end of the previous section (see Fig. 4). Our goal is to
impose the following boundary exactly via the necessary moment constraints:

$\displaystyle-\frac{\partial B_{x}}{\partial y}+\frac{B_{x}}{c_{w}a}$
$\displaystyle=$ $\displaystyle 0,\quad
B_{y}=B_{0}\quad\quad\mbox{at}\;\;y=-a,$ (51a) $\displaystyle+\frac{\partial
B_{x}}{\partial y}+\frac{B_{x}}{c_{w}a}$ $\displaystyle=$ $\displaystyle
0,\quad B_{y}=B_{0}\quad\quad\mbox{at}\;\;y=+a.$ (51b)

Unlike the previous approach, here $\bm{x}_{b}$ is located on directly on the
lattice node and not on a link, i.e., we consider the boundary to be ‘wet’.
Decomposing the tensor component $\Lambda_{yx}$ by means of its equilibrium
and non-equilibrium parts, i.e.,
$\Lambda_{yx}=\Lambda_{yx}^{(0)}+\delta_{t}\Lambda_{yx}^{(1)}$, and evaluating
this at the wall boundary where
$\Lambda_{yx}^{(0)}=(u_{y}B_{x}-B_{y}u_{x})(\bm{x}_{b})=0$ via the no-slip
condition, we get $\Lambda_{yx}=\delta_{t}\Lambda_{yx}^{(1)}$. From Eq. (18),
which follows from a C-E analysis as given in the derivation in Sec. VI, we
have $\Lambda_{yx}^{(1)}=-\tau_{m}c_{sm}^{2}\partial_{y}B_{x}$. Thus, we can
express the boundary node evaluation of the tensor component $\Lambda_{yx}$ as

$\Lambda_{yx}(\bm{x}_{b})=-\tau_{m}c_{sm}^{2}\delta_{t}\frac{\partial
B_{x}}{\partial y}(\bm{x}_{b}),$ (52)

which implies that the wall normal gradient of the induced magnetic field at
the boundary can be written in terms of the first moment involving the
velocity in the $y$ direction of the distribution function in the $x$
direction $g_{\alpha x}$ at $\bm{x}_{b}$ as follows:

$\frac{\partial B_{x}}{\partial
y}(\bm{x}_{b})=-\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}\Lambda_{yx}(\bm{x}_{b})=-\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}\sum_{\alpha=0}^{4}e_{\alpha
y}g_{\alpha x}(\bm{x}_{b}),$ (53)

On the other hand, the zeroth moment of $g_{\alpha x}$ at the boundary node
gives the magnetic field component $B_{x}$ as

$B_{x}(\bm{x}_{b})=\sum_{\alpha=0}^{4}g_{\alpha x}(\bm{x}_{b}).$ (54)

Then, using Eqs. (53) and (54) in the Shercliff boundary condition for the
bottom wall in Eq. (51a), we get

$\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}\sum_{\alpha=0}^{4}e_{\alpha
y}g_{\alpha x}(\bm{x}_{b})+\frac{1}{c_{w}a}\sum_{\alpha=0}^{4}g_{\alpha
x}(\bm{x}_{b})=0,$

which simplifies to

$\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}(g_{2x}-g_{4x})+\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{2x}+g_{3x}+g_{4x})=0\quad\quad\mbox{at}\;\;y=-a.$
(55)

On the other hand $B_{y}=B_{0}$ at $\bm{x}_{b}$ (see Eq. (51a)) implies

$g_{0y}+g_{1y}+g_{2y}+g_{3y}+g_{4y}=B_{0}\quad\quad\mbox{at}\;\;y=-a.$ (56)

Then, solving for the unknowns $g_{2x}$ and $g_{2y}$ from Eqs. (55) and (56),
respectively, we finally get the moment-based closures for the distribution
functions at the bottom wall as

$\displaystyle g_{2x}$ $\displaystyle=$
$\displaystyle\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}+\frac{1}{c_{w}a}\right]^{-1}\left\\{\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}-\frac{1}{c_{w}a}\right]g_{4x}-\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{3x})\right\\}\quad\quad\mbox{at}\;\;y=-a,$
(57a) $\displaystyle g_{2y}$ $\displaystyle=$ $\displaystyle
B_{0}-(g_{0y}+g_{1y}+g_{3y}+g_{4y}).$ (57b)

Similarly, at $y=+a$, from Eq. (51b) and using Eqs. (53) and (54), we get

$-\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}(g_{2x}-g_{4x})+\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{2x}+g_{3x}+g_{4x})=0\quad\quad\mbox{at}\;\;y=+a,$
(58)

and $B_{y}=B_{0}$ implies that

$g_{0y}+g_{1y}+g_{2y}+g_{3y}+g_{4y}=B_{0}\quad\quad\mbox{at}\;\;y=+a.$ (59)

Rearranging Eqs. (58) and (59) finally yield the incoming distribution
functions from the top wall for the moment-based formulation as

$\displaystyle g_{4x}$ $\displaystyle=$
$\displaystyle\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}+\frac{1}{c_{w}a}\right]^{-1}\left\\{\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}-\frac{1}{c_{w}a}\right]g_{2x}-\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{3x})\right\\}\quad\quad\mbox{at}\;\;y=+a,$
(60a) $\displaystyle g_{4y}$ $\displaystyle=$ $\displaystyle
B_{0}-(g_{0y}+g_{1y}+g_{2y}+g_{3y}).$ (60b)

As in the previous section, we now analyze the limiting cases involving the
wall conductance ratio $c_{w}$, and in what follows, and will express below
just the $x$ component of the distribution function which are dependent in
this quantity. The $y$ component expressions given above are based on the
imposed magnetic field $B_{0}$ and are invariant under changes in $c_{w}$. For
the lower boundary at $y=-a$, from Eq. (57a), we get

$\displaystyle c_{w}=0:\quad g_{2x}$ $\displaystyle=$
$\displaystyle-(g_{0x}+g_{1x}+g_{3x}+g_{4x})\quad\quad\quad\quad\mbox{(Perfectly
insulating wall)},$ $\displaystyle c_{w}\rightarrow\infty:g_{2x}$
$\displaystyle=$ $\displaystyle
g_{4x}\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\;\mbox{(Perfectly conducting
wall)}.$

Similarly, degenerate cases at the upper boundary ($y=+a$) follow from Eq.
(60a) as

$\displaystyle c_{w}=0:\quad g_{4x}$ $\displaystyle=$
$\displaystyle-(g_{0x}+g_{1x}+g_{2x}+g_{3x})\quad\quad\quad\quad\mbox{(Perfectly
insulating wall)},$ $\displaystyle c_{w}\rightarrow\infty:g_{4x}$
$\displaystyle=$ $\displaystyle
g_{2x}\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\;\mbox{(Perfectly conducting
wall)}.$

We note here that the results given above for the insulated wall condition
match with those presented in Ref. [43]. For completeness, the moment-based
implementation of the velocity boundary conditions needed for the coupled MHD
solution is briefly summarized in A. Moreover, B.2 discusses the extension of
the above moment-based boundary scheme for MHD with moving conducting walls.
It may be noted that the link-based boundary scheme as well as the moment-
based boundary formulation derived in the last two sections can be implemented
to resolve the Hartmann layers as well as the Shercliff or side layers in MHD
bounded by conducting walls, which arise in many practical problems especially
under high magnetic field strengths, via using appropriate local grid
refinement approaches.

##  VIII Results and Discussion

###  VIII.1 MHD bounded by conducting walls with stationary boundaries: Body
force-driven flow

We will first consider a case study involving the Hartmann flow through a
channel bounded by thin conducting walls to perform a numerical investigations
of the two new boundary closures for the LKS for representing the Shercliff BC
that were discussed earlier. A schematic of this flow problem is shown in Fig.
4 and its description is given in Sec. VI.1. For this liquid metal fluid of
kinematic viscosity $\nu$, density $\rho$, electrical conductivity $\sigma$
flowing due a body force of magnitude $F_{x}$ in the presence of an external
magnetic field $B_{0}$ and bounded by plates with conductance wall ratio
$c_{w}$ and separated by a distance $2a$, the analytical solutions for the
profiles of the velocity and the induced magnetic fields were derived by Chang
and Lundgren [21], and are given by

$\displaystyle u^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{u(y)}{U_{ref}}=\frac{1}{\mbox{Ha}}\left[\frac{c_{w}+1}{\mbox{Ha}\;c_{w}+\tanh(\mbox{Ha})}\right]\left[1-\frac{\cosh(\mbox{Ha}\xi)}{\cosh(\mbox{Ha})}\right],$
$\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{B_{x}(y)}{B_{ref}}=-\frac{\xi}{\mbox{Ha}}+\frac{1}{\mbox{Ha}}\left[\frac{c_{w}+1}{\mbox{Ha}\;c_{w}+\tanh(\mbox{Ha})}\right]\frac{\sinh(\mbox{Ha}\xi)}{\cosh(\mbox{Ha})},$

where $\xi=y/a$ and $-1\leq\xi\leq 1$. See also Ref. [10]. Here, the reference
velocity $U_{ref}$, the reference scale for the magnetic field $B_{ref}$, as
well as the Hartmann number Ha are given by

$U_{ref}=\frac{a^{2}}{\rho\nu}F_{x},\quad\quad
B_{ref}=a^{2}F_{x}\mu_{m}\left(\frac{\sigma}{\rho\nu}\right)^{1/2},\quad\quad\mbox{Ha}=B_{0}a\left(\frac{\sigma}{\rho\nu}\right)^{1/2},$

where the latter represents the ratio of the electromagnetic or the Lorentz
force to the viscous force. We consider periodic boundary condition along the
flow, i.e., $x$ direction, and the no-slip velocity boundary conditions as
well as the conducting wall condition for the induced magnetic field in the
$x$ direction and the external magnetic field in the $y$ direction at the
bottom and top walls are imposed using the link-based and the moment-based
boundary schemes discussed in the previous sections. For validating these two
boundary schemes, we considered two different values of the Harmann number
$\mbox{Ha}=5$ and $\mbox{Ha}=20$, and in each case, the wall conductance ratio
$c_{w}$ is specified to have a range of values given by
$c_{w}=\\{0,0.1,0.5,1.0,10^{4}\\}$. In this regard, the computational domain
is resolved using $3\times 101$ grid nodes, and the relaxation times of the LB
schemes used to compute the velocity field and magnetic field, respectively,
are chosen as $\tau=0.57$ and $\tau_{m}=0.75$. In each of the chosen set of
parameters, the simulations are run until steady state.

Figures 5 and 6 present comparisons between the computed profiles of the
magnetic field obtained using the link-based and moment-based boundary
schemes, respectively, with the analytical solution at the two Hartmann
numbers and for the range of the wall conductance ratio given above. For
clarity, the profiles along only the top half of the channel are shown and it
may be noted that in the bottom half of the channel, the magnetic field
profiles are antisymmetric to that in its top half. When $c_{w}$ is zero,
corresponding to the insulated wall case, as expected, the induced magnetic
field is seen to zero. As $c_{w}$ increases, the finite wall conductivity
effects become more prominent, with non-zero values of the induced magnetic
field on the wall and accompanied by progressively decreasing magnitudes of
its slope normal to the wall are observed. While both the boundary condition
implementations are found to work for any values of the wall conductance ratio
(i.e., $0\leq c_{w}\leq\infty$), it was noticed that when $c_{w}=10^{4}$, the
wall becomes essentially perfectly conducting with zero wall-normal gradient
for the induced magnetic field. For the smaller Ha case, variations in the
magnetic fields are seen to occur for the most part of the channel; by
contrast, for the higher Ha case, such variations are confined within the
relatively thinner Hartmann layer. It is evident that there is an excellent
agreement with the analytical solution for the structures of the magnetic
field for both the boundary schemes for all the range of values considered for
the electrical properties of the walls for both the Hartmann numbers.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 5: The magnetic field profiles $B_{x}(y)$ across the top half of the
channel for the Hartmann flow bounded by conducting walls for different values
of the wall conductance ratio $c_{w}=0,0.1,0.5,1$, and $10^{4}$ at Hartmann
numbers (a) $\mbox{Ha}=5$ and (b) $\mbox{Ha}=20$ computed via the link-based
boundary scheme (lines) compared with the analytical solution of Chang and
Lundren (1961) [21] (symbols).

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 6: The magnetic field profiles $B_{x}(y)$ across the top half of the
channel for the Hartmann flow bounded by conducting walls for different values
of the wall conductance ratio $c_{w}=0,0.1,0.5,1$, and $10^{4}$ at Hartmann
numbers (a) $\mbox{Ha}=5$ and (b) $\mbox{Ha}=20$ computed via the moment-based
boundary scheme (lines) compared with the analytical solution of Chang and
Lundren (1961) [21] (symbols).

Moreover, in Figs. 7 and 8 the respective comparisons for the velocity field
obtained using link-based and moment-based boundary implementations are shown.
Since the Lorentz force couples such strong variations in the induced magnetic
fields at different values of the wall conductance ratio, it is expected that
the latter will also have significant effect on the velocity field across the
channel. In general, the higher the wall conductance ratio, the smaller are
the magnitudes of the velocity field for any given Ha. This is consistent with
the fact that walls with finite conductivities ($c_{w}\neq 0$) supports a
tangential flow of current within them, whose interaction with the external
magnetic field generates larger Lorentz force that resists the fluid motion
when compared to the insulated wall case ($c_{w}=0$), where no current flow
exists within the wall. Moreover, as Ha is increased, as the Hartmann layer
becomes thinner, there is a significant flattening of the velocity profiles
within the core of the flow. Again, the analytical predictions of these
variations in the velocity field profiles at different $c_{w}$ and Ha are
found to be very well reproduced by both the link-based and moment-based
boundary schemes constructed in the previous sections.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 7: The velocity field profiles $u(y)$ across the top half of the
channel for the Hartmann flow bounded by conducting walls for different values
of the wall conductance ratio $c_{w}=0,0.1,0.5,1$, and $10^{4}$ at Hartmann
numbers (a) $\mbox{Ha}=5$ and (b) $\mbox{Ha}=20$ computed via the link-based
boundary scheme (lines) compared with the analytical solution of Chang and
Lundren (1961) [21] (symbols).

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 8: The velocity field profiles $u(y)$ across the top half of the
channel for the Hartmann flow bounded by conducting walls for different values
of the wall conductance ratio $c_{w}=0,0.1,0.5,1$, and $10^{4}$ at Hartmann
numbers (a) $\mbox{Ha}=5$ and (b) $\mbox{Ha}=20$ computed via the moment-based
boundary scheme (lines) compared with the analytical solution of Chang and
Lundren (1961) [21] (symbols).

Next, let us estimate the order of accuracy under grid refinement for both the
boundary implementation schemes. In this regard, we define the following error
estimates between the computed and analytical solutions of the velocity and
magnetic fields under the second norm:

$\displaystyle\mbox{Error}_{u}$ $\displaystyle=$
$\displaystyle||u_{c}-u_{a}||_{2}=\left[\frac{\sum_{i}(u_{c,i}-u_{a,i})^{2}}{N}\right]^{1/2},$
$\displaystyle\mbox{Error}_{B_{x}}$ $\displaystyle=$
$\displaystyle||B_{xc}-B_{xa}||_{2}=\left[\frac{\sum_{i}(B_{xc,i}-B_{xa,i})^{2}}{N}\right]^{1/2},$

where $u_{c}$ and $u_{a}$ represent the computed and analytical solutions of
the velocity field, while $B_{xc}$ and $B_{xa}$ denote the computed and
analytical results for the magnetic field. The subscript $i$ and $\sum_{i}$
denote the summation is carried out across the width of the channel resolved
by $N$ number of grid nodes. For consistent convergence to the incompressible
limit, we employ the usual diffusive scaling, where for fixed relaxation times
$\tau$ and $\tau_{m}$, the scales of the velocity and magnetic fields are
varied in proportion to the spatial discretization length scales given by
$1/N$. In other words, the compressibility errors decrease at the same rate as
the spatial truncation errors in this limit process to represent the
incompressible MHD solutions. In particular, for each of the boundary scheme,
we consider $\mbox{Ha}=20$ and fix the relaxation times $\tau=0.9$ and
$\tau_{m}=0.95$, and vary the grid resolution by choosing $N=50,100$ and $200$
for the following different values of the wall conductance ratio:
$c_{w}=0,0.1,0.5$, and $1.0$. The errors between the computed solutions
obtained using the link-based boundary scheme and analytical solutions for the
magnetic field and the velocity field at different values of the wall
conductance ratios for various choices of the grid nodes are plotted in Figs.
9(a) and 9(b), respectively. Similar errors versus the number of grid nodes
for the magnetic field and velocity field in the case of the moment-based
boundary scheme are shown in Figs. 10(a) and 10(b), respectively. In both
these figures, a reference slope of $-2$ is also displayed. Clearly, both the
link-based and moment-based boundary schemes for LB methods for MHD show
second order accuracy in computing solutions for the magnetic and velocity
fields for the Hartmann flow with conducting walls. These are also consistent
with the second order accurate results arising from an error analysis made for
the special case involving the insulated walls in Ref. [43].

(a)

(b)

Figure 9: Grid convergence of the (a) magnetic field $B_{x}$ and (b) velocity
field $u$ using the LB method for MHD using the link-based boundary scheme as
a function of the number of grid nodes in the wall normal direction $N$ at
Hartmann number $\mbox{Ha}=20$ for different values of the wall ratio
($c_{w}=0,0.1,0.5,1$) with $\tau=0.9$ and $\tau_{m}=0.95$.

(a)

(b)

Figure 10: Grid convergence of the (a) magnetic field $B_{x}$ and (b) velocity
field $u$ using the LB method for MHD using the moment-based boundary scheme
as a function of the number of grid nodes in the wall normal direction $N$ at
Hartmann number $\mbox{Ha}=20$ for different values of the wall ratio
($c_{w}=0,0.1,0.5,1$) with $\tau=0.9$ and $\tau_{m}=0.95$.

###  VIII.2 MHD bounded by conducting walls with a moving boundary: Shear-
driven flow

We will now present a case study that provides a validation of the extension
of the boundary implementation schemes for moving walls given in B, where B.1
focuses on the link-based formulation and B.2 on the moment-based scheme. In
this regard, we consider the Hartmann-Couette flow of an electrically
conducting fluid with conductivity $\sigma$ enclosed between lower and upper
plates having the wall conductance ratios of $c_{wl}$ and $c_{wu}$,
respectively, with a spacing of $2a$ between them. The shear flow is set up by
the motion of the top plate at a constant velocity $U_{w}$ in the $x$
direction, with the bottom plate being at rest, and it is subjected to an
imposed magnetic field $B_{0}$ in the $y$ direction. For this general MHD flow
configuration involving contrasts in the electrical properties of the two
bounding plates, we have presented a new derivation of the analytical solution
in C, which will be utilized as a benchmark reference against which the
numerical results obtained using our two boundary implementation approaches
will be compared. The simulations are set up with a grid resolution of
$3\times 101$ grid nodes, where the top plate velocity is chosen to be
$U_{w}=0.05$ and using the relaxation parameters $\tau=0.6$ and
$\tau_{m}=0.65$ for the LB schemes for solution of the velocity field and
magnetic field, respectively. Two different values of the Hartmann number
$\mbox{Ha}=5$ and $20$, based on which the imposed magnetic field $B_{0}$ is
determined, and different combinations of the wall conductance ratios
$(c_{wu},c_{wl})$ for the bounding plates are considered, and the computations
in each case are run until the results converge to a steady state.

Figure 11 presents a comparison between the computed magnetic field profiles
obtained using the link-based boundary scheme and the analytical solution for
two values of the Hartmann number ($\mbox{Ha}=5$ and $\mbox{Ha}=20$) and for
the following combinations of the values of the upper and lower wall
conductance ratios: $(c_{wu},c_{wl})=(0,0)$, $(c_{wu},c_{wl})=(0.5,0.5)$, and
$(c_{wu},c_{wl})=(0,0.5)$. Clearly, computed solutions are in excellent
agreement with the analytical solutions for all the choices of wall
conductance ratios and the Hartmann numbers. In general, the strengths of the
induced magnetic fields are seen to increase for cases where one or both the
walls are conducting when compared to the case where both the walls are
electrically insulated. Moreover, as the Hartmann number is increased for a
given pair of the conductance ratios for the top and bottom walls, the induced
magnetic field becomes larger in magnitude, and a flattening effect of their
profiles is observed and accompanied by the wall boundary layer becoming
thinner.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 11: The magnetic field profiles $B_{x}(y)$ across the channel for the
shear-driven Hartmann-Couette flow bounded by upper (u) and lower (l)
conducting walls for different pairs of values of the wall conductance ratios
$(c_{wu},c_{wl})=(0,0)$, $(c_{wu},c_{wl})=(0.5,0.5)$, and
$(c_{wu},c_{wl})=(0,0.5)$ at Hartmann numbers (a) $\mbox{Ha}=5$ and (b)
$\mbox{Ha}=20$ computed via the link-based boundary scheme (lines) compared
with the analytical solution (symbols) given in Eqs. (91) and (92) in C.

The comparisons made for the corresponding velocity profiles for the same
choices of $(c_{wu},c_{wl})$ and Ha in Fig. 12 show that the numerical
solutions based on the link-based scheme for a moving wall are again in very
good agreement with the reference exact solution derived in B.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 12: The velocity field profiles $u(y)$ across the channel for the
shear-driven Hartmann-Couette flow bounded by upper (u) and lower (l)
conducting walls for different pairs of values of the wall conductance ratios
$(c_{wu},c_{wl})=(0,0)$, $(c_{wu},c_{wl})=(0.5,0.5)$, and
$(c_{wu},c_{wl})=(0,0.5)$ at Hartmann numbers (a) $\mbox{Ha}=5$ and (b)
$\mbox{Ha}=20$ computed via the link-based boundary scheme (lines) compared
with the analytical solution (symbols) given in Eqs. (91) and (92) in C.

Notice that imposing an external magnetic field in a shear flow results in
velocity profiles with a pair of inflection points away from each wall, which
can be contrasted with the linear velocity profile arising in the absence of
an imposed magnetic field. The layer formed between the wall and such an
inflection point is seen to become narrower as Ha increases. Moreover, unlike
the case for the induced magnetic field profiles in Fig. 11, the velocity
profiles shown in Fig. 12 are invariant if the values of the conductance ratio
of the top and bottom walls are the same, and variations in them are observed
only when there is a contrast in the conductance ratios at the top and bottom
walls. Thus, the velocity profiles for the cases $(c_{wu},c_{wl})=(0,0)$ and
$(c_{wu},c_{wl})=(0.5,0.5)$ are identical, which differ from that arising from
the choice $(c_{wu},c_{wl})=(0,0.5)$ that involves a contrast in the
electrical properties of the bounding walls. This interesting finding could be
useful in modulating the structure of shear-driven MHD flows based on
manipulating the contrasts in the conductance ratios of the confining walls.

Simulations performed using the same set of parameters as mentioned above
using the moment-based scheme for moving walls given in B.2 demonstrate that
the resulting induced magnetic field and the velocity field compare very well
with the analytical solutions as shown in Figs. 13 and 14.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 13: The magnetic field profiles $B_{x}(y)$ across the channel for the
shear-driven Hartmann-Couette flow bounded by upper (u) and lower (l)
conducting walls for different pairs of values of the wall conductance ratios
$(c_{wu},c_{wl})=(0,0)$, $(c_{wu},c_{wl})=(0.5,0.5)$, and
$(c_{wu},c_{wl})=(0,0.5)$ at Hartmann numbers (a) $\mbox{Ha}=5$ and (b)
$\mbox{Ha}=20$ computed via the moment-based boundary scheme (lines) compared
with the analytical solution (symbols) given in Eqs. (91) and (92) in C.

(a) $\mbox{Ha}=5$

(b) $\mbox{Ha}=20$

Figure 14: The velocity field profiles $u(y)$ across the channel for the
shear-driven Hartmann-Couette flow bounded by upper (u) and lower (l)
conducting walls for different pairs of values of the wall conductance ratios
$(c_{wu},c_{wl})=(0,0)$, $(c_{wu},c_{wl})=(0.5,0.5)$, and
$(c_{wu},c_{wl})=(0,0.5)$ at Hartmann numbers (a) $\mbox{Ha}=5$ and (b)
$\mbox{Ha}=20$ computed via the moment-based boundary scheme (lines) compared
with the analytical solution (symbols) given in Eqs. (91) and (92) in C.

Figure 15 demonstrates the ability of the moment-based boundary scheme in
accurately representing the analytical predictions for a case with fixing
$\mbox{Ha}=5$ and maintaining one of the walls insulated (i.e., $c_{wl}=0$),
and varying the wall conductance ratio of the other wall across a wide range
of values (i.e., $c_{wu}=0,0.1,0.5,1.0$ and $10^{4}$). Notice that there are
dramatic variations in the structures of both the induced magnetic fields and
velocity fields as the top wall’s electrical property varies from being
insulated ($c_{wu}=0$) to its close to being a perfect conductor
($c_{wu}=10^{4}$). In particular, there is a qualitative change in the
velocity profile as $c_{wu}$ is increased with the inflection points shifting
towards the moving top wall and ultimately vanishing when it is perfectly
conducting where a monotonic variation is observed. For completeness, shown in
Fig. 16 are comparisons between the results obtained using the link-based
boundary scheme and the analytical solution for the induced magnetic field and
velocity profiles for this MHD shear flow problem using the same conditions
given below, which confirms excellent agreement between them. Clearly, the
complex interactions among the magnetic and velocity fields, as well as the
electrical properties of the walls in such MHD flows are well reproduced by
our new boundary schemes for LB approaches.

(a)

(b)

Figure 15: The (a) magnetic field $B_{x}(y)$ and (b) velocity field $u(y)$
across the channel for the shear-driven Hartmann-Couette flow bounded by upper
(u) and lower (l) conducting walls at a fixed $\mbox{Ha}=5$ and lower wall
conductance ratio $c_{wl}=0.0$ (insulated) with varying values of the upper
wall conductance ratio $c_{wu}=0,0.1,0.5,1.0$ and $10^{4}$ computed via the
moment-based boundary scheme (lines) compared with the analytical solution
(symbols) given in Eqs. (91) and (92) in C.

(a)

(b)

Figure 16: The (a) magnetic field $B_{x}(y)$ and (b) velocity field $u(y)$
across the channel for the shear-driven Hartmann-Couette flow bounded by upper
(u) and lower (l) conducting walls at a fixed $\mbox{Ha}=5$ and lower wall
conductance ratio $c_{wl}=0.0$ (insulated) with varying values of the upper
wall conductance ratio $c_{wu}=0,0.1,0.5,1.0$ and $10^{4}$ computed via the
link-based boundary scheme (lines) compared with the analytical solution
(symbols) given in Eqs. (91) and (92) in C.

##  IX Summary, Conclusions and Outlook

The electrical properties of the bounding solid wall surfaces in MHD play an
important role in determining the structure of the induced magnetic field and
the velocity field. The finite electrical conductivities of the walls can be
effectively represented using the Shercliff thin wall boundary condition,
which is parameterized by a wall conductance ratio. In this paper, we have
presented two new boundary implementation schemes for such a conducting wall
condition for the LB algorithm for the magnetic field induction. These include
a link-based approach and an on-node moment-based method which provide the
necessary closure relations for the incoming distribution functions for both
stationary and moving conducting walls. Numerical simulations involving body
force-driven as well as shear-driven MHD flows validate the accuracy of both
the boundary schemes for a wide range of the wall conductance ratio. For the
latter shear-driven MHD flow case, we have derived new analytical solutions
under a more general situation involving different conductance ratios for
different bounding wall surfaces, which exhibit interesting flow physics
dependent on the contrasts in such properties. The numerical predications
based on the LB boundary schemes are also shown to be in good agreement with
exact solutions for such more general cases. In addition, a second order grid
convergence of the induced magnetic field as well as the velocity field using
both the proposed boundary formulations is demonstrated.

We note here that the moment-based boundary scheme preserves the necessary
constraints the moments at the boundary, which coincides with the lattice
nodes, exactly, but are limited to straight boundaries. On the other hand, the
present link-based boundary formulation is applicable for curved boundaries
provided the boundary location is considered half-way between the lattice
nodes. Both the boundary schemes are readily applicable in three-dimensions
and with other more general collision models.

Several extensions of this work are envisaged, which are subjects for future
investigations. These include extending the link-based approach for curved
boundaries with general cut-link distances via suitable interpolations for
performing LB simulations of MHD with conducting walls in complex geometries.
Moreover, it may be beneficial to formulate link-based boundary schemes
involving multiple relaxation times for the magnetic field [50], or to develop
a new two-relaxation time based LB scheme [55] for the magnetic induction
equation for simulations of wall-bounded MHD flows. Also, adapting the LB
methods based on central moments [56, 57] that have been extended for
inhomogeneous and anisotropic flows [58, 59, 60], could enable efficient and
robust simulations of MHD flows associated with boundary layers. In addition,
performing MHD simulations by solving the magnetic induction equations within
both the liquid and the solid directly, with generally different electrical
conductivities in each domain, and then imposing the appropriate conjugate
continuity conditions for the electric, magnetic and velocity fields at their
boundary locations may be useful to consider in certain situations. With
significant increase in focus on fusion energy related research currently and
in the near future, and complemented by other applications, there is a major
potential for LB approaches for simulating and studying wall-confined MHD flow
physics with general conducting walls under a wide range of conditions.

## Acknowledgements

The second author would like to acknowledge the support of the US National
Science Foundation (NSF) for research under Grant CBET-1705630 and for support
of the development of a computer cluster infrastructure through the project
“CC* Compute: Accelerating Science and Education by Campus and Grid Computing”
under Award 2019089.

##  Appendix A Moment-based Velocity Boundary Condition Formulation

Defining the moment of order $(m+n)$ as
$\Pi_{mn}=\sum_{\alpha=0}^{8}f_{\alpha}e_{\alpha x}^{m}e_{\alpha y}^{n}$, the
standard D2Q9 lattice supports the following nine independent moments
$\Pi_{00}$, $\Pi_{10}$, $\Pi_{01}$, $\Pi_{20}$, $\Pi_{02}$, $\Pi_{11}$,
$\Pi_{21}$, $\Pi_{12}$, and $\Pi_{22}$, out of which the first six influence
the hydrodynamics directly and the remaining three are the kinetic or ghost
moments. Let’s now first consider the bottom wall at $y=-a$, for which the
unknown distribution functions are $f_{2}$, $f_{5}$, and $f_{6}$ (see Fig.
2(b)), requiring the specification of three independent moment relations. In
this regard, we present the groupings of the moments based on linearly
independent combinations of the unknowns as follows:

Combinations of unknown distributions | Moments
---|---
$f_{2}+f_{5}+f_{6}$ | $\Pi_{00},\boxed{\Pi_{01}},\Pi_{02}$
$f_{5}-f_{6}$ | $\boxed{\Pi_{10}},\Pi_{11},\Pi_{12}$
$f_{5}+f_{6}$ | $\boxed{\Pi_{20}},\Pi_{21},\Pi_{22}$

We choose one moment from each category or row to establish a linearly
independent system of equations for the unknowns. When appropriate, we choose
only the lower order, i.e., the hydrodynamic moments [40]. In this case, we
constrain the wall normal component of momentum, and wall parallel components
of the momentum and momentum flux based on the conditions at the stationary
boundary, which are boxed in the above listings. Thus, invoking the no-slip
boundary conditions on $\bm{x}_{b}$, i.e., $\bm{u}=(u_{x},u_{y})=\bm{0}$, and
taking into account $\Pi_{20}^{(1)}(\bm{x}_{b})=0$, since via a C-E analyais
it is related to $\partial_{x}u_{x}$ which vanishes on the boundary, we
specify

$\Pi_{10}(\bm{x}_{b})=0,\quad\quad\Pi_{01}(\bm{x}_{b})=0,\quad\quad\Pi_{20}(\bm{x}_{b})=\Pi_{20}^{(0)}(\bm{x}_{b})=c_{s}^{2}\rho(\bm{x}_{b}).$
(61)

Here, we need an estimate for the boundary density $\rho(\bm{x}_{b})$. Writing
$\rho(\bm{x}_{b})=\Pi_{00}(\bm{x}_{b})=\sum_{K}f_{\alpha}(\bm{x}_{b})+\sum_{U}f_{\alpha}(\bm{x}_{b})$,
where $K=\\{0,1,3,4,7,8\\}$ and $U=\\{2,5,6\\}$ are the sets of the known and
unknown directions, respectively. From $\Pi_{01}(\bm{x}_{b})=0$, it follows
that $\sum_{U}f_{\alpha}(\bm{x}_{b})=\sum_{K_{1}}f_{\alpha}(\bm{x}_{b})$,
where $K_{1}=\\{4,7,8\\}$ is a subset of the known directions. Based on these
last two equations, it follows that
$\rho=\sum_{K}f_{\alpha}(\bm{x}_{b})+2\sum_{K_{1}}f_{\alpha}(\bm{x}_{b})$. So,
we can exactly estimate the boundary density as

$\rho(\bm{x}_{b})=f_{0}(\bm{x}_{b})+f_{1}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+2\left[f_{4}(\bm{x}_{b})+f_{7}(\bm{x}_{b})+f_{8}(\bm{x}_{b})\right].$
(62)

Then, writing the left sides of each of the three sub-equations in Eq. (61)
based on their definitions in terms of the distribution functions given at the
beginning of this section and evaluated at the boundary, we get

$\displaystyle\Pi_{10}(\bm{x}_{b})$ $\displaystyle=$ $\displaystyle
f_{1}(\bm{x}_{b})-f_{3}(\bm{x}_{b})+\boxed{f_{5}(\bm{x}_{b})}-\boxed{f_{6}(\bm{x}_{b})}-f_{7}(\bm{x}_{b})+f_{8}(\bm{x}_{b})=0,$
(63a) $\displaystyle\Pi_{01}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\boxed{f_{2}(\bm{x}_{b})}-f_{4}(\bm{x}_{b})+\boxed{f_{5}(\bm{x}_{b})}+\boxed{f_{6}(\bm{x}_{b})}-f_{7}(\bm{x}_{b})-f_{8}(\bm{x}_{b})=0,$
(63b) $\displaystyle\Pi_{20}(\bm{x}_{b})$ $\displaystyle=$ $\displaystyle
f_{1}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+\boxed{f_{5}(\bm{x}_{b})}+\boxed{f_{6}(\bm{x}_{b})}+f_{7}(\bm{x}_{b})+f_{8}(\bm{x}_{b})=c_{s}^{2}\rho(\bm{x}_{b}),$
(63c)

where the variables within the boxes are the unknown quantities to be
determined, while the rest are known. Solving the above three moment-closure
equations given in Eq. (63) for the unknown or incoming distribution functions
shown within the boxes, we finally get

$\displaystyle f_{2}(\bm{x}_{b})$ $\displaystyle=$ $\displaystyle-
c_{s}^{2}\rho(\bm{x}_{b})+f_{1}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+f_{4}(\bm{x}_{b})+2\left[f_{7}(\bm{x}_{b})+f_{8}(\bm{x}_{b})\right],$
$\displaystyle f_{5}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}c_{s}^{2}\rho(\bm{x}_{b})-\left[f_{1}(\bm{x}_{b})+f_{8}(\bm{x}_{b})\right],$
$\displaystyle f_{6}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}c_{s}^{2}\rho(\bm{x}_{b})-\left[f_{3}(\bm{x}_{b})+f_{7}(\bm{x}_{b})\right].$

Similar approach can be used for obtaining the necessary relations for the
incoming distribution functions at the top wall at $y=+a$, which yield

$\rho(\bm{x}_{b})=f_{0}(\bm{x}_{b})+f_{1}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+2\left[f_{2}(\bm{x}_{b})+f_{5}(\bm{x}_{b})+f_{6}(\bm{x}_{b})\right],$

and

$\displaystyle f_{4}(\bm{x}_{b})$ $\displaystyle=$ $\displaystyle-
c_{s}^{2}\rho(\bm{x}_{b})+f_{1}(\bm{x}_{b})+f_{2}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+2\left[f_{5}(\bm{x}_{b})+f_{6}(\bm{x}_{b})\right],$
$\displaystyle f_{7}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}c_{s}^{2}\rho(\bm{x}_{b})-\left[f_{3}(\bm{x}_{b})+f_{6}(\bm{x}_{b})\right],$
$\displaystyle f_{8}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}c_{s}^{2}\rho(\bm{x}_{b})-\left[f_{1}(\bm{x}_{b})+f_{5}(\bm{x}_{b})\right].$

##  Appendix B Boundary Schemes for the Lattice Kinetic Equation for the
Magnetic Field for Moving Walls with Finite Electrical Conductivities

###  B.1 Link-based implementation

Let us consider that the boundary is moving at a velocity $U_{w}$ parallel to
the tangential direction $\hat{t}$ of the wall (see Fig. 1). Then, the
boundary conditions on the fluid velocity components normal and tangential
directions are $u_{n}(\bm{x}_{b})=0$ and $u_{t}(\bm{x}_{b})=U_{w}$,
respectively. In the following, for convenience, we take $i$ and $j$
directions to correspond to the tangential and normal directions,
respectively, i.e., $i=\hat{t}$ and $j=\hat{n}$. The imposed magnetic field
$B_{0}$ is along the normal direction, i.e., $B_{n}=B_{0}$. Then, the
equilibrium part of the anti-symmetric tensor $\Lambda_{ji}$ on moving
boundary is given by
$\Lambda_{ji}^{(0)}(\bm{x}_{b})=u_{n}(\bm{x}_{b})B_{t}(\bm{x}_{b})-B_{n}(\bm{x}_{b})u_{t}(\bm{x}_{b})=-B_{n}(\bm{x}_{b})u_{t}(\bm{x}_{b})$.
As a result, for a moving wall, from the above it follows that the equilibrium
part of the $\Lambda_{ji}$ has a non-zero contribution given by

$\Lambda_{ji}^{(0)}(\bm{x}_{b})=-B_{0}U_{w}.$ (64)

Using this expression in the two equations in Eq. (39) and simplifying under
the standard assumptions made in Sec. VI, we get

$\displaystyle g_{\alpha i}(\bm{x}_{b})+g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$ $\displaystyle 2W_{\alpha}B_{i}(\bm{x}_{b}),$ (65a)
$\displaystyle g_{\alpha i}(\bm{x}_{b})-g_{\bar{\alpha}i}(\bm{x}_{b})$
$\displaystyle\approx$ $\displaystyle-2W_{\alpha}\frac{e_{\alpha
n}}{c_{sm}^{2}}B_{0}U_{w}-2W_{\alpha}\tau_{m}\delta_{t}e_{\alpha_{k}}\partial_{k}B_{i}(\bm{x}_{b}).$
(65b)

Then, following the steps given in Sec. VI and using the above updated
equations for the moving walls, the two equations in Eq. (46) modify to

$\displaystyle\frac{1}{2}\left[g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)\right]$
$\displaystyle=$ $\displaystyle W_{\alpha}B_{i}(\bm{x}_{b}),$ (66a)
$\displaystyle-\frac{1}{2\left(\tau_{m}-\frac{1}{2}\right)\delta_{t}}\left[g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)+2W_{\alpha}\frac{e_{\alpha
n}}{c_{sm}^{2}}B_{0}U_{w}\right]$ $\displaystyle=$ $\displaystyle
W_{\alpha}\delta_{t}e_{\alpha j}\partial_{j}B_{i}(\bm{x}_{b}).$ (66b)

When these two equations in the above Eq. (66) are substituted in the
Shercliff boundary condition (see Eq. (47) given earlier) and taking into
account the weighting factor $\theta$ defined in Eq. (49), we obtain the
constraint relation between the unknown incoming distribution function
$g_{\alpha i}(\bm{x}_{f},t+\delta_{t})$ and the known outgoing post-collision
distribution function $\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)$ augmented for
the effect of the moving wall as

$(1+\theta)g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})+(1-\theta)\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)+2\theta
W_{\alpha}\frac{e_{\alpha n}}{c_{sm}^{2}}B_{0}U_{w}=0.$ (67)

This equation (Eq. (67)) can then be rewritten as follows providing the
required link-based closure for the distribution function for the magnetic
field involving a moving wall:

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=\left(\frac{\theta-1}{\theta+1}\right)\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)-\underline{\left(\frac{2\theta}{\theta+1}\right)W_{\alpha}\frac{e_{\alpha
n}}{c_{sm}^{2}}B_{0}U_{w}},$ (68)

where the underlined term expresses the effect of the shearing motion of the
wall on the magnetic distribution function. The limiting or degenerate cases
of this result can then be summarized as follows:

(a) Insulated wall: $c_{w}=\theta=0$, which corresponds to

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=-\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t),$

(b) Perfectly conducting wall: $c_{w}\;\mbox{or}\;\theta\rightarrow\infty$,
which corresponds to

$g_{\alpha
i}(\bm{x}_{f},t+\delta_{t})=+\tilde{g}_{\bar{\alpha}i}(\bm{x}_{f},t)-2W_{\alpha}\frac{e_{\alpha
n}}{c_{sm}^{2}}B_{0}U_{w},$

Clearly, when the wall is insulated, its motion does not influence the
incoming magnetic distribution function. On the other hand, the largest effect
of the wall motion is introduced when the wall is perfectly conducting. Let’s
now illustrate the above result given in Eq. (68) for the Hartmann-Couette
flow between two parallel plates separated by a length $2a$ and subjected to
an imposed magnetic field $B_{0}$ normal to the walls, where the top plate at
$y=a$ is moving at a velocity $U_{w}$, while the bottom on at $y=-a$ is at
rest. Focusing on the top wall, and using $\alpha=4$, $\bar{\alpha}=2$,
$e_{\alpha n}=e_{4y}=-1$ and $W_{\alpha}=1/6$, we get

$\displaystyle y=+a:\quad g_{4x}$ $\displaystyle=$
$\displaystyle\frac{\theta-1}{\theta+1}\tilde{g}_{2x}+\frac{\theta}{3(\theta+1)}\frac{B_{0}U_{w}}{c_{sm}^{2}},$
$\displaystyle g_{4y}$ $\displaystyle=$
$\displaystyle-\tilde{g}_{4y}+2W_{4}B_{0},$

which replaces the corresponding results given at the end of Sec. VI, while
the closure for the stationary bottom wall remains the same as before.

It may be noted that the incoming distribution functions for the velocity
field can be expressed via the standard half-way bounce back scheme augmented
with a momentum correction for the moving wall [37] given by
$f_{\alpha}(\bm{x}_{f},t+\delta_{t})=\tilde{f}_{\bar{\alpha}}(\bm{x}_{f},t)+2w_{\alpha}e_{\alpha
t}U_{w}/c_{s}^{2}$.

###  B.2 Moment-based implementation

As in the previous formulation given above, consider the Hartmann-Couette flow
bounded by conducting walls, where the top plate at $y=a$ is moving at a
velocity $U_{w}$. Then, the boundary conditions on the components of the
velocity and magnetic fields at the moving boundary is given by

$u_{x}=U_{w},\quad u_{y}=0,\quad+\frac{\partial B_{x}}{\partial
y}+\frac{B_{x}}{c_{w}a}=0,\quad B_{y}=B_{0}\quad\quad\mbox{at}\;\;y=+a.$ (69)

As before (see Sec. VII), we first write the decomposition of the tensor
component $\Lambda_{yx}$ as
$\Lambda_{yx}=\Lambda_{yx}^{(0)}+\delta\Lambda_{yx}^{(1)}$. Here, unlike in
the previous case for the moving boundary at $\bm{x}_{b}$, its equilibrium
part becomes

$\Lambda_{yx}^{(0)}(\bm{x}_{b})=u_{y}(\bm{x}_{b})B_{x}(\bm{x}_{b})-B_{y}(\bm{x}_{b})u_{x}(\bm{x}_{b})=-B_{y}(\bm{x}_{b})u_{x}(\bm{x}_{b})=-B_{0}U_{w}.$
(70)

For the non-equilibrium part, we have
$\Lambda_{yx}^{(1)}(\bm{x}_{b})=-\tau_{m}c_{sm}^{2}\partial_{y}B_{x}(\bm{x}_{b})$.
Thus, it follow that

$\Lambda_{yx}(\bm{x}_{b})=-B_{0}U_{w}-\delta_{t}\tau_{m}c_{sm}^{2}\frac{\partial
B_{x}}{\partial y}(\bm{x}_{b}),$ (71)

or

$\frac{\partial B_{x}}{\partial
y}(\bm{x}_{b})=-\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}\left[\Lambda_{yx}(\bm{x}_{b})+B_{0}U_{w}\right].$
(72)

Next, writing the definitions of the relevant moments of $g_{\alpha x}$ for
$B_{x}$ and $\Lambda_{yx}$ (see Sec. VII) along with the above Eq. (72) in the
Shercliff boundary condition $\partial B_{x}/\partial y+B_{x}/(c_{w}a)=0$ at
$\bm{x}_{b}$, we get

$-\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}\left[(g_{2x}-g_{4x})+B_{0}U_{w}\right]+\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{2x}+g_{3x}+g_{4x})=0.$
(73)

Then, solving for the unknown incoming distribution function $g_{4x}$, we
finally obtain the required moment-based closure at the boundary as

$g_{4x}=\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}+\frac{1}{c_{w}a}\right]^{-1}\left\\{\left[\frac{1}{\tau_{m}c_{sm}^{2}\delta_{t}}-\frac{1}{c_{w}a}\right]g_{2x}-\frac{1}{c_{w}a}(g_{0x}+g_{1x}+g_{3x})+\underline{\frac{B_{0}U_{w}}{\tau_{m}c_{sm}^{2}\delta_{t}}}\right\\},$
(74)

where the underlined term represents the effect of the moving wall. It may be
noted that the expression of the unknown distribution function $g_{4y}$ based
on maintaining the constraint $\sum_{\alpha=0}^{4}g_{\alpha y}=B_{0}$ remains
the same as that given in Sec. VII. The limiting special cases of Eq. (74) are
then given as follows:

$\displaystyle c_{w}=0:\quad g_{4x}$ $\displaystyle=$
$\displaystyle-(g_{0x}+g_{1x}+g_{2x}+g_{3x})\quad\quad\quad\quad\mbox{(Perfectly
insulating wall)},$ $\displaystyle c_{w}\rightarrow\infty:g_{4x}$
$\displaystyle=$ $\displaystyle
g_{2x}+B_{0}U_{w}\qquad\qquad\qquad\qquad\quad\;\;\mbox{(Perfectly conducting
wall)}.$

On the other hand, the incoming distribution functions for the velocity field
for representing the moving boundary with a velocity $U_{w}$ at $y=a$ can be
obtained by imposing the following hydrodynamic moment constraints, amending
those given in Eq. (61) in A:

$\Pi_{10}(\bm{x}_{b})=0,\quad\quad\Pi_{01}(\bm{x}_{b})=\rho(\bm{x}_{b})U_{w},\quad\quad\Pi_{20}(\bm{x}_{b})=\Pi_{20}^{(0)}(\bm{x}_{b})=\rho(\bm{x}_{b})(c_{s}^{2}+U_{w}^{2}).$

Then, following the approach presented in A, we obtain the required closure
relations that read as

$\displaystyle f_{4}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle-\rho(\bm{x}_{b})(c_{s}^{2}+U_{w}^{2})+f_{1}(\bm{x}_{b})+f_{2}(\bm{x}_{b})+f_{3}(\bm{x}_{b})+2\left[f_{5}(\bm{x}_{b})+f_{6}(\bm{x}_{b})\right],$
$\displaystyle f_{7}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\rho(\bm{x}_{b})(c_{s}^{2}+U_{w}^{2}-U_{w})-\left[f_{3}(\bm{x}_{b})+f_{6}(\bm{x}_{b})\right],$
$\displaystyle f_{8}(\bm{x}_{b})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\rho(\bm{x}_{b})(c_{s}^{2}+U_{w}^{2}+U_{w})-\left[f_{1}(\bm{x}_{b})+f_{5}(\bm{x}_{b})\right],$

where $\rho(\bm{x}_{b})$ is computed from the same expression as that given at
the end of A.

##  Appendix C Derivation of the Analytical Solution for Hartmann-Couette Flow
Bounded by Plates with Different Wall Conductance Ratios

We will now present a systematic derivation of the analytical solution for the
magnetic and velocity fields in a Hartmann-Couette flow bounded by plates with
contrasts in the wall conductance ratios between them. It should be noted that
exact solutions involving different wall conductance ratios on different
bounding wall surfaces for the coupled MHD system are generally rare. The main
objective here is to develop such an analytical solution and use it as a
benchmark case for performing a validation of the new boundary schemes
involving MHD shear flows under more general conditions on the electrical
properties of the walls in Sec. VIII. Nevertheless, we will point out later in
this section that the latter leads to an interesting way to modulate the
resulting structures of magnetic and velocity fields in a shear flow of liquid
metals and the associated peculiar flow physics.

We consider an electrically conducting fluid with conductivity $\sigma$ and
bounded by the lower and upper plates, with the wall conductance ratios of
$c_{wl}$ and $c_{wu}$, respectively, and separated by a distance of $2a$. The
fluid motion is set up due to shear generated by the steady motion of the top
plate at a velocity $U_{w}$ in the $x$ direction under an external magnetic
field $B_{0}$ applied normal to the plates in the $y$ direction, while the
bottom plate is taken to be stationary. The origin of the coordinate system is
placed midway between the two plates as in the other canonical problem shown
in Fig. 4. Then, under the assumption of one-dimensional, steady, laminar
Stokes flow, the governing MHD equations given in Eqs. (1b) and (1c),
respectively, written for the respective $x$ component simplify to

$\rho\nu\frac{d^{2}u_{x}}{dy^{2}}+\frac{B_{0}}{\mu_{m}}\frac{dB_{x}}{dy}=0.$
(75)

and

$\frac{1}{\sigma\mu_{m}}\frac{d^{2}B_{x}}{dy^{2}}+B_{0}\frac{du_{x}}{dy}=0.$
(76)

The above two equations are subject to the following no-slip velocity boundary
conditions for the stationary bottom wall and moving top wall at a velocity
$U_{w}$, as well as the Shercliff conducting wall conditions for the $x$
component of the magnetic field with $c_{wl}$ and $c_{wu}$ for the wall
conductance ratios at the lower and upper walls, respectively. These can be
written as

$\displaystyle u_{x}$ $\displaystyle=$ $\displaystyle
0,\quad\quad\quad-\frac{\partial B_{x}}{\partial
y}+\frac{B_{x}}{c_{wl}a}=0,\quad\quad\mbox{at}\;\;y=-a,$ (77a) $\displaystyle
u_{x}$ $\displaystyle=$ $\displaystyle U_{w},\quad\quad+\frac{\partial
B_{x}}{\partial y}+\frac{B_{x}}{c_{wu}a}=0,\quad\quad\mbox{at}\;\;y=+a.$ (77b)

Then, for convenience, defining the following dimensionless variables natural
for this flow problem

$\xi=\frac{y}{a},\qquad\qquad\mbox{Ha}=B_{0}a\left(\frac{\sigma}{\rho\nu}\right)^{1/2},\qquad\qquad
u_{x}^{*}=\frac{u_{x}}{U_{w}},\qquad\qquad
B_{x}^{*}=\frac{B_{x}}{U_{w}(\rho\nu\sigma)^{1/2}\mu_{m}},$ (78)

Equations (75) and (76) can be rewritten in their more compact non-dimensional
forms as follows:

$\displaystyle\frac{d^{2}u_{x}^{*}}{d\xi^{2}}+\mbox{Ha}\frac{dB_{x}^{*}}{d\xi}$
$\displaystyle=$ $\displaystyle 0,$ (79a)
$\displaystyle\frac{d^{2}B_{x}^{*}}{d\xi^{2}}+\mbox{Ha}\frac{du_{x}^{*}}{d\xi}$
$\displaystyle=$ $\displaystyle 0.$ (79b)

and the boundary conditions given in Eq. (77) become

$\displaystyle u_{x}^{*}$ $\displaystyle=$ $\displaystyle
0,\quad\quad-\frac{\partial
B_{x}^{*}}{\partial\xi}+\frac{B_{x}^{*}}{c_{wl}}=0,\quad\quad\mbox{at}\;\;\xi=-1,$
(80a) $\displaystyle u_{x}^{*}$ $\displaystyle=$ $\displaystyle
1,\quad\quad+\frac{\partial
B_{x}^{*}}{\partial\xi}+\frac{B_{x}^{*}}{c_{wu}}=0,\quad\quad\mbox{at}\;\;\xi=+1.$
(80b)

In order to solve the above set of coupled differential equations for
$u_{x}^{*}(\xi)$ and $B_{x}^{*}(\xi)$, by exploiting their special structures,
we introduce the following transformation variables (sometimes referred to as
the Elsasser variables [61]):

$V_{1}^{*}=u_{x}^{*}+B_{x}^{*},\quad\quad V_{2}^{*}=u_{x}^{*}-B_{x}^{*}.$ (81)

As a result, combining Eqs. (79a) and (79b) via taking their sum and
difference and using Eq. (81), we get the following intermediate de-coupled
system:

$\displaystyle\frac{d^{2}V_{1}^{*}}{d\xi^{2}}+\mbox{Ha}\frac{dV_{1}^{*}}{d\xi}$
$\displaystyle=$ $\displaystyle 0,$ (82)
$\displaystyle\frac{d^{2}V_{2}^{*}}{d\xi^{2}}-\mbox{Ha}\frac{dV_{2}^{*}}{d\xi}$
$\displaystyle=$ $\displaystyle 0.$ (83)

These two equations can be readily solved, which can be written as

$V_{1}^{*}(\xi)=m+n\exp(-\mbox{Ha}\xi),\quad\quad
V_{2}^{*}(\xi)=m^{\prime}+n^{\prime}\exp(+\mbox{Ha}\xi),$ (84)

where $m$, $m^{\prime}$, $n$ and $n^{\prime}$ are arbitrary integration
constants. We can then recover the solutions in terms of the original
variables $u_{x}^{*}(\xi)$ and $B_{x}^{*}(\xi)$ via using
$u_{x}^{*}=(V_{1}^{*}+V_{2}^{*})/2$ and $B_{x}^{*}=(V_{1}^{*}-V_{2}^{*})/2$,
which yield

$\displaystyle u_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle(p+p^{\prime})+q\exp(-\mbox{Ha}\xi)+q^{\prime}\exp(+\mbox{Ha}\xi),$
(85a) $\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle(p-p^{\prime})+q\exp(-\mbox{Ha}\xi)-q^{\prime}\exp(+\mbox{Ha}\xi),$
(85b)

where the redefined integration constants are given by $p=m/2$,
$p^{\prime}=m^{\prime}/2$, $q=n/2$ and $q^{\prime}=n^{\prime}/2$. Rewriting
the last two equations (Eqs. (85a) and (85b)) further in terms of the
hyperbolic functions for the bounded flow problem under consideration, and
using $r=q+q^{\prime}$ and $r^{\prime}=q-q^{\prime}$ for convenience, we
finally get the general solutions for the velocity and magnetic fields as

$\displaystyle u_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle(p+p^{\prime})+r\cosh(\mbox{Ha}\xi)-r^{\prime}\sinh(\mbox{Ha}\xi),$
(86a) $\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle(p-p^{\prime})+r^{\prime}\cosh(\mbox{Ha}\xi)-r\sinh(\mbox{Ha}\xi).$
(86b)

The four constants $p$, $p^{\prime}$, $r$ and $r^{\prime}$ can now be obtained
from the four boundary conditions given above. Now, using the velocity field
given above in the no-slip velocity boundary conditions at the bottom and top
plates given in Eq. (80), it follows that

$r=\frac{1/2-(p+p^{\prime})}{\cosh(\mbox{Ha})},\quad\quad
r^{\prime}=-\frac{1}{2\sinh(\mbox{Ha})}.$ (87)

Then, using the magnetic field given in Eq. (86b) and substituting the
Shercliff boundary conditions at the bottom and top plates, respectively (see
Eq. (80)), and using the expressions for the $r$ and $r^{\prime}$ given above
in Eq. (87) and simplifying, we get the following coupled set of two equations
for $p$ and $p^{\prime}$ as

$\displaystyle-\lambda_{1}(p+p^{\prime})+(p-p^{\prime})$ $\displaystyle=$
$\displaystyle\gamma_{1},$ (88a)
$\displaystyle-\lambda_{2}(p+p^{\prime})+(p-p^{\prime})$ $\displaystyle=$
$\displaystyle\gamma_{2},$ (88b)

where

$\displaystyle\lambda_{1}$ $\displaystyle=$ $\displaystyle
c_{wl}\mbox{Ha}+\tanh(\mbox{Ha}),\quad\quad\gamma_{1}=\frac{1}{2}\left[\frac{1}{\tanh(\mbox{Ha})}-\tanh(\mbox{Ha})\right],$
(89a) $\displaystyle\lambda_{2}$ $\displaystyle=$ $\displaystyle
c_{wu}\mbox{Ha}+\tanh(\mbox{Ha}),\quad\quad\gamma_{2}=\frac{1}{2}\left[\frac{1}{\tanh(\mbox{Ha})}+\tanh(\mbox{Ha})\right]+c_{wu}\mbox{Ha}.$
(89b)

Then, solving for $(p+p^{\prime})$ and $(p-p^{\prime})$ from the two equations
given in Eq. (88), we obtain

$(p+p^{\prime})=\frac{\gamma_{2}-\gamma_{1}}{\lambda_{1}+\lambda_{2}},\quad\quad(p-p^{\prime})=\frac{\lambda_{1}\gamma_{2}+\lambda_{2}\gamma_{1}}{\lambda_{1}+\lambda_{2}}.$
(90)

Using these results for $(p+p^{\prime})$ and $(p-p^{\prime})$ (see Eq. (90))
in Eqs. (86a) and (86b), and along with the expressions for $r$ and
$r^{\prime}$ (see Eq. (87)) and simplifying and rearranging, we finally obtain
the analytical solutions for the velocity and magnetic fields in the Hartmann-
Couette flow bounded by plates with contrasts in electrical conductivities,
which can be expressed as follows:

$\displaystyle u_{x}^{*}(\xi)$ $\displaystyle=$ $\displaystyle
K_{1}+\left(\frac{1}{2}-K_{1}\right)\frac{\cosh(\mbox{Ha}\xi)}{\cosh(\mbox{Ha})}+\frac{1}{2}\frac{\sinh(\mbox{Ha}\xi)}{\sinh(\mbox{Ha})},$
(91a) $\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$ $\displaystyle
K_{2}-\left(\frac{1}{2}-K_{1}\right)\frac{\sinh(\mbox{Ha}\xi)}{\cosh(\mbox{Ha})}-\frac{1}{2}\frac{\cosh(\mbox{Ha}\xi)}{\sinh(\mbox{Ha})},$
(91b)

where

$\displaystyle K_{1}$ $\displaystyle=$
$\displaystyle\frac{\tanh(\mbox{Ha})+c_{wu}\mbox{Ha}}{2\tanh(\mbox{Ha})+(1+\phi)c_{wu}\mbox{Ha}},$
(92a) $\displaystyle K_{2}$ $\displaystyle=$
$\displaystyle\frac{1+((1+\phi)/2)c_{wu}\mbox{Ha}\left[\tanh(\mbox{Ha})+1/\tanh(\mbox{Ha})\right]+(c_{wu}\mbox{Ha})^{2}\phi}{2\tanh(\mbox{Ha})+(1+\phi)c_{wu}\mbox{Ha}},$
(92b) $\displaystyle\phi$ $\displaystyle=$
$\displaystyle\frac{c_{wl}}{c_{wu}}.$ (92c)

Here, the value of the factor $\phi$ represents the contrast between the wall
conductance ratio at the bottom and top plates and thus the coefficients
$K_{1}$ and $K_{2}$ parameterize the influence of electrical properties of the
bounding walls. Equations (91a) and (91b) represent new analytical solutions
for the shear driven MHD flow between two plates with different electrical
conductivities and can be used for validating computational techniques for MHD
flows involving moving walls. Note that for the special case with the same
conductance ratios for both the walls ($c_{wl}=c_{wu}=c_{w}$), i.e., $\phi=1$,
the coefficient $K_{1}=1/2$ and the above analytical solutions simplify to

$\displaystyle u_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+\frac{\sinh(\mbox{Ha}\xi)}{\sinh(\mbox{Ha})}\right],$
(93a) $\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1+c_{w}\mbox{Ha}\left[\tanh(\mbox{Ha})+1/\tanh(\mbox{Ha})\right]+(c_{w}\mbox{Ha})^{2}}{2\left[\tanh(\mbox{Ha})+c_{w}\mbox{Ha}\right]}-\frac{1}{2}\frac{\cosh(\mbox{Ha}\xi)}{\sinh(\mbox{Ha})}.$
(93b)

This leads to the important conclusion that while the induced magnetic field
$B_{x}^{*}$ is always affected by the value of the wall conductance ratios,
whether it is the same or different for the two walls, by contrast, the
velocity field $u_{x}^{*}$ is sensitive to the electrical properties of the
walls only if there is a contrast between the values of the wall conductance
ratio at the bottom and top plates. This finding provides a potential
mechanism to control the flow and transport in MHD shear flows by adjusting
the electrical properties of the bounding plates. Also, as a further
simplification, considering both the walls to be insulated, i.e., $c_{w}=0$,
when we have $K_{2}=1/(2\tanh(\mbox{Ha}))$, the structure of the velocity and
magnetic fields reduce further to

$\displaystyle u_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+\frac{\sinh(\mbox{Ha}\xi)}{\sinh(\mbox{Ha})}\right],$
(94a) $\displaystyle B_{x}^{*}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2\tanh(\mbox{Ha})}\left[1-\frac{\cosh(\mbox{Ha}\xi)}{\cosh(\mbox{Ha})}\right].$
(94b)

## References

  * Shercliff [1956] J. Shercliff, The flow of conducting fluids in circular pipes under transverse magnetic fields, Journal of Fluid Mechanics 1, 644 (1956).
  * Dellar [2002] P. J. Dellar, Lattice kinetic schemes for magnetohydrodynamics, Journal of Computational Physics 179, 95 (2002).
  * Hartmann [1937] J. Hartmann, Hg–dynamics i. Theory of the laminar flow of an electrically conductive liquid in a homogeneous magnetic field, Det Kgl Danske Videnskabernes Selskkab Math–fys Medd 15, 1 (1937).
  * Hartmann and Lazarus [1937] J. Hartmann and F. Lazarus, Hg–dynamics ii. Experimental investigations on the flow of mercury in a homogeneous magnetic field, Det Kgl Danske Videnskabernes Selskkab Math–fys Medd 15, 1 (1937).
  * Alfvén [1942] H. Alfvén, Existence of electromagnetic-hydrodynamic waves, Nature 150, 405 (1942).
  * Lundquist [1949] S. Lundquist, Experimental investigations of magneto-hydrodynamic waves, Physical Review 76, 1805 (1949).
  * Elsasser [1955] W. M. Elsasser, Hydromagnetism. i. a review, American Journal of Physics 23, 590 (1955).
  * Molokov  _et al._ [2007] S. S. Molokov, R. Moreau, and H. K. Moffatt,  _Magnetohydrodynamics: Historical evolution and trends_ , Vol. 80 (Springer Science & Business Media, 2007).
  * Davidson [1999] P. Davidson, Magnetohydrodynamics in materials processing, Annual Review of Fluid Mechanics 31, 273 (1999).
  * Müller and Bühler [2013] U. Müller and L. Bühler,  _Magnetofluiddynamics in channels and containers_ (Springer Science & Business Media, 2013).
  * Morley  _et al._ [2000] N. B. Morley, S. Smolentsev, L. Barleon, I. R. Kirillov, and M. Takahashi, Liquid magnetohydrodynamics—recent progress and future directions for fusion, Fusion Engineering and design 51, 701 (2000).
  * Smolentsev  _et al._ [2010] S. Smolentsev, R. Moreau, L. Bühler, and C. Mistrangelo, MHD thermofluid issues of liquid-metal blankets: phenomena and advances, Fusion Engineering and Design 85, 1196 (2010).
  * Rhodes  _et al._ [2018] T. J. Rhodes, S. Smolentsev, and M. Abdou, Magnetohydrodynamic pressure drop and flow balancing of liquid metal flow in a prototypic fusion blanket manifold, Physics of Fluids 30, 057101 (2018).
  * Smolentsev [2021] S. Smolentsev, Physical background, computations and practical issues of the magnetohydrodynamic pressure drop in a fusion liquid metal blanket, Fluids 6, 110 (2021).
  * Mistrangelo  _et al._ [2021] C. Mistrangelo, L. Bühler, S. Smolentsev, V. Klüber, I. Maione, and J. Aubert, MHD flow in liquid metal blankets: Major design issues, MHD guidelines and numerical analysis, Fusion Engineering and Design 173, 112795 (2021).
  * Shercliff [1965] J. A. Shercliff,  _A Textbook of Magnetohydrodynamics_ (Pergamon Press, Oxford, 1965).
  * Moreau [1990] R. J. Moreau,  _Magnetohydrodynamics_ , Vol. 3 (Springer Science & Business Media, 1990).
  * Davidson [2001] P. A. Davidson,  _An introduction to magnetohydrodynamics_ (Cambridge University Press, 2001).
  * Blishchik  _et al._ [2021] A. Blishchik, M. van der Lans, and S. Kenjereš, An extensive numerical benchmark of the various magnetohydrodynamic flows, International Journal of Heat and Fluid Flow 90, 108800 (2021).
  * Blishchik and Kenjereš [2021] A. Blishchik and S. Kenjereš, Observation of a novel flow regime caused by finite electric wall conductance in an initially turbulent magnetohydrodynamic duct flow, Physical Review E 104, L013101 (2021).
  * Chang and Lundgren [1961] C. C. Chang and T. S. Lundgren, Duct flow in magnetohydrodynamics, Zeitschrift für angewandte Mathematik und Physik ZAMP 12, 100 (1961).
  * Hunt [1965] J. Hunt, Magnetohydrodynamic flow in rectangular ducts, Journal of fluid mechanics 21, 577 (1965).
  * Yu [1970] C. Yu, Convective magnetohydrodynamic flow in a vertical channel, Applied Scientific Research 22, 127 (1970).
  * Singh and Lal [1984] B. Singh and J. Lal, Finite element method for unsteady MHD flow through pipes with arbitrary wall conductivity, International journal for numerical methods in fluids 4, 291 (1984).
  * Sterl [1990] A. Sterl, Numerical simulation of liquid-metal MHD flows in rectangular ducts, Journal of Fluid Mechanics 216, 161 (1990).
  * Mistrangelo [2006] C. Mistrangelo,  _Three-dimensional MHD flow in sudden expansions_ (Forschungszentrum Karlsruhe, 2006).
  * Vantieghem  _et al._ [2009] S. Vantieghem, X. Albets-Chico, and B. Knaepen, The velocity profile of laminar MHD flows in circular conducting pipes, Theoretical and computational fluid dynamics 23, 525 (2009).
  * Mistrangelo and Bühler [2017] C. Mistrangelo and L. Bühler, Magnetohydrodynamic flows in liquid metal blankets for fusion reactors, PAMM 17, 115 (2017).
  * Arlt  _et al._ [2017] T. Arlt, J. Priede, and L. Bühler, The effect of finite-conductivity Hartmann walls on the linear stability of Hunt’s flow, Journal of Fluid Mechanics 822, 880 (2017).
  * McNamara and Zanetti [1988] G. R. McNamara and G. Zanetti, Use of the Boltzmann equation to simulate lattice-gas automata, Physical review letters 61, 2332 (1988).
  * Higuera and Jiménez [1989] F. J. Higuera and J. Jiménez, Boltzmann approach to lattice gas simulations, EPL (Europhysics Letters) 9, 663 (1989).
  * He and Luo [1997] X. He and L.-S. Luo, Theory of the lattice boltzmann method: From the Boltzmann equation to the lattice boltzmann equation, Physical Review E 56, 6811 (1997).
  * Dellar [2013a] P. J. Dellar, An interpretation and derivation of the lattice Boltzmann method using strang splitting, Computers & Mathematics with Applications 65, 129 (2013a).
  * Benzi  _et al._ [1992] R. Benzi, S. Succi, and M. Vergassola, The lattice Boltzmann equation: theory and applications, Physics Reports 222, 145 (1992).
  * Lallemand  _et al._ [2021] P. Lallemand, L.-s. Luo, M. Krafczyk, and W.-A. Yong, The lattice boltzmann method for nearly incompressible flows, Journal of Computational Physics 431, 109713 (2021).
  * Krüger  _et al._ [2017] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The lattice Boltzmann method, Springer International Publishing 10, 978 (2017).
  * Ladd [1994] A. J. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. part 1. theoretical foundation, Journal of fluid mechanics 271, 285 (1994).
  * Bouzidi  _et al._ [2001] M. Bouzidi, M. Firdaouss, and P. Lallemand, Momentum transfer of a boltzmann-lattice fluid with boundaries, Physics of fluids 13, 3452 (2001).
  * Noble  _et al._ [1995] D. R. Noble, S. Chen, J. G. Georgiadis, and R. O. Buckius, A consistent hydrodynamic boundary condition for the lattice Boltzmann method, Physics of Fluids 7, 203 (1995).
  * Bennett [2010] S. Bennett,  _A lattice Boltzmann model for diffusion of binary gas mixtures_ , Ph.D. thesis, University of Cambridge (2010).
  * Bennett  _et al._ [2012] S. Bennett, P. Asinari, and P. J. Dellar, A lattice Boltzmann model for diffusion of binary gas mixtures that includes diffusion slip, International journal for numerical methods in fluids 69, 171 (2012).
  * Reis and Dellar [2012] T. Reis and P. J. Dellar, Lattice Boltzmann simulations of pressure-driven flows in microchannels using navier–maxwell slip boundary conditions, Physics of Fluids 24, 112001 (2012).
  * Dellar [2013b] P. J. Dellar, Moment-based boundary conditions for lattice boltzmann magnetohydrodynamics, in  _Numerical Mathematics and Advanced Applications 2011_ (Springer, 2013) pp. 83–90.
  * Chen  _et al._ [1991] S. Chen, H. Chen, D. Martnez, and W. Matthaeus, Lattice Boltzmann model for simulation of magnetohydrodynamics, Physical Review Letters 67, 3776 (1991).
  * Martínez  _et al._ [1994] D. O. Martínez, S. Chen, and W. H. Matthaeus, Lattice Boltzmann magnetohydrodynamics, Physics of plasmas 1, 1850 (1994).
  * Chen  _et al._ [1992] S. Chen, D. O. Martinez, W. Matthaeus, and H. Chen, Magnetohydrodynamics computations with lattice gas automata, Journal of statistical physics 68, 533 (1992).
  * Succi  _et al._ [1991] S. Succi, M. Vergassola, and R. Benzi, Lattice Boltzmann scheme for two-dimensional magnetohydrodynamics, Physical Review A 43, 4521 (1991).
  * Mendoza and Munoz [2008] M. Mendoza and J. Munoz, Three-dimensional lattice Boltzmann model for magnetic reconnection, Physical Review E 77, 026713 (2008).
  * Pattison  _et al._ [2008] M. Pattison, K. Premnath, N. Morley, and M. Abdou, Progress in lattice Boltzmann methods for magnetohydrodynamic flows relevant to fusion applications, Fusion Engineering and Design 83, 557 (2008).
  * Dellar [2009] P. J. Dellar, Moment equations for magnetohydrodynamics, Journal of Statistical Mechanics: Theory and Experiment 2009, P06003 (2009).
  * Dellar [2013c] P. J. Dellar, Lattice Boltzmann magnetohydrodynamics with current-dependent resistivity, Journal of Computational Physics 237, 115 (2013c).
  * De Rosis  _et al._ [2021] A. De Rosis, R. Liu, and A. Revell, One-stage simplified lattice Boltzmann method for two-and three-dimensional magnetohydrodynamic flows, Physics of Fluids 33, 085114 (2021).
  * Chapman  _et al._ [1990] S. Chapman, T. G. Cowling, and D. Burnett,  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ (Cambridge university press, 1990).
  * Qian  _et al._ [1992] Y.-H. Qian, D. d’Humières, and P. Lallemand, Lattice BGK models for navier-stokes equation, EPL (Europhysics Letters) 17, 479 (1992).
  * Ginzburg  _et al._ [2008] I. Ginzburg, F. Verhaeghe, and D. d’Humieres, Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions, Communications in computational physics 3, 427 (2008).
  * Geier  _et al._ [2006] M. Geier, A. Greiner, and J. G. Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Physical Review E 73, 066705 (2006).
  * Premnath and Banerjee [2011] K. N. Premnath and S. Banerjee, On the three-dimensional central moment lattice Boltzmann method, Journal of Statistical Physics 143, 747 (2011).
  * Yahia and Premnath [2021] E. Yahia and K. N. Premnath, Central moment lattice Boltzmann method on a rectangular lattice, Physics of Fluids 33, 057110 (2021).
  * Yahia  _et al._ [2021] E. Yahia, W. Schupbach, and K. N. Premnath, Three-dimensional central moment lattice Boltzmann method on a cuboid lattice for anisotropic and inhomogeneous flows, Fluids 6, 326 (2021).
  * Yahia [2021] E. Yahia,  _Efficient central moment lattice Boltzmann algorithms based on rectangular and cuboid lattices for inhomogeneous and anisotropic flows_ , Ph.D. thesis, University of Colorado Denver (2021).
  * Elsasser [1950] W. M. Elsasser, The hydromagnetic equations, Physical Review 79, 183 (1950).