# Magneto-Stokes Flow in a Shallow Free-Surface Annulus

Cy S. David1<EMAIL_ADDRESS>Eric W. Hester2  Yufan Xu1,3  Jonathan M.
Aurnou1  1Department of Earth, Planetary, and Space Sciences, University of
California, Los Angeles, CA 90095, USA 2Department of Mathematics, University
of California, Los Angeles, CA 90095, USA 3Princeton Plasma Physics
Laboratory, Princeton, NJ 08540, USA

###### Abstract

We analyse a magnetohydrodynamic flow inspired by the kinematic reversibility
of viscous Taylor-Couette flows. The system considered here shares the
cylindrical-annular geometry of the Taylor-Couette cell, but uses applied
electromagnetic forces to drive “magneto-Stokes” flow in a shallow, free-
surface layer of electrolyte. An analytical solution is presented and
validated with coupled laboratory and numerical experiments. The dominant
balance of Lorentz forcing and basal viscous drag reproduces the kinematic
reversibility observed by G.I. Taylor with precise electromagnetic control.
Induced fluid deformation may be undone by simply reversing the polarity of
electric current through the system. We illustrate this analogy with theory
and experiment, and we draw a further connection to potential flow using the
Hele-Shaw approximation. The stability and controllability of the magneto-
Stokes system make it an attractive tool for investigating shear flows in a
variety of settings from industrial to astrophysical. In addition, the set-
up’s simplicity and robustness make magneto-Stokes flow a good candidate for
PIV calibration and for educational demonstrations of magnetohydrodynamics,
boundary layers, and flow transition.

###### keywords:

Low-Reynolds-number flows, MHD and Electrohydrodynamics, Hele-Shaw flows

##  1 Introduction

The year 2023 marks the 100th anniversary of G.I. Taylor’s seminal publication
(Taylor, 1923) on flow confined between two rotating concentric cylinders,
termed Taylor-Couette flow in his honour. We celebrate the centennial of this
work by revisiting what is arguably the most famous experiment on the Taylor-
Couette device: Taylor’s film demonstration (Taylor, 1967) of kinematic
reversibility—the remarkable time-independence of Stokes flows that allows
deformations to be undone by a reversal of the applied forces.

We consider kinematically-reversible Taylor-Couette-like flow in a bench-top
device with no moving solid components, using electromagnetic body forces
rather than viscous traction to drive motion. In this magnetohydrodynamic
(MHD) modification of the Taylor-Couette cell, the usually rotating sidewalls
of the cylindrical annulus are fixed and made electrically-conducting. The
base is kept electrically-insulating, while the lid is removed to allow free-
surface flow. An applied axial magnetic field and radial electric current
drive azimuthal flow of an electrolyte, for which (1) magnetic induction is
small compared to magnetic diffusion (low magnetic Reynolds number, Rm), (2)
magnetic drag is small relative to viscous forces (low Hartmann number, Ha),
and (3) inertia is small compared to viscous forces (low Reynolds number,
$\Rey$). Under these conditions, the azimuthal momentum balance is dominated
by the Lorentz force and viscous drag due to the channel sidewalls and base.
We term the resulting circulatory motion “annular magneto-Stokes flow.”

Similar MHD flows through cylindrical annular ducts with conducting sidewalls
and axial magnetic field have been considered since the works of Early & Dow
(1950), Anderson et al. (1959), and Braginsky (1959). Early modelling efforts
focus on the high-Hartmann number limit for analytical convenience (Hunt &
Stewartson, 1965) and relevance to liquid metal flows (Baylis & Hunt, 1971).
Later works have numerically surveyed a broader range of Hartmann numbers in
closed annular ducts (Poyé et al., 2020), analysed the stability of Hartmann
layers (Moresco & Alboussire, 2004), and studied the effects of modified
electric boundary conditions (Stelzer et al., 2015a, b) in similar
cylindrical-annular MHD flows.

In contrast to the above studies relevant to large-scale liquid metal flows,
recent applications of MHD to microfluidic systems have generated interest in
low-Hartmann number solutions (West et al., 2003; Khal’zov & Smolyakov, 2006)
appropriate for electrolyte flows. Digilov (2007) and Pérez-Barrera et al.
(2015, 2016) provide analytical solutions for the steady annular magneto-
Stokes velocity profile over radius, $u_{\theta}(r)$; Valenzuela-Delgado et
al. (2018a, b) use a (semi-analytical) Galerkin approximation to predict
steady, axisymmetric flow over radius and depth, $u_{\theta}(r,z)$.

We extend and attempt to unify these efforts by providing the fully-analytical
representation of $u_{\theta}(r,z,t)$, distinguishing deep and shallow-layer
regimes of annular magneto-Stokes flow and analysing the time-dependent
behaviour of small, but finite Reynolds number flows. These theoretical
developments are detailed in Section 2. The solution for shallow-layer flow is
validated with laboratory and numerical experiments in Sections 3 and 4. In
Section 5, we make explicit the analogy between magneto-Stokes flow in a
shallow-layer annulus and viscous Taylor-Couette flow, replicating Taylor’s
famous kinematic reversibility result in the laboratory. We then draw a
connection to another classical property of Stokes flow, involving the Hele-
Shaw approximation and potential flow around an obstacle. Applications to
engineering, astrophysical fluid dynamics, and education are discussed in
Section 6.

##  2 Theory

Figure 1: (a) Diagram of the magneto-Stokes system. A power supply controls
the electric current $I$ through the fluid layer. A configuration of radially
outwards ($+\boldsymbol{e}_{r}$) current density $\boldsymbol{J}$ and
downwards ($-\boldsymbol{e}_{z}$) magnetic field $\boldsymbol{B}$ produces an
azimuthal ($+\boldsymbol{e}_{\theta}$) Lorentz force $\boldsymbol{F_{L}}$ on
the fluid. (b) Photograph of the channel used in laboratory experiments, with
flow visualised by blue dye. The channel rests atop a wooden case of permanent
magnets, which is replaced by a solenoid electromagnet (not pictured) for
Cases I-IV and KR discussed in Section 3.

###  2.1 Axisymmetric governing equations

We consider a free-surface layer of conducting fluid of depth $h$ in the
annular gap between two cylindrical electrodes of radius $r_{i}$ and $r_{o}$
($r_{i}<r_{o}$). A controlled current $I$ runs through the fluid from inner to
outer electrode, and the entire annulus is subject to a vertical, imposed
magnetic field $\boldsymbol{B}=-B_{0}\boldsymbol{e_{z}}$. Figure 1(a) shows a
schematic of the annular channel with the imposed magnetic field and current.
For a low-conductivity fluid like saltwater and small $B_{0}$, the magnetic
field is quasi-static (e.g., Davidson 2016; Knaepen & Moreau 2008; Favier et
al. 2011; Verma 2017) and the total Lorentz force on the fluid is

$\boldsymbol{F_{L}}=\boldsymbol{J}\times\boldsymbol{B},$ (1)

where the current density $\boldsymbol{J}$ may be expressed using Ohm’s Law,

$\boldsymbol{J}=\sigma(\boldsymbol{E}+\boldsymbol{u}\times\boldsymbol{B}),$
(2)

in which $\sigma$ is the electrical conductivity of the fluid,
$\boldsymbol{E}$ is the imposed electric field, and $\boldsymbol{u}$ is the
fluid velocity. Substituting (2) into (1) and adopting cylindrical coordinates
($r$,$\theta$,$z$) shows that the Lorentz force consists of a driving
component resulting from the applied fields and a magnetic drag component that
damps fluid motion across magnetic field lines:

$\boldsymbol{F_{L}}=\sigma(\boldsymbol{E}\times\boldsymbol{B})-\sigma
B_{0}^{2}(u_{r}\boldsymbol{e}_{r}+u_{\theta}\boldsymbol{e}_{\theta}).$ (3)

Letting $U$ denote a characteristic velocity scale, the magnitude of the
magnetic drag ($\sim\sigma B_{0}^{2}U$) may be compared to that of the viscous
drag ($\sim\varrho\nu U/h^{2}$) by means of the Hartmann number:

$\textit{Ha}=\sqrt{\frac{\sigma B_{0}^{2}U}{\varrho\nu
U/h^{2}}}=B_{0}h\sqrt{\frac{\sigma}{\varrho\nu}},$ (4)

where $\varrho$ and $\nu$ are the density and kinematic viscosity of the
fluid, respectively.

In our experiments, $\textit{Ha}\lesssim 10^{-2}$, and magnetic drag is thus
four orders of magnitude smaller than viscous forces. In the low Ha limit, the
only component of current density $\boldsymbol{J}$ that makes a significant
contribution to the Lorentz force is $\sigma\boldsymbol{E}$, which may be
determined purely from the electric boundary conditions. A DC power supply can
control the voltage across the sidewalls such that the total current $I$ is
set to a desired value at each time $t$. In this case, current rather than
voltage is used as a control parameter, and current density is appropriately
expressed as

$\boldsymbol{J}=\sigma\boldsymbol{E}=\frac{I(t)}{2\pi rh}\boldsymbol{e_{r}},$
(5)

neglecting any fringing of electric field lines due to the finite fluid depth
(i.e., assuming that $\partial_{z}\boldsymbol{E}=0$).

The resulting Lorentz force,

$\boldsymbol{F_{L}}=\frac{B_{0}I(t)}{2\pi rh}\boldsymbol{e_{\theta}},$ (6)

drives circulatory flow governed by the incompressible equations of motion:

$\partial_{t}\boldsymbol{u}+\boldsymbol{u}\bcdot\bnabla\boldsymbol{u}=-\frac{1}{\varrho}\bnabla
p+\nu\nabla^{2}\boldsymbol{u}+\frac{1}{\varrho}\boldsymbol{F_{L}},\quad\bnabla\bcdot\boldsymbol{u}=0,$
(7)

where $p$ is the deviation in pressure from the static pressure field $\varrho
g(h-z)$.

Symbol |  Description
---|---
$r_{i},r_{o}$ |  Inner and outer sidewall radii [m]
$h$ |  Depth of fluid [m]
$g$ |  Gravitational acceleration [m/s2]
$\nu$ |  Kinematic viscosity [m2/s]
$\varrho$ |  Density [kg/m3]
$\sigma$ |  Electrical conductivity [S/m]
$\mu_{0}$ |  Magnetic permeability of free space [kg m s-2 A-2]
$I_{0}$ |  Maximum electric current amplitude [A]
$B_{0}$ |  Magnetic flux density (magnitude) [T]
$d$ |  Diameter of obstacle (Section 5.2) [m]
$\delta_{95,i},\delta_{95,o}$ |  95% thickness of inner and outer sidewall boundary layers [m], defined in (27)
Table 1: Dimensional parameter definitions. Values are given in Section 3.

We scale radial distances $r$ by the outer sidewall radius $r_{o}$, vertical
distances $z$ by the fluid depth $h$, and electric current $I$ by its maximum
value $I_{0}$, defining the nondimensional quantities $\rho=r/r_{o}$,
$\zeta=z/h$, and $\kappa=I/I_{0}$. A velocity scale $U$ may be found by
balancing the basal viscous drag and Lorentz forces at the surface mid-gap,
$z=h$, $r=(r_{i}+r_{o})/2$:

$2\frac{\nu U}{h^{2}}\sim\frac{B_{0}I_{0}}{2\pi h\varrho(r_{i}+r_{o})/2}.$ (8)

The factor of $2$ before ${\nu U}/{h^{2}}$ reflects the effect of the free-
surface boundary condition, and is included to ensure that scaled velocities
approach unity at the surface mid-gap for $h\ll r_{o}$. The resulting surface
mid-gap velocity scale,

$U=\frac{B_{0}I_{0}h}{2\pi\nu\varrho(r_{i}+r_{o})},$ (9)

is used to produce the nondimensional velocity
$\boldsymbol{\upsilon}=\upsilon_{\rho}\boldsymbol{e_{r}}+\upsilon_{\theta}\boldsymbol{e_{\theta}}+(h/r_{o})\upsilon_{\zeta}\boldsymbol{e_{z}}$,
with components $\upsilon_{\rho}=u_{r}/U$, $\upsilon_{\theta}=u_{\theta}/U$,
and $\upsilon_{\zeta}=(r_{o}/h)u_{z}/U$. An inertial scale $\varrho{U}^{2}$ is
used to nondimensionalise reduced pressure as $\Pi=p/(\varrho{U}^{2})$. Time
is nondimensionalised as $\tau=t/T$ by a surface mid-gap advective timescale:

$T=\frac{r_{i}+r_{o}}{U},$ (10)

where order-unity numerical factors have been dropped for notational
convenience. In sum, our nondimensionalisation makes the mapping

$(\boldsymbol{u},p,I)(r,\theta,z,t)\mapsto(\boldsymbol{\upsilon},\Pi,\kappa)(\rho,\theta,\zeta,\tau).$
(11)

We introduce an a priori control Reynolds number $\Rey$ that is proportional
to the ratio of viscous to advective timescales:

$\Rey=\frac{h^{2}/\nu}{T}=\frac{B_{0}I_{0}h^{3}}{2\pi\nu^{2}\varrho(r_{i}+r_{o})^{2}}.$
(12)

The motivation for our choice of Reynolds number and velocity scale will
become apparent upon solving the shallow limit of the governing equations.
Using our definition of $\Rey$, we may define a magnetic Reynolds number Rm,
which compares advective effects to magnetic rather than viscous diffusion:

$\textit{Rm}=\frac{h^{2}/\eta}{T}=\frac{\nu}{\eta}\Rey,$ (13)

where $\eta=(\mu_{0}\sigma)^{-1}$ is the fluid’s magnetic diffusivity and
$\mu_{0}$ is the magnetic permeability of free space. For $\textit{Rm}\ll 1$,
magnetic diffusion dominates over advection, and the quasi-static description
of the Lorentz force (6) is valid (e.g., Davidson 2016; Knaepen & Moreau 2008;
Favier et al. 2011; Verma 2017). As defined, $\textit{Rm}\lesssim 10^{-10}$
for our experiments.

Also relevant are the radius ratio $\mathcal{R}$ and aspect ratio
$\mathcal{H}$, defined as

$\mathcal{R}=r_{i}/r_{o}\quad\text{and}\quad\mathcal{H}=h/r_{o}.$ (14)

Symbols for all dimensional parameters are listed in table 1. Definitions of
scales and nondimensional parameters are collected in table 2.

Symbol | Definition | Eqn. | Description
---|---|---|---
$U$ | ${B_{0}hI_{0}}/[{2\pi\nu\varrho(r_{i}+r_{o})}]$ | (9) | Shallow-layer velocity scale [m/s]
$U_{\text{deep}}$ | ${B_{0}I_{0}(r_{o}-r_{i})^{2}}/[{8\pi h\nu\varrho(r_{i}+r_{o})}]$ | (30) | Deep-layer velocity scale [m/s]
$T$ | $(r_{i}+r_{o})/U$ | (10) | Advective timescale [s]
$T_{\text{su}}$ | ${4h^{2}}/({\pi^{2}\nu})$ | (36) | Spin-up timescale [s]
$T_{\text{circ}}$ | ${\pi(r_{i}+r_{o})}/{U}$ | (37) | Channel circulation timescale [s]
[2pt/2pt] |  |  |
Fr | $U/\sqrt{gh}$ | (18) | Froude number
$\Rey$ | $h^{2}/(\nu T)={B_{0}I_{0}h^{3}}/[{2\pi\varrho\nu^{2}(r_{i}+r_{o})^{2}}]$ | (12) | Control Reynolds number
Rm | $\nu\mu_{0}\sigma_{0}\Rey$ | (13) | Magnetic Reynolds number
Ha | $B_{0}h\sqrt{\sigma_{0}/(\varrho\nu)}$ | (4) | Hartmann number
$\mathcal{H}$ | $h/r_{o}$ | (14) | Fluid layer aspect ratio
$\mathcal{R}$ | $r_{i}/r_{o}$ | (14) | Sidewall radius ratio
$\mathcal{D}$ | $d/r_{o}$ | (51) | nondimensional obstacle diameter
$\Rey_{d}$ | $Ud/\nu$ | (53) | Obstacle-scale Reynolds number
$\Delta_{i}$ | $\delta_{95,i}/(r_{o}-r_{i})$ | (28) | Relative inner boundary layer thickness
$\Delta_{o}$ | $\delta_{95,o}/(r_{o}-r_{i})$ | (28) | Relative outer boundary layer thickness
Table 2: Scales and nondimensional parameters. All quantities below the dashed
line are nondimensional.

The resulting nondimensional axisymmetric equations can be expressed in scaled
cylindrical coordinates ($\rho$,$\theta$,$\zeta$) as

$\displaystyle\Rey\left\\{\partial_{\tau}\upsilon_{\rho}+(\mathcal{R}+1)\left[(\boldsymbol{\upsilon}_{\bot}\bcdot{\bnabla}_{\bot})\upsilon_{\rho}-\frac{\upsilon_{\theta}^{2}}{\rho}+\partial_{\rho}\Pi\right]\right\\}$
$\displaystyle=\frac{1}{\rho}{\nabla}^{2}_{\bot}(\rho\upsilon_{\rho}),$ (15a)
$\displaystyle\Rey\left\\{\partial_{\tau}\upsilon_{\theta}+(\mathcal{R}+1)\left[(\boldsymbol{\upsilon}_{\bot}\bcdot{\bnabla}_{\bot})\upsilon_{\theta}+\frac{\upsilon_{\rho}\upsilon_{\theta}}{\rho}\right]\right\\}$
$\displaystyle=\frac{1}{\rho}{\nabla}^{2}_{\bot}(\rho\upsilon_{\theta})+\frac{\mathcal{R}+1}{\rho}\kappa(\tau),$
(15b)
$\displaystyle\Rey\left\\{\partial_{\tau}\upsilon_{\zeta}+({\mathcal{R}+1})\left[(\boldsymbol{\upsilon}_{\bot}\bcdot{\bnabla}_{\bot})\upsilon_{\zeta}+\frac{1}{\mathcal{H}^{2}}{\partial_{\zeta}\Pi}\right]\right\\}$
$\displaystyle=\frac{1}{\rho}{\nabla}^{2}_{\bot}(\rho\upsilon_{\zeta})+\mathcal{H}^{2}\frac{\upsilon_{\zeta}}{\rho},$
(15c)
$\displaystyle\partial_{\rho}(\rho\upsilon_{\rho})+{\rho}\partial_{\zeta}\upsilon_{\zeta}$
$\displaystyle=0,$ (15d)

where

$\boldsymbol{\upsilon}_{\bot}=\upsilon_{\rho}\boldsymbol{e_{r}}+\upsilon_{\zeta}\boldsymbol{e_{z}},\qquad{\bnabla}_{\bot}=\boldsymbol{e_{r}}\partial_{\rho}+\boldsymbol{e_{z}}\partial_{\zeta},\quad\text{and}\quad{\nabla}^{2}_{\bot}=\mathcal{H}^{2}\rho\partial_{\rho}\frac{1}{\rho}\partial_{\rho}+\partial_{\zeta}^{2}.$
(16)

We assume vertical hydrostasy, $0={\partial_{\zeta}\Pi}$, and cyclostrophic
balance, ${\upsilon_{\theta}^{2}}/{\rho}=\partial_{\rho}\Pi$, achieved via
deflection of the free surface. Under these conditions, the meridional flow
$\boldsymbol{\upsilon}_{\bot}$ vanishes, leaving a linear equation for
axisymmetric magneto-Stokes flow:

$\Rey\partial_{\tau}(\rho\upsilon_{\theta})={\nabla}^{2}_{\bot}(\rho\upsilon_{\theta})+({\mathcal{R}+1})\kappa(\tau).$
(17)

The ratio of the deflection the free surface $\sim U^{2}/g$ to the depth of
the fluid $h$ scales with the square of the Froude number

$\textit{Fr}={U}/{\sqrt{gh}}.$ (18)

In our validation experiments, $\textit{Fr}\lesssim 10^{-2}$, and the free
surface may be treated as flat (e.g., Greenspan & Howard 1963). Thus, we seek
solutions to (17) on the rectangular domain
$(\rho,\zeta)\in(\mathcal{R},1)\times(0,1)$ with boundary conditions

$\upsilon_{\theta}(\mathcal{R},\zeta,\tau)=\upsilon_{\theta}(1,\zeta,\tau)=\upsilon_{\theta}(\rho,0,\tau)=\left[\partial_{\zeta}\upsilon_{\theta}\right]_{\zeta=1}=0.$
(19)

###  2.2 Spin up from rest

We consider solutions to (17) with boundary conditions (19) that develop from
an initially quiescent fluid ($\upsilon_{\theta}=0$ at $\tau=0$) once constant
electric current is applied ($\kappa=1$ for $\tau>0$). The initial value
problem (IVP) may be solved by splitting $\upsilon_{\theta}$ into stationary
and time-dependent components:
$\upsilon_{\theta}(\rho,\zeta,\tau)=\bar{\upsilon}_{\theta}(\rho,\zeta)+\upsilon_{\theta}^{\prime}(\rho,\zeta,\tau)$,
where

${\nabla}^{2}_{\bot}(\rho\bar{\upsilon}_{\theta})=-({\mathcal{R}+1}),$ (20)

and

$\Rey\partial_{\tau}(\rho\upsilon_{\theta}^{\prime})={\nabla}^{2}_{\bot}(\rho\upsilon_{\theta}^{\prime}),\quad\upsilon_{\theta}^{\prime}=-\bar{\upsilon}_{\theta}\;\text{at}\;\tau=0,$
(21)

each under the same boundary conditions as those on $\upsilon_{\theta}$.

####  2.2.1 Stationary solution

We proceed by finding the stationary solution, which obeys (20) and satisfies
$\bar{\upsilon}_{\theta}(\mathcal{R},\zeta)=\bar{\upsilon}_{\theta}(1,\zeta)=\bar{\upsilon}_{\theta}(\rho,0)=\left[\partial_{\zeta}\bar{\upsilon}_{\theta}\right]_{\zeta=1}=0$.
The solution may be expanded using the basis of vertical eigenfunctions:

$\bar{\upsilon}_{\theta}(\rho,\zeta)=\sum_{n=1}^{\infty}u_{n}(\rho)\sin(k_{n}\zeta),\quad
k_{n}=\pi(n-1/2),$ (22)

where $k_{n}$ is the $n$th wavenumber in the vertical direction. Expanding
both sides of (20) with the same basis and making the transformation
$\rho=s\mathcal{H}/k_{n}$,
$u_{n}(\rho)={2(\mathcal{R}+1)}{k_{n}^{-2}\mathcal{H}^{-1}}v_{n}(s)$ yields

$s^{2}\frac{\mathrm{d}^{2}v_{n}}{\mathrm{d}s^{2}}+s\frac{\mathrm{d}v_{n}}{\mathrm{d}s}-(s^{2}+1)v_{n}(s)=-s.$
(23)

The ordinary differential equation (23) has solutions of the form
$v_{n}(s)=s^{-1}-A_{n}\rmI_{1}(s)-B_{n}\rmK_{1}(s)$, where $\rmI_{1}(s)$ and
$\rmK_{1}(s)$ are the modified Bessel functions of the first and second kind,
respectively. Constants $A_{n}$ and $B_{n}$ are determined from the radial
boundary conditions $u_{n}(\mathcal{R})=v_{n}(\mathcal{R}k_{n}/\mathcal{H})=0$
and $u_{n}(1)=v_{n}(k_{n}/\mathcal{H})=0$ as

$A_{n}=\frac{\mathcal{H}}{k_{n}\mathcal{R}}\left[\frac{\mathcal{R}\rmK_{1}\left({k_{n}\mathcal{R}}/{\mathcal{H}}\right)-\rmK_{1}\left({k_{n}}/{\mathcal{H}}\right)}{\rmI_{1}\left({k_{n}}/{\mathcal{H}}\right)\rmK_{1}\left({k_{n}\mathcal{R}}/{\mathcal{H}}\right)-\rmK_{1}\left({k_{n}}/{\mathcal{H}\\\
}\right)\rmI_{1}\left({k_{n}\mathcal{R}}/{\mathcal{H}}\right)}\right],$ (24a)

$B_{n}=\frac{\rmI_{1}\left({k_{n}}/{\mathcal{H}}\right)-\mathcal{R}\rmI_{1}\left({k_{n}\mathcal{R}}/{\mathcal{H}}\right)}{\mathcal{R}\rmK_{1}\left({k_{n}\mathcal{R}}/{\mathcal{H}}\right)-\rmK_{1}\left({k_{n}}/{\mathcal{H}}\right)}A_{n}.$
(24b)

Finally, changing back to original (nondimensional) coordinates $(\rho,\zeta)$
and substituting back into (22) yields

$\bar{\upsilon}_{\theta}(\rho,\zeta)=\left(\frac{\mathcal{R}+1}{2}\right)\frac{2\zeta-\zeta^{2}}{\rho}-\sum_{n=1}^{\infty}\frac{2(\mathcal{R}+1)}{k_{n}^{2}\mathcal{H}}\left[A_{n}\rmI_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)+B_{n}\rmK_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)\right]\sin(k_{n}\zeta).$
(25)

The first term in (25) is equal to the asymptotic solution of (20) in the
shallow-layer limit $\mathcal{H}\to 0$:

$\lim_{\mathcal{H}\to
0}\bar{\upsilon}_{\theta}=\left(\frac{\mathcal{R}+1}{2}\right)\frac{2\zeta-\zeta^{2}}{\rho},$
(26)

which contains the quadratic dependence on height typical of creeping channel
flows (e.g., Milne-Thomson 1938) and inherits the inverse dependence on radius
of the Lorentz force (since $\lVert\boldsymbol{J}\rVert\propto 1/r$). Note
that we have chosen our velocity scale $U$ such that $\bar{\upsilon}_{\theta}$
has unit value at the surface mid-gap ($\zeta=1$, $\rho=\mathcal{R}/2+1/2$).

For $\mathcal{H}>0$, the series in (25) produces boundary layers at both
sidewalls. The size of these boundary layers may be described by their 95%
thicknesses $\delta_{95,i}$, $\delta_{95,o}$, which are defined such that

$\bar{\upsilon}_{\theta}=0.95\lim_{\mathcal{H}\to
0}\bar{\upsilon}_{\theta}\quad\text{at}\quad\rho=(r_{i}+\delta_{95,i})/r_{o},\;(r_{o}-\delta_{95,o})/r_{o}.$
(27)

Figure 2(a) shows the stationary solution given by (25) for a channel with
geometric ratios $\mathcal{H}=0.05$ and $\mathcal{R}=0.25$. Contours
correspond to different vertical positions $\zeta$ within the fluid, and
dashed yellow lines indicate the 95% thicknesses of the sidewall boundary
layers. Figure 2(b) plots the terms of (20) computed from the stationary
solution (25) at the free surface. Basal viscous drag
($\partial_{\zeta}^{2}\bar{\upsilon}_{\theta}$) and the Lorentz force dominate
the momentum balance in the bulk, and sidewall viscous drag
($\mathcal{H}^{2}\rho\partial_{\rho}[{\rho}^{-1}\partial_{\rho}\bar{\upsilon}_{\theta}]$)
is only significant in the boundary layers.

Figure 2: (a) Stationary magneto-Stokes flow solution given by (25) for a
channel with geometric ratios $\mathcal{H}=0.05$ and $\mathcal{R}=0.25$.
Labelled contours trace the solution profile at different heights $\zeta$
above the channel base. A dash-dotted grey line shows the solution in the
shallow limit ($\mathcal{H}\to 0$), and yellow dashed lines correspond to the
95% thicknesses of the sidewall boundary layers. (b) Terms in the azimuthal
momentum balance (20) evaluated at $\zeta=1$, computed from the same solution.
Lorentz, basal viscous drag, and sidewall viscous drag forces are computed as
$(\mathcal{R}+1)/\rho$, $\partial_{\zeta}^{2}\bar{\upsilon}_{\theta}$, and
$\mathcal{H}^{2}\rho\partial_{\rho}({\rho}^{-1}\partial_{\rho}\bar{\upsilon}_{\theta})$,
respectively.

The size of each boundary layer relative to the channel width is captured by

$\Delta_{i}=\delta_{95,i}/(r_{o}-r_{i})\quad\text{and}\quad\Delta_{o}=\delta_{95,o}/(r_{o}-r_{i}),$
(28)

which scale with the geometric ratios $\mathcal{H}$, $\mathcal{R}$ as

$\Delta_{i}\approx\Delta_{o}\approx
10(0.22+\mathcal{R}^{2.5})\mathcal{H}^{1.0}$ (29)

for $7.5\times 10^{-3}\leq\mathcal{H}\leq 0.06$ and $0.01\leq\mathcal{R}\leq
0.7$. All numerical factors and powers in (29) are fit to values of
$\Delta_{o}$ computed from (25). Since $\mathcal{R}$ is constrained to lie
between zero and unity, the dominant behaviour of (29) is a linear dependence
on the aspect ratio, $\mathcal{H}$. The boundary layers shrink in proportion
to $\mathcal{H}$ as basal viscous drag dominates over sidewall viscous
effects.

If the layer depth is increased while keeping inertial forces small ($\Rey\ll
1$), the forces in the bulk become dominated by the balance of Lorentz and
sidewall viscous drag forces; the basal viscous drag becomes subdominant. As
$\mathcal{H}\to\infty$ and $\Rey\to 0$, the sidewall boundary layers merge
($\Delta_{i}+\Delta_{o}>1$) and the flow becomes independent of depth apart
from a thin viscous boundary layer above the channel base. Thus, the condition
$\Delta_{i}+\Delta_{o}=1$ traces a curve in $\mathcal{R}$,$\mathcal{H}$-space
below which magneto-Stokes flow is considered “shallow” and above which it is
considered “deep” (assuming $\Rey\ll 1$). Figure 3(a) uses the scaling (29) to
separate these regimes in the space of aspect and radius ratios. Blue, gray,
and purple dots in figure 3(a) correspond to shallow, transitional, and deep
flows in channels of the same radius ratio $\mathcal{R}=0.25$, whose predicted
radial profiles are plotted in figure 3(b). Each of the plotted velocity
profiles is scaled by its magnitude at the surface mid-gap,
$u_{\theta,\text{mid-gap}}=U\bar{\upsilon}(\rho=\mathcal{R}/2+1/2,\zeta=1)$.
The velocity scale $U$ is not used to normalise the profiles in figure 3(b)
because it is only appropriate for the shallow regime, in which
$u_{\theta,\text{mid-gap}}$ grows with increasing aspect ratio
$\mathcal{H}\equiv h/r_{o}$ as basal viscous drag ($\sim 2{\nu U}/{h^{2}}$)
becomes weaker relative to the Lorentz force ($\sim B_{0}I_{0}/[\pi
h\varrho(r_{i}+r_{o})]$). In the deep regime, the opposite is true:
$u_{\theta,\text{mid-gap}}$ shrinks with increasing $\mathcal{H}$ as the
Lorentz force ($\sim B_{0}I_{0}/[\pi h\varrho(r_{i}+r_{o})]$) becomes weaker
relative to the sidewall viscous drag ($\sim 8\nu
U_{\text{deep}}/[r_{o}-r_{i}]^{2}$). This deep-layer behaviour is captured by
the velocity scale

$U_{\text{deep}}=\frac{B_{0}I_{0}(r_{o}-r_{i})^{2}}{8\pi
h\nu\varrho(r_{i}+r_{o})}=\left(\frac{1-\mathcal{R}}{2\mathcal{H}}\right)^{2}U.$
(30)

Figure 3: (a) Regime diagram for annular magneto-Stokes flow in the space of
radius ratios $\mathcal{R}$ and aspect ratios $\mathcal{H}$. Renderings of
cylindrical annuli beside the plot depict (exaggerated) channel geometries
implied by axes values. Gradated blue, white, and purple tones correspond to
increasing values of $\Delta_{i}+\Delta_{o}$ predicted by (29). A dashed grey
line indicates the boundary $\Delta_{i}+\Delta_{o}=1$ between shallow and deep
layer magneto-Stokes regimes. Points labelled with roman numerals correspond
to laboratory cases discussed in Sections 3 and 4. (b) A comparison of
magneto-Stokes flows in channels of varying aspect ratio $\mathcal{H}$ at a
fixed radius ratio $\mathcal{R}=0.25$. Plotted are radial profiles of surface
azimuthal velocity predicted from theory (25) and scaled by surface values at
the channel centre, $u_{\theta,\text{mid-
gap}}=U\bar{\upsilon}(\rho=\mathcal{R}/2+1/2,\zeta=1)$. Solid blue, dashed
grey, and solid purple curves correspond to points of the same colour in the
regime diagram. Dash-dotted blue and purple curves show the corresponding
shallow (26) and deep (33) solutions, respectively.

Though the deep limit ($\mathcal{H}\to\infty$) is not the focus of the present
study, obtaining an asymptotic steady-state solution is straightforward after
re-scaling velocity as $w_{\theta}=u_{\theta}/U_{\text{deep}}$. Then, the
stationary component of $w_{\theta}$ is
$\bar{w}_{\theta}={4\mathcal{H}^{2}}/{(1-\mathcal{R})^{2}}\bar{\upsilon}_{\theta}$,
and (20) may be rewritten as

$\rho\partial_{\rho}\frac{1}{\rho}\partial_{\rho}(\rho\bar{w}_{\theta})+\frac{1}{\mathcal{H}^{2}}\partial_{\zeta}^{2}(\rho\bar{w}_{\theta})=-\frac{4(1+\mathcal{R})}{(1-\mathcal{R})^{2}}.$
(31)

As $\mathcal{H}\to\infty$,

$\rho\frac{\mathrm{d}}{\mathrm{d}\rho}\left[\frac{1}{\rho}\frac{\mathrm{d}}{\mathrm{d}\rho}(\rho\bar{w}_{\theta})\right]=-\frac{4(1+\mathcal{R})}{(1-\mathcal{R})^{2}}.$
(32)

Then, the steady solution for an infinitely deep channel is

$\lim_{\mathcal{H}\to\infty}U\bar{\upsilon}_{\theta}/U_{\text{deep}}=\bar{w}_{\theta}(\rho)=\frac{2\left(1-\rho^{2}\right)\mathcal{R}^{2}\log(\mathcal{R})-2\left(1-\mathcal{R}^{2}\right)\rho^{2}\log\left(\rho\right)}{(1-\mathcal{R})^{3}\rho},$
(33)

which is plotted in figure 3(b) as a purple dash-dotted curve. For the
dimensional form of (33), see Gleeson et al. (2004). The maximum
nondimensional velocity is of order unity: $\max_{\rho}w_{\theta}\in(0.74,1)$
for $\mathcal{R}\in(0,1)$, justifying the choice of velocity scale,
$U_{\text{deep}}$.

####  2.2.2 Time-dependent component

The initial value problem (21) may be solved to arbitrary accuracy via
Fourier-Bessel series, as is done in Appendix A. Here, we provide an
approximate solution, which we will show suffices for predicting shallow-layer
laboratory flows. We retain the influence of sidewall viscous forces on the
stationary solution, but argue that, in a shallow layer, they have a
negligible effect on how fast the fluid spins up from rest once voltage is
applied to the electrodes. Under this assumption, we take the shallow-layer
limit $\mathcal{H}\to 0$ in (21), yielding

$\Rey\partial_{\tau}\upsilon_{\theta}^{\prime}\approx\partial_{\zeta}^{2}\upsilon_{\theta}^{\prime},\quad\upsilon_{\theta}^{\prime}=-\bar{\upsilon}_{\theta}\;\text{at}\;\tau=0,$
(34)

so that

$\upsilon_{\theta}^{\prime}(\rho,\zeta,\tau)\approx-\sum_{n=1}^{\infty}u_{n}(\rho)\sin(k_{n}\zeta)\exp\left(-\frac{k_{n}^{2}}{\Rey}\tau\right).$
(35)

The gravest mode ($n=1$) spins up slowest at a timescale

$T_{\text{su}}=\frac{4\Rey}{\pi^{2}}T=\frac{4h^{2}}{\pi^{2}\nu}.$ (36)

For comparison, we may define a circulation timescale equal to the time it
takes for a fluid parcel at the surface mid-gap ($\zeta=1$,
$\rho=\mathcal{R}/2+1/2$) to circulate through the channel once:

$T_{\text{circ}}=\frac{2\pi(r_{i}+r_{o})/2}{U}=\pi T.$ (37)

The ratio of timescales is proportional to the a priori Reynolds number:

$\Rey=\frac{\pi^{3}}{4}\frac{T_{\text{su}}}{T_{\text{circ}}}.$ (38)

Shallow layers with $\Rey\to 0$ spin-up instantaneously ($T_{\text{su}}\ll
T_{\text{circ}}$), as expected for Stokes flow. For $\Rey>0$, there is a time-
dependence in the flow following any changes in the electromagnetic forcing
(e.g., turning the current on or off).

For simplicity, we retain only the time-dependent behaviour of the slowest
mode, such that

$\upsilon_{\theta}^{\prime}(\rho,\zeta,\tau)\approx-\left\\{\sum_{n=1}^{\infty}u_{n}(\rho)\sin(k_{n}\zeta)\right\\}\exp\left(-\frac{\pi^{2}}{4\Rey}\tau\right).$
(39)

The approximate solution for a flow spun up from rest by a constant
electromagnetic forcing is formed as the sum of the stationary (25) and time-
dependent (39) solutions:

$\displaystyle\upsilon_{\theta}(\rho,\zeta,\tau)=\bar{\upsilon}_{\theta}+\upsilon_{\theta}^{\prime}$
$\displaystyle\approx$
$\displaystyle\left\\{\astrut\sum_{n=1}^{\infty}\frac{2(\mathcal{R}+1)}{k_{n}^{2}\mathcal{H}}\left[\frac{\mathcal{H}}{k_{n}\rho}-A_{n}\rmI_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)-B_{n}\rmK_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)\right]\sin(k_{n}\zeta)\astrut\right\\}$
(40)
$\displaystyle\mbox{}\times\left[1-\exp\left(-\frac{\pi^{2}}{4\Rey}\tau\right)\right],$

where $k_{n}$, $A_{n}$, and $B_{n}$ are defined in (22, 24). This solution is
approximate only in its time-dependent behaviour, and approaches the exact
solution as $\tau/\Rey\to\infty$. In the following sections, we use laboratory
experiments and numerical solution of the IVP to show that this error is as
small as the errors in the model assumptions, making (40) a sufficient
description of our actual shallow-layer magneto-Stokes flows.

##  3 Experimental methods

We validate the approximate solution (40) via four laboratory experiments
(Cases I–IV), matched with numerical solution of the IVP described in Section
2.2. The results of these validation cases are discussed in Section 4.
Additional laboratory experiments (Cases KR, HS-I–III) which demonstrate
kinematic reversibility and a connection to Hele-Shaw obstacle flows, are
detailed in Section 5.

###  3.1 Laboratory experiments

Validation experiments (Cases I–IV) are performed using an open-top annular
channel consisting of a 17.5 cm-radius steel outer cylindrical sidewall, clear
acrylic base, and a removable stainless steel inner cylindrical sidewall,
which may be replaced with cylinders of different radii. The channel is placed
inside the solenoidal electromagnet bore of UCLA’s RoMag device (Xu et al.,
2022), and a DC bench power supply provides a controlled current $I$ between
the channel sidewalls. An $80.0\pm 0.05$ g/L NaCl:H2O solution is used as the
working fluid for all cases; 0.1 mL of dish detergent is added for every litre
of solution to reduce surface tension. The fluid is estimated to have
electrical conductivity $\sigma=12.3\pm 0.1$ S/m, kinematic viscosity
$\nu=(1.10\pm 0.05)\times 10^{-6}$ m2/s, and density $\varrho=1059.1\pm 0.7$
kg/m3, using the salinity-based models of Park (1964), Isdale et al. (1972),
and Isdale & Morris (1972), respectively.

| Case
---
|  $I_{0}$ (a)
---
($\pm 0.003$ A)
|  $B_{0}$ (Gauss)
---
|  $h$ (cm)
---
($\pm 0.05$ cm)
|  $r_{i}$ (cm)
---
($\pm 0.03$ cm)
|  $r_{o}$ (cm)
---
($\pm 0.03$ cm)
|  $d$ (cm)
---
($\pm 0.03$ cm)
|  $U$ (cm/s)
---
I | 0.081 | $300\pm 2\hphantom{0}$ | 0.4 |  04.44 | 17.52 | — | $0.56\pm 0.08$
II | 0.099 | $300\pm 6\hphantom{0}$ | 0.4 |  07.62 | 17.52 | — | $0.63\pm 0.09$
III | 0.039 | $200\pm 4\hphantom{0}$ | 0.8 |  07.62 | 17.52 | — | $0.36\pm 0.04$
IV | 0.099 | $350\pm 10$ | 0.4 | 10.14 | 17.52 | — | $0.7\pm 0.1$
KR | 0.099 | $350\pm 10$ | 0.4 | 10.14 | 17.52 | — | $0.7\pm 0.1$
HS-I | 0.060 | $190\pm 10$ | 0.2 |  03.75 |  09.84 | 1.59 | $0.24\pm 0.06$
HS-II | 0.801 | $190\pm 10$ | 0.4 |  03.75 |  09.84 | 1.59 | $6.3\pm 0.9$
HS-III | 1.011 | $190\pm 10$ | 0.4 |  04.44 | 17.52 | 4.41 | $4.9\pm 0.7$
Table 3: Dimensional experimental parameters and predicted velocity $U$ at
surface mid-gap, computed using (9). Error in values of $U$ reflect the
propagation of measurement uncertainty of the control parameters. Obstacle
diameter $d$ is only relevant to Cases HS-I–III detailed in Section 5.2.

Inner radius, fluid height, electric current, and magnetic field strength are
varied across Cases I–IV; values of these dimensional parameters are reported
in Table 3. Cases I–IV span three different radius ratios $\mathcal{R}$ =
0.25, 0.44, 0.58 and two different aspect ratios $\mathcal{H}$ = 0.023, 0.046.
Figure 3(a) shows where each of the four cases (blue, orange, red, and green
points) lies in the shallow-layer magneto-Stokes region of
$\mathcal{H},\mathcal{R}$-space, discussed in Section 2.2.1. Fluid depth is
kept at or above 0.2 cm to minimise relative differences in depth due to
capillary effects. Voltage across the electrodes is maintained under $\sim$1.5
V, which prevents formation of unwanted bubbles due to electrolysis of water.
Under this constraint, electric current and magnetic field strength are held
between 0.04–0.1 A and 200–350 Gauss, respectively, to keep $\Rey<1$ and
$\textit{Fr}^{2}\ll 1$. The low electrical conductivity of the working fluid
ensures that $\textit{Ha}\lesssim 10^{-2}$ and $\textit{Rm}\lesssim 10^{-10}$.
Values of these nondimensional control parameters are reported in Table 4.

Case | $\mathcal{R}$ | $\mathcal{H}$ | $\Rey$ | Fr | Ha | Rm | $\mathcal{D}$ | $\Rey_{d}$
---|---|---|---|---|---|---|---|---
I | 0.25 | 0.023 | 0.380 | 0.030 | 0.012 | $6.6\times 10^{-12}$ | — | —
II | 0.44 | 0.023 | 0.360 | 0.032 | 0.012 | $6.3\times 10^{-12}$ | — | —
III | 0.44 | 0.046 | 0.780 | 0.012 | 0.016 | $1.3\times 10^{-11}$ | — | —
IV | 0.58 | 0.023 | 0.350 | 0.034 | 0.014 | $6.1\times 10^{-12}$ | — | —
KR | 0.58 | 0.023 | 0.350 | 0.034 | 0.014 | $6.1\times 10^{-12}$ | — | —
HS-I | 0.38 | 0.020 | 0.058 | 0.016 | $3.8\times 10^{-3}$ | $1.0\times 10^{-12}$ | 0.16 | 32
HS-II | 0.38 | 0.041 | 6.200 | 0.300 | $7.6\times 10^{-3}$ | $1.1\times 10^{-10}$ | 0.16 | 850
HS-III | 0.25 | 0.023 | 3.000 | 0.230 | $7.6\times 10^{-3}$ | $5.2\times 10^{-11}$ | 0.25 | 1800
Table 4: Nondimensional experimental parameters. The radius ratio
$\mathcal{R}$, aspect ratio $\mathcal{H}$, control Reynolds number $\Rey$,
Froude number Fr, Hartmann number Ha, and magnetic Reynolds number Rm are
defined in Section 2. Nondimensional obstacle diameter $\mathcal{D}$ and
obstacle Reynolds number $\Rey_{d}$ are only relevant to Cases HS-I–III and
are defined in Section 5.

The procedure for Cases I–IV is as follows: before the start of each case, a
streak of buoyant blue dye is drawn on the surface of the quiescent fluid,
connecting inner and outer sidewalls along a diameter of the annulus. At time
$t=0$, voltage is applied across the sidewalls such that a constant current is
maintained. From time $t=0$ to $t=0.5T_{\text{circ}}$, an overhead camera
records the advection of the buoyant dye streak over the fluid surface. The
emplacement of the dye streak at the beginning of each case often leaves
behind small spots of dye on the channel’s base. These spots are digitally
masked out of each video so that the surface dye streak can be uniquely
identified using basic image processing methods.

Video recordings are processed with a custom phase-boundary-tracking program
to isolate and track the leading edge of the advected dye streak. Each video
is first corrected for minor perspective distortion using inner and outer
sidewall boundaries as markers. The sidewalls are identified via Canny edge
detection and Hough circle transform (Canny, 1986; Yuen et al., 1990),
implemented in the OpenCV Python library (Bradski, 2000). Each transformed
video frame is then thresholded based on luminance and colour to isolate the
bright blue streak of dye, the edges of which are identified via Canny edge
detection with morphological closing. Uncertainty in position of the dye
streak edges is computed as the change in estimated position under a 10%
adjustment of color threshold values. The leading and trailing edges of the
advected dye streak are identified in each frame by comparing angular
positions, accounting for branch cuts in azimuth.

The (Lagrangian) angular position $\theta(r,t)$ of the leading dye streak edge
may then be used to compute radial profiles of azimuthal velocity over time. A
profile of $\theta$ over $r$ is constructed at every 0.5 $T_{\text{su}}$ from
$t=0$ to $t=0.5T_{\text{circ}}$. This provides enough temporal resolution of
the spin-up, while ensuring that the angular distance travelled by the dye
streak over each time step is greater than the uncertainty in dye streak
position by at least a factor of 3. At each time, the angular position
$\theta$ of the leading edge is determined at 15 points across the channel (in
$r$), excluding the menisci at the inner and outer sidewalls where dye appears
to spread rapidly via adhesion instead of fluid advection. A time-series of
azimuthal velocity $u_{\theta}(r)$ is then estimated via second-order central
difference of $\theta(r,t)$ over time, and uncertainty in the dye streak
position is propagated to error in $u_{\theta}(r)$.

###  3.2 Direct numerical simulations

Cases I–IV are matched with direct numerical simulations (DNS) of flow
governed by (17) with conditions (19), $\upsilon_{\theta}=0$ at $\tau=0$, and
$\kappa=1$ for $\tau>0$. We employ the Dedalus pseudospectral framework (Burns
et al., 2020), using 201 Chebyshev modes in $\rho$ and 100 cosine modes in
$\zeta$. The no-slip condition is enforced at $\zeta=0$ using the volume-
penalty method (Hester et al., 2021), which adds a spatially-masked linear
damping term $-\Gamma({\zeta}/{\delta}){\upsilon_{\theta}}/{\tau_{\text{VP}}}$
to the momentum equation (17):

$\partial_{\tau}\upsilon_{\theta}-\frac{1}{\Rey}\left(\mathcal{H}^{2}\partial_{\rho}^{2}\upsilon_{\theta}+\partial_{\zeta}^{2}\upsilon_{\theta}\right)+\frac{\upsilon_{\theta}}{\tau_{\text{VP}}}=\frac{\mathcal{H}^{2}}{\Rey\rho}\left(\partial_{\rho}\upsilon_{\theta}-\frac{\upsilon_{\theta}}{\rho}\right)+\frac{\mathcal{R}+1}{\Rey\rho}+\left[1-\Gamma\left(\frac{\zeta}{\delta}\right)\right]\frac{\upsilon_{\theta}}{\tau_{\text{VP}}},$
(41)

where $\Gamma(x)=[1-\rm{erf}(\sqrt{\pi}x)]/2$ is a masking function and
$\tau_{\text{VP}}$ is the volume-penalty damping timescale nondimensionalised
by $T$. The masking function is smoothed over a vertical length scale $\delta$
(set here to $\delta=0.05$), which is used to determine an appropriate value
for $\tau_{\text{VP}}$. This is effected by requiring the damping length-scale
$\sqrt{\tau_{\text{VP}}/\Rey}$ to be proportional to the smoothing scale:
$\sqrt{\tau_{\text{VP}}/\Rey}=\delta/\delta^{*}$. The proportionality constant
$\delta^{*}$ is set to the optimal value found by Hester et al. (2021),
$\delta^{*}=3.11346786$; this choice of $\delta^{*}$ eliminates the
displacement length error associated with the mask, $\Gamma(x)$. See Hester et
al. (2021) for further details on optimising the volume-penalty method.

In all simulations, the system is integrated from $\tau=0$ to
$\tau=1.1T_{\text{circ}}/T=1.1\pi$ over $2\times 10^{4}$ time stpdf using the
second-order semi-implicit backwards difference (SBDF2) scheme (Wang & Ruuth
2008, Eqn. 2.8). The terms on the left-hand side of (41) are time-integrated
implicitly while those on the right-hand side are treated explicitly.

Reported results match those obtained at half their spatial resolution as well
as the analytical solution at $t=3T_{\text{su}}$. These computations were
inexpensive and thus preferable to computing the full analytic solution of the
IVP (see Appendix A) for each case. The code used for these simulations and
for the dye-tracking described in Section 3.1 is available online
(https://github.com/cysdavid/magnetoStokes).

##  4 Results

Figure 4 shows three photographs of laboratory Case I, taken when power is
turned on (a), after 0.25 circulation times (b), and after 0.5 circulation
times (c). The time-integrated analytical solution (40) and DNS results are
overlain in magenta and grey, respectively. Laboratory, analytical, and
numerical results match well in the bulk, differing most within the boundary
layers (dashed yellow lines). A closeup (b, inset) of the inner boundary layer
shows the laboratory dye streak trailing behind the analytical and numerical
solutions, resulting from different rates of spin-up within the boundary layer
rather than different final velocities. The bulk flow also exhibits finite
spin-up time effects. In each panel of figure 4, a magenta dot is placed on
the analytical curve at $r=(r_{i}+r_{o})/2$. For flow that has fully spun up,
$T_{\text{circ}}$ is the time it takes for this dot to make one revolution. In
figure 4(b), the magenta dot has travelled slightly less than a quarter
revolution from $t=0$ to $t=0.25T_{\text{circ}}$, a product of the finite
$\Rey$ value in our experiments.

Figure 4: Snapshots of a free surface dye track (blue) from laboratory Case I,
(a) when power is turned on, (b) after $\sim$5 spin up times, and (c) after
$\sim$10 spin up times. Time integrations of the approximate analytical
solution (40) and DNS result are overlain as magenta and grey curves,
respectively. Dotted yellow circles correspond to the 95% thickness of each
sidewall boundary layer as predicted from (40). The red cord near the 12
o’clock position in each photograph is the electrical wire leading from the
power supply to the inner electrode.

Figure 5: Scaled azimuthal velocity at $r=(r_{i}+r_{o})/2$, $z=h$ vs. time for
Cases I–IV. For each case, time is nondimensionalised by the spin-up
timescale, $T_{\text{su}}$. Theoretical curves are computed using (40). Error
bars and shaded envelopes represent $\pm 1$ standard deviation propagated from
uncertainty in the dye tracking velocimetry algorithm and from measurement of
input parameters.

Figure 6: Radial profiles of scaled azimuthal velocity at the free surface for
Cases I-IV at (a) 1 spin-up time after rest, (b) 2 spin-up times, and (c) 3
spin-up times. Solid theoretical curves are computed using (40). DNS are shown
via dashed curves. Also plotted is the $1/\rho$ profile (dash-dotted grey
curve) corresponding to the shallow ($\mathcal{H}\to 0$), long-time solution
(26) at $\zeta=1$. Error bars and shaded envelopes represent $\pm 1$ standard
deviation propagated from uncertainty in the dye tracking velocimetry
algorithm and from measurement of controlled parameters. Rendered cylindrical
annuli in the lower legend depict the (exaggerated) channel geometry of each
case.

Figure 5 shows the time evolution of scaled azimuthal velocity at the surface
mid-gap ($\rho=\mathcal{R}/2+1/2$, $\zeta=1$) for cases I–IV. Included are
laboratory data (points), DNS results (dashed curves), and analytical
predictions (solid curves) given by (40). Bars on the data points correspond
to error introduced by the dye-tracking algorithm and from measurement
uncertainty propagated through the computed velocity scale $U$. The latter
also introduces uncertainty to the analytical and numerical solutions, which
are scaled by $U$. Shaded regions show the $\pm 1$ standard deviation envelope
around the theoretical predictions. Time is nondimensionalised by the scale
relevant to spin-up, $T_{\text{su}}$. Data, DNS, and theory collapse onto the
same curve given by (40), that tends toward $u_{\theta}/U=1$ as
$t/T_{\text{su}}\to\infty$.

Figure 6 plots radial profiles of scaled azimuthal velocity for all four
cases. Included are laboratory data (points), DNS results (dashed curves), and
the approximate analytical solution given by (40) (solid curves) after 1 spin-
up time (a), 2 spin-up times (b), and 3 spin-up times (c). The scaled velocity
profiles evolve towards the $1/\rho$ curve (gray dash-dotted line) over time,
apart from the boundary layers. A closeup of the inner boundary layer for the
two cases with $\mathcal{R}=0.44$ (panel c, inset) shows a clear separation
between the profiles for $\mathcal{H}=0.023$ (Case II) and $\mathcal{H}=0.046$
(Case III) and excellent agreement with theory.

A slight separation between data, DNS, and theory at Case I’s peak in velocity
(near $\rho=0.3$) for $t/T_{\text{su}}\leq 2$ can be seen in figure 6(a,b).
The DNS result (dashed blue line) spins up faster than the approximate
solution (solid blue line), while the laboratory flow (blue points) spins up
slower. The former is a result of the sidewall viscous drag’s effect on spin-
up. As evident in the full solution in Appendix A, features with higher
spatial frequency in the radial direction (e.g., the sharp peak near
$\rho=0.3$) spin up faster than lower frequency features. This effect is
neglected in the approximate solution (40) but is retained in the DNS.

The mismatch between laboratory data and the approximate solution near
$\rho=0.3$ in figure 6(a,b) may be an effect of nonlinear advection, which we
neglect in our model by assuming that the free surface instantaneously adjusts
to the azimuthal flow at all times. For the finite Reynolds number flow in
Case I ($\Rey=0.38$), it is possible that the dynamic adjustment of the free-
surface during spin-up induces a transient meridional circulation,
$\boldsymbol{\upsilon}_{\bot}=\upsilon_{\rho}\boldsymbol{e_{r}}+\upsilon_{\zeta}\boldsymbol{e_{z}}$.
A different laboratory case (Case KR) discussed in Section 5 exhibits slight
outwards radial motion of surface floaters, supporting the existence of a weak
meridional circulation at finite Reynolds numbers. This secondary flow could
advect momentum away from the peak in azimuthal velocity, $\upsilon_{\theta}$,
during spin-up. The free surface may then equilibrate such that the radial
pressure gradient balances centrifugation, shutting off the meridional
circulation and permitting the azimuthal flow to reach the steady-state
profile (25) predicted by linear theory. This hypothesised mechanism is likely
only significant for channels with small inner cylinders ($\mathcal{R}\lesssim
0.25$), for which the peak surface velocity may be more than twice the value
at the mid-gap. The role of surface tension near the meniscus at the inner
cylinder may further complicate the transient flow behaviour.

##  5 Reproducing properties of classical Stokes flow

For a shallow layer ($\mathcal{H}\ll 1$) under weak electromagnetic forcing
($\Rey\ll 1$), annular magneto-Stokes flow recovers the remarkable, classical
properties of viscous flow—with the precise control afforded by the MHD
forcing. Using the shallow magneto-Stokes system, we demonstrate the kinematic
reversibility result of Taylor (1967) and the analogy to 2D inviscid flows
made by Hele-Shaw (1898).

###  5.1 Kinematic reversibility

Figure 7: Snapshots from laboratory demonstrations of kinematic reversibility
in identical flows with $\Rey=0.35$, visualised with dye (top row) and with
triangular surface floaters (bottom row) (a) when power is turned on; (b)
after driving flow forward for 0.5 channel circulation times, switching off
power, and allowing flow to decelerate over 5 spin-up times; and (c) after
driving flow in the reverse direction for 0.5 channel circulation times,
switching off power, and allowing flow to decelerate over 5 spin-up times.
Gray shaded regions in the top row and pink triangles in the bottom row
indicate the initial locations of dye streaks and floaters.

In the aforementioned limits ($\mathcal{H}\to 0$, $\Rey\to 0$), the sidewall
boundary layers vanish along with the time-evolution of the flow. The surface
velocity profile is then identical to Taylor-Couette flow with inner and outer
sidewall rotation rates given by:

$\Omega_{i}=\frac{B_{0}I(t)h}{4\pi\varrho\nu
r_{i}^{2}}\quad\text{and}\quad\Omega_{o}=\frac{B_{0}I(t)h}{4\pi\varrho\nu
r_{o}^{2}}.$ (42)

On the timescale of circulatory motion $T_{\text{circ}}$ through the annular
channel, vertical viscous forces are transmitted instantaneously and balance
the Lorentz force at all times. Thus, the flow loses its dependence on the
initial conditions and retains only a parametric dependence on time through
the applied electric current, $I(t)$. A reversal of current polarity is
analogous to a reversal of sidewall rotation in the Taylor-Couette cell. Both
induce a reversal of fluid deformation, but the precise control over $I(t)$
coupled with the shallow geometry of the magneto-Stokes system permits
macroscale reversibility in low-viscosity fluids, like saltwater. In contrast,
bench-top Taylor-Couette cells typically require viscous syrup for feasible
low Reynolds flows (Taylor, 1967).

Figure 7 shows laboratory images for two quasi-reversible magneto-Stokes flows
with $\Rey=0.35$ and identical geometry, $\mathcal{H}$, $\mathcal{R}$ (Case
KR). The top row depicts a dye-visualised flow, while the bottom row shows a
flow with triangular floaters made from cut paper. In both experiments, fluid
is spun up, allowed to circulate, and then spun down. Then, current polarity
is switched, and this sequence is performed in reverse.

Panel (a) shows the fluid surface at rest just as current is switched on to
$+0.1$ A. The flow is driven over $0.5T_{\text{circ}}$, and then current is
switched off. To account for finite-Reynolds number effects, the fluid is
allowed to decelerate over $5T_{\text{su}}$, shown in panel (b). A reverse
current, $-0.1$ A, is then applied, fluid is driven over $0.5T_{\text{circ}}$
(towards initial material positions), and then power is switched off. The
right panel (c) shows snapshots at $5T_{\text{su}}$ after current is turned
off. Gray shaded regions in the top row of images indicate the initial
locations of dye streaks. The dye streaks roughly return to their initial
positions, with smearing likely as a result of dense dye sinking through
vertically-sheared flow. The effects of slowly sinking dye may have been
observed in Case KR but not Cases I–IV due to the former’s longer duration;
Case KR was 2.6 times longer ($\Delta t=180$ s) than the longest validation
case (Case III, $\Delta t=69$ s). It is also possible a weak meridional
circulation may further smear the dye streaks via transport through the
vertically and radially-sheared azimuthal flow.

The floaters demonstrate a cleaner result: the inset in panel (c) shows a
closeup of the floaters’ initial (pink overlay) and final (grey, actual image)
positions. The difference between initial and final positions is less than
0.5% of the total distance travelled by each floater. Notably, the final
orientation of each triangular floater closely matches initial conditions,
indicating that flow in Case KR is quasi-reversible. We observe only a small
radially-outward shift of floater positions, which we attribute to a weak
meridional circulation.

###  5.2 Magneto-Stokes flow around an obstacle: Comparison to potential flow

Another interesting feature of viscous flows only arises in thin layers, when
viscous forces due to shear across the shallow (vertical) direction are
stronger than those in the extended (horizontal) direction (cf. Goldburg et
al. 1997), such that flow in the extended direction is irrotational as well as
nondivergent. Streamlines of these low-Reynolds-number flows are identical to
those predicted by potential theory appropriate for inviscid, “infinite-
Reynolds-number” flow. Hele-Shaw (1898) first made this connection, using
pressure-driven flow between two glass plates to visualised potential flow
streamlines around an obstacle.

Similar potential flow streamlines may be found in shallow-layer magneto-
Stokes flow around small obstacles placed in the annular channel. The quasi-
static Lorentz force plays an analogous role to pressure in the Hele-Shaw
cell, producing an irrotational, nondivergent free-stream flow. In fact, for
any weakly-conducting fluid immersed in a uniform magnetic field
($\boldsymbol{B}=-B_{0}\boldsymbol{e_{z}}$) and a quasi-static electric field
($\bnabla\times\boldsymbol{E}=0$) without free charges
($\bnabla\bcdot\boldsymbol{E}=0$) or significant gradients along the direction
of the magnetic field ($\partial_{z}\boldsymbol{E}=0$), the resulting Lorentz
force $\boldsymbol{J}\times\boldsymbol{B}$ is nondivergent and irrotational
within the fluid domain:

$\bnabla\times(\boldsymbol{J}\times\boldsymbol{B})=-\sigma
B_{0}\partial_{z}\boldsymbol{E}-\sigma\boldsymbol{B}(\bnabla\bcdot\boldsymbol{E})-\boldsymbol{J}\bcdot\bnabla\boldsymbol{B}+\boldsymbol{J}(\bnabla\bcdot\boldsymbol{B})=0$
(43)

and

$\bnabla\bcdot(\boldsymbol{J}\times\boldsymbol{B})=\sigma\boldsymbol{B}\bcdot(\bnabla\times\boldsymbol{E})-\boldsymbol{J}\bcdot(\bnabla\times\boldsymbol{B})=0,$
(44)

where we have used Ohm’s Law, $\boldsymbol{J}=\sigma\boldsymbol{E}$. The
Lorentz force may then be regarded as the gradient of a harmonic function
(e.g., Aris 1989), $\Phi$:

$\boldsymbol{F_{L}}=\boldsymbol{J}\times\boldsymbol{B}=\bnabla\Phi,\quad\nabla^{2}\Phi=0.$
(45)

For the magneto-Stokes annulus,

$\Phi=\frac{I(t)B_{0}}{2\pi h}\theta,$ (46)

where we have again neglected any fringing of electric field lines due the
finite depth of the conducting fluid. The angular coordinate $\theta$ is
permitted to be multivalued such that it is differentiable everywhere apart
from $r=0$. The presence of any obstacles within the channel may induce a
small correction to (46), which we neglect for analytical convenience by
assuming the obstacle and fluid have equal conductivity.

Flow around such obstacles is then governed by (7) with (45, 46). For
$\Rey,\mathcal{H}\ll 1$, the nondimensional form of (7) becomes

$0=(\mathcal{R}+1)\bnabla_{\text{2D}}\left[\kappa(\tau)\theta-P\right]-\frac{1}{\mathcal{H}}\partial_{\zeta}P\boldsymbol{e_{z}}+\partial_{\zeta}^{2}\boldsymbol{\upsilon},\quad\partial_{\rho}(\rho\upsilon_{\rho})+\partial_{\theta}\upsilon_{\theta}+\rho\partial_{\zeta}\upsilon_{\zeta}=0,$
(47)

in which all quantities are nondimensionalised as in Section 2.1 apart from
pressure, for which an MHD scale is used: $p={I_{0}B_{0}}/({2\pi h})P$. The
vector definitions
$\boldsymbol{\upsilon}=\upsilon_{\rho}\boldsymbol{e_{r}}+\upsilon_{\theta}\boldsymbol{e_{\theta}}+\mathcal{H}\upsilon_{\zeta}\boldsymbol{e_{z}}$
and
$\bnabla_{\text{2D}}=\boldsymbol{e_{r}}\partial_{\rho}+{\rho}^{-1}\boldsymbol{e_{\theta}}\partial_{\theta}$
are employed here for brevity.

For a shallow layer ($\mathcal{H}\ll 1$) we may make a Hele-Shaw approximation
under which velocity is purely horizontal and varies parabolically with depth:

$\boldsymbol{\upsilon}=(2\zeta-\zeta^{2})\boldsymbol{\upsilon}_{\text{2D}}(\rho,\theta),\quad\boldsymbol{\upsilon}_{\text{2D}}\bcdot\boldsymbol{e_{z}}=0$
(48)

Then, using (47), the horizontal velocity may be constructed from potential
functions:

$\boldsymbol{\upsilon}_{\text{2D}}=\frac{\mathcal{R}+1}{2}\bnabla_{\text{2D}}\left[\kappa(\tau)\theta-P(\rho,\theta)\right],$
(49)

where the reduced pressure satisfies Laplace’s equation

$\frac{1}{\rho}\partial_{\rho}(\rho\partial_{\rho}P)+\frac{1}{\rho^{2}}\partial_{\theta}^{2}P=0,$
(50)

and thus may be constructed to enforce the no-penetration condition on the
lateral boundaries via classical techniques of potential theory (e.g., Lamb
1906; Milne-Thomson 1938).

As an example, we consider steady flow ($I=I_{0}$) around a cylindrical
obstacle of diameter $d$ centred at $\theta=0$ and $r=r_{c}$, where
$d<(r_{o}-r_{i})$ and $r_{i}<r_{c}<r_{o}$. The obstacle’s influence on $P$
then depends on the geometric ratios

$\mathcal{C}=r_{c}/r_{o}\quad\text{and}\quad\mathcal{D}=d/r_{o}.$ (51)

An approximate expression for $P$ may be found by using the Milne-Thomson
theorem (Milne-Thomson, 1938) to embed a small circular streamline within the
free-stream. The resulting nondimensional velocity may be expressed as

$\boldsymbol{\upsilon}_{\text{2D}}\approx\frac{\mathcal{R}+1}{2}\bnabla_{\text{2D}}\left\\{\theta+\arctan\left[\frac{\mathcal{D}^{2}\rho\sin\theta}{\mathcal{C}\left(4\mathcal{C}^{2}-\mathcal{D}^{2}+4\rho^{2}\right)+\rho\left(\mathcal{D}^{2}-8\mathcal{C}^{2}\right)\cos\theta}\right]\right\\}.$
(52)

This potential flow solution satisfies no-penetration boundary conditions
exactly on the circular obstacle and approximately on the annulus’ sidewalls.
If the obstacle diameter is $<30\%$ of the channel width, the root mean square
angle between streamlines and sidewalls is $<3^{\circ}$, making (52)
sufficiently accurate to compare to laboratory data.

Figure 8: (a) Dye-visualised laboratory flow around a circular obstacle (white
disk) of nondimensional diameter $\mathcal{D}=0.16$ (Case HS-I). (b) Periodic
vortex shedding around an obstacle with $\mathcal{D}=0.16$ (Case HS-II). (c)
Aperiodic vortex shedding around an obstacle with $\mathcal{D}=0.25$ (Case HS-
II). (d) Closeup photograph of Case HS-I with potential flow streamlines (52)
overlain in black. (e,f) Tracer concentration and nondimensional vorticity
$\omega=[\partial_{\rho}(\rho\upsilon_{\theta})-\partial_{\theta}\upsilon_{\rho}]/\rho$
for a numerical simulation matching Case HS-I. Overlain in black are closely-
spaced streamlines of (52). Gray curves are fictitious streamlines that lie
within the solid obstacle.

Figure 8(a) shows an overhead photograph of a laboratory experiment (Case
HS-I), in which a circular obstacle ($\mathcal{D}=0.16$) is fixed at the
channel mid-gap ($\mathcal{C}=\mathcal{R}/2+1/2$). The obstacle-scale Reynolds
number,

$\Rey_{d}=Ud/\nu,$ (53)

is kept low ($\Rey_{d}=32$) so that flow around the obstacle remains laminar.
A closeup image of Case HS-I in panel (d) shows the streaking of dye (blue)
around the obstacle (white disk) overlain with streamlines of (52) (black
curves). Gray curves within the obstacle (panels a, d–f) show the doublet
component of the potential flow solution that forces streamlines around the
circular boundary.

In the laboratory, the obstacle is electrically insulating and likely produces
a small correction to the Lorentz force (46), which is not accounted for in
(52). However, this effect appears to be minor as the dye streaklines largely
follow the black streamlines in (d).

Case HS-I was matched with a fully nonlinear 2D numerical solution of (6),
(7), and (48) for the Hele-Shaw flow field
$\boldsymbol{\upsilon}_{\text{2D}}(\rho,\theta)$ using Dedalus (code available
at https://github.com/cysdavid/magnetoStokes). The no-slip condition on the
obstacle boundary was approximated using the volume penalty method (Hester et
al., 2021). Panels (e) and (f) plot tracer concentration and vertical
vorticity
$\omega=[\partial_{\rho}(\rho\upsilon_{\theta})-\partial_{\theta}\upsilon_{\rho}]/\rho$
from the matched numerical result. Potential flow streamlines computed from
(52) are overlain in black, and largely coincide with numerical tracer
streaks. The vorticity field reveals non-interacting boundary layers that
produce small deviations from the potential flow solution.

At higher $\Rey_{d}$, these boundary layers detach from the obstacle, as
expected. Panels (b) and (c) of figure 8 show snapshots of laboratory cases
HS-II ($\Rey_{d}=850$) and HS-III ($\Rey_{d}=1800$), respectively. Both flows
exhibit separation, a recirculation zone, and periodic (b) or aperiodic (c)
vortex-shedding of the wake. The transition between the laminar and vortex-
shedding regimes shown in panels (a) and (b) is be controlled by varying the
electric current through the system.

Note that Cases HS-I–III were performed using an array of permanent magnets
enclosed in a wooden container (see figure 1b) instead of RoMag’s solenoid
electromagnet. The tabletop magnet array produces similar results as the
solenoid electromagnet, but permits close-up imaging of the flow.

##  6 Discussion

We analyse magneto-Stokes flow of a thin, free-surface electrolyte layer
within a macroscopic annular channel with conducting sidewalls. An approximate
solution for initially-quiescent flow is validated via DNS and laboratory
experiments within a solenoid bore. We find consistent results using an
inexpensive array of permanent magnets.

In the shallow-layer regime, the bulk flow varies quadratically with depth and
inversely with distance from the annular axis. Near the sidewalls, viscous
forces produce boundary layers that scale with the channel aspect ratio. Flow
spins up slightly slower within the inner boundary layer than predicted by our
model, possibly due to the production of a small meridional circulation as the
free surface deforms. This error from the approximate analytical solution (40)
is likely negligible for all channels but those with the smallest inner
cylinders ($\mathcal{R}\lesssim 0.25$).

The ability to precisely shape magneto-Stokes flow by controlling the channel
geometry and electromagnetic forcing yields a range of possible applications.
For example, our analytical solution makes magneto-Stokes flow a convenient
calibration tool for particle tracking velocimetry (PTV) (e.g., Valenzuela-
Delgado et al. 2018a) and particle image velocimetry (PIV).

The theoretical framework outlined here may also be applied to microfluidic
systems that use the Lorentz force to pump biological and industrial
electrolytes (e.g., Pamme 2006; Pérez-Barrera et al. 2016; Ortiz-Pérez et al.
2017). A major difficulty associated with micro-scale flows is the inefficient
mixing of chemical species at low Reynolds numbers (Lee et al., 2011). West et
al. (2003) address this problem with a miniaturised magneto-Stokes annulus
that enhances mixing between laterally-separated phases by stretching the
phase boundary in the azimuthal direction and thus allowing molecular
diffusion to occur over a larger interfacial area. The leading edge of the dye
spiral in figure 4 is an example of such a phase boundary. The circulatory
flow increases the interfacial area for diffusive flux by adding turns to the
dye spiral.

West et al. (2003) and Gleeson et al. (2004) model micromixing within a narrow
magneto-Stokes annulus ($[1-\mathcal{R}]/\mathcal{H}\to 0$) for which basal
viscous stresses are neglected and the flow is effectively “deep.” However,
the present analysis suggests that the shallow-layer regime may permit even
greater mixing between laterally-separated phases in a micro-scale magneto-
Stokes annulus. One strategy could involve using a wider channel
($\mathcal{R}\ll 1$) with a smaller aspect ratio $\mathcal{H}$ (increasing
$I_{0}B_{0}$ such that $U$ is similar). The sidewall boundary layers would
shrink ($\Delta_{i}\sim\mathcal{H}$), and we expect any phase boundaries to be
rapidly sheared out near the inner cylinder by the $1/r$ azimuthal flow,
increasing the interfacial area over which diffusion can occur.

The magneto-Stokes annulus may also be useful in laboratory simulations of
astrophysical flows. In particular, the $u_{\theta}\sim 1/r$ shallow-layer
magneto-Stokes bulk flow shares the linear centrifugal stability property
($\partial_{r}[u_{\theta}r]\geq 0$) of protoplanetary disks, which have
Keplerian ($u_{\theta}\sim 1/\sqrt{r}$) velocity profiles (Ji & Goodman,
2023). Accretion at the centre of these disks necessitates outwards transport
of angular momentum, which is challenging to explain in the absence of
centrifugal instability (Fromang & Lesur, 2019).

Laboratory studies in Taylor-Couette devices and MHD annular channels have
proven essential to testing proposed mechanisms of angular momentum transport,
which include nonlinear shear-flow instability and magneto-rotational
instability (Ji & Goodman 2023 and ref. therein). Less studied is the vertical
shear instability (VSI), which may occur in accretion disks where radial
thermal gradients induce vertical shear ($\partial_{z}[u_{\theta}]\neq 0$)
through the balance of centrifugation and baroclinic torques (Barker & Latter,
2015). The vertical shear and quasi-Keplerian profile achieved within the
magneto-Stokes annulus may provide an opportunity to study VSI in the
laboratory. A key challenge would be maintaining these features as properties
of the base flow while reducing viscous effects enough to permit instability.

The magneto-Stokes annulus also has promise as a cost-effective device for
fluid dynamics education. Only a handful of MHD educational demonstrations
have been developed in the last five decades (e.g., Heiser & Shercliff 1965;
Melcher & Warren 1965; LeMaster & Shepley 1970; Font & Dudley 2004; Digilov
2007), and most require liquid metal and industrial equipment. In contrast,
the magneto-Stokes system uses saltwater and may be constructed from household
materials for under $100. This simple and robust setup illustrates essential
concepts of both MHD and Stokes flows, including kinematic reversibility and
quasi-2D inviscid behavior under the Hele-Shaw approximation. Abstract
principles of potential flow theory can be illustrated without the need for
expensive Hele-Shaw cells. Additionally, the system’s electromagnetic forcing
makes flow control convenient and precise; educators can demonstrate the
transition from laminar to turbulent regimes of flow past an obstacle by
adjusting a dial on the power supply.

In sum, the magneto-Stokes annulus provides easily generated MHD flows that
permit robust comparisons between theory, DNS, and experiment.

[Supplementary data]A supplementary Python package and accompanying Jupyter
notebook that implement our theoretical predictions are available at
https://github.com/cysdavid/magnetoStokes

[Acknowledgements]The authors thank Ernest Gomis and Andrea Chlarson for their
early laboratory work on the magneto-Stokes annulus. C.S.D. thanks Ellin Zhao
for her advice on image-processing methods, which greatly improved the phase-
boundary tracking code used in this study.

[Funding]This research was supported by the National Science Foundation (EAR
1620649, EAR 1853196).

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The laboratory and DNS data that support our
findings are openly available at https://doi.org/10.5061/dryad.rr4xgxdfc.
Additionally, our dye-tracking velocimetry program and DNS codes may be freely
accessed at https://github.com/cysdavid/magnetoStokes. Also in this GitHub
repository is a Mathematica notebook that reproduces the analytical solution
found in this study.

[Author ORCIDs]C.S. David, https://orcid.org/0009-0006-9774-1471; E.W. Hester,
https://orcid.org/0000-0003-1651-9141; Y. Xu,
https://orcid.org/0000-0001-9123-124X; J.M. Aurnou,
https://orcid.org/0000-0002-8642-2962

[Author contributions]C.S.D. developed the model and found the solution
presented here. E.W.H. wrote the DNS codes for both the axisymmetric standard
cases and the non-axisymmetric Hele-Shaw obstacle flow; C.S.D ran these codes.
C.S.D. and Y.X. carried out laboratory experiments, with Y.X. providing key
insights into flow visualisation. J.M.A. conceived the idea for this study and
guided its development. C.S.D. prepared the first draft of this manuscript,
and all authors contributed to the final version.

## Appendix A

In Section 2.2.2, we provide an approximate solution to (21) with boundary
conditions $u^{\prime}_{\theta}=0$ at $\rho=\mathcal{R}$, $\rho=1$, $\zeta=0$
and $\partial_{\zeta}u^{\prime}_{\theta}=0$ at $\zeta=1$, and initial
condition $u^{\prime}_{\theta}=-\bar{\upsilon}_{\theta}$ at $\tau=0$. Here, we
supply the full solution using a general Fourier-Bessel series of first order:

$\upsilon_{\theta}^{\prime}(\rho,\zeta,\tau)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}C_{mn}R_{m}(\rho)\sin\left(k_{n}\zeta\right)\exp\left({-\frac{\mathcal{H}^{2}\mu_{m}^{2}+k_{n}^{2}}{\Rey}\tau}\right),$
(54)

where each radial eigenfunction is a combination of first-order Bessel
functions of first and second kind:

$R_{m}(\rho)=\frac{\rmJ_{1}\left(\mu_{m}\rho\right)}{\rmJ_{1}(\mu_{m})}-\frac{\rmY_{1}\left(\mu_{m}\rho\right)}{\rmY_{1}(\mu_{m})},$
(55)

and each eigenvalue $\mu_{m}$ is the $m$-th smallest positive root of

$\rmY_{1}(\mu)\rmJ_{1}(\mu\mathcal{R})-\rmJ_{1}(\mu)\rmY_{1}(\mu\mathcal{R})=0.$
(56)

Table 5 reports the first ten eigenvalues $\mu_{m}$ for
$\mathcal{R}=0.25,\;0.5,\;0.75$ computed to five decimal places.

The coefficients $C_{mn}$ may be found by applying the initial condition:

$\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}C_{mn}R_{m}(\rho)\sin\left(k_{n}\zeta\right)=\upsilon_{\theta}^{\prime}(\rho,\zeta,0)=-\bar{\upsilon}_{\theta}(\rho,\zeta)=-\sum_{n=1}^{\infty}u_{n}(\rho)\sin\left(k_{n}\zeta\right),$
(57)

where

$u_{n}(\rho)=\frac{2(\mathcal{R}+1)}{k_{n}^{2}\mathcal{H}}\left[\frac{\mathcal{H}}{k_{n}\rho}-A_{n}\rmI_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)-B_{n}\rmK_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)\right]$
(58)

as in Section 2. Then,

$C_{mn}=-\frac{\langle u_{n},R_{m}\rangle}{\langle
R_{m},R_{m}\rangle}\quad\text{where}\quad\langle
f,g\rangle=\int_{\mathcal{R}}^{1}\rho f(\rho)g(\rho)d\rho.$ (59)

Using the properties of Bessel functions, the inner products in (59) may be
computed as follows:

$\begin{split}\langle
R_{m},R_{m}\rangle=\rho^{2}\biggr{[}&\frac{\rmJ_{1}\left(\mu\rho\right){}^{2}-\rmJ_{0}\left(\mu\rho\right)\rmJ_{2}\left(\mu\rho\right)}{\rmJ_{1}(\mu){}^{2}}+\\\
&\frac{\rmJ_{2}\left(\mu\rho\right)\rmY_{0}\left(\mu\rho\right)-2\rmJ_{1}\left(\mu\rho\right)\rmY_{1}\left(\mu\rho\right)+\rmJ_{0}\left(\mu\rho\right)\rmY_{2}\left(\mu\rho\right)}{2\rmJ_{1}(\mu)\rmY_{1}(\mu)}\biggr{]}\Biggr{\lvert}_{\rho=\mathcal{R}}^{\rho=1}\end{split}$
(60)

and

$\begin{split}{\langle
u_{n},R_{m}\rangle}=\frac{2(\mathcal{R}+1)}{k_{n}^{2}\mathcal{H}}\biggr{[}&\frac{\mathcal{H}}{k_{n}}\langle
1/\rho,R_{m}(\rho)\rangle\\\
&-A_{n}\langle\rmI_{1}(k_{n}\rho/\mathcal{H}),R_{m}(\rho)\rangle-
B_{n}\langle\rmK_{1}(k_{n}\rho/\mathcal{H}),R_{m}(\rho)\rangle\biggr{]}.\end{split}$
(61)

Within (61), the first inner product evaluates to

$\langle
1/\rho,R_{m}(\rho)\rangle=\frac{1}{\mu_{m}}\left(\frac{1-\rmJ_{0}\left(\rho\mu_{m}\right)}{\rmJ_{1}\left(\mu_{m}\right)}+\frac{\rmY_{0}\left(\rho\mu_{m}\right)}{\rmY_{1}\left(\mu_{m}\right)}\right)\Biggr{\lvert}_{\rho=\mathcal{R}}^{\rho=1},$
(62)

the second evaluates to

$\begin{split}\langle\rmI_{1}(k_{n}\rho/\mathcal{H}),R_{m}(\rho)\rangle&=\frac{\mathcal{H}\rho}{\rmJ_{1}(\mu_{m})\rmY_{1}(\mu_{m})\left(\mathcal{H}^{2}\mu_{m}^{2}+k_{n}^{2}\right)}\\\
&\times\biggr{\\{}k_{n}\rmI_{2}\left(\frac{k_{n}\rho}{\mathcal{H}}\right)\left[\rmY_{1}(\mu_{m})\rmJ_{1}(\mu_{m}\rho)-\rmJ_{1}(\mu_{m})\rmY_{1}(\mu_{m}\rho)\right]\\\
&+\mathcal{H}\mu_{m}\rmI_{1}\left(\frac{k_{n}\rho\
}{\mathcal{H}}\right)\left[\rmY_{1}(\mu_{m})\rmJ_{2}(\mu_{m}\rho)-\rmJ_{1}(\mu_{m})\rmY_{2}(\mu_{m}\rho)\right]\biggr{\\}}\Biggr{\lvert}_{\rho=\mathcal{R}}^{\rho=1},\end{split}$
(63)

and the third evaluates to

$\begin{split}\langle\rmK_{1}(k_{n}\rho/\mathcal{H}),R_{m}(\rho)\rangle&=\frac{\mathcal{H}\rho}{\rmJ_{1}(\mu_{m})\rmY_{1}(\mu_{m})\left(\mathcal{H}^{2}\mu_{m}^{2}+k_{n}^{2}\right)}\\\
&\times\biggr{\\{}k_{n}\rmK_{2}\left(\frac{k_{n}\rho}{\mathcal{H}}\right)\left[\rmJ_{1}(\mu_{m})\rmY_{1}(\mu_{m}\rho)-\rmY_{1}(\mu_{m})\rmJ_{1}(\mu_{m}\rho)\right]\\\
&+\mathcal{H}\mu_{m}\rmK_{1}\left(\frac{k_{n}\rho}{\mathcal{H}}\right)\left[\rmY_{1}(\mu_{m})\rmJ_{2}(\mu_{m}\rho)-\rmJ_{1}(\mu_{m})\rmY_{2}(\mu_{m}\rho)\right]\biggr{\\}}\Biggr{\lvert}_{\rho=\mathcal{R}}^{\rho=1}.\end{split}$
(64)

$m$ | $\mathcal{R}=0.25$ | $\mathcal{R}=0.5$ | $\mathcal{R}=0.75$
---|---|---|---
1 |  04.44751 |  06.39316 |  012.60560
2 |  08.53692 | 12.62470 |  025.15256
3 | 12.67896 | 18.88893 |  037.71235
4 | 16.84163 | 25.16241 |  050.27542
5 | 21.01395 | 31.43971 |  062.83981
6 | 25.19147 | 37.71895 |  075.40485
7 | 29.37208 | 43.99932 |  087.97028
8 | 33.55467 | 50.28038 | 100.53594
9 | 37.73861 | 56.56191 | 113.10176
10 | 41.92350 | 62.84378 | 125.66768
Table 5: Radial eigenvalues $\mu_{m}$ associated with the full time-dependent
solution (65). Each $\mu_{m}$ is computed as the $m$-th smallest positive root
of (56).

The full solution to the IVP posed in Section 2.2 may be formed by summing
(25) and (54):

$\displaystyle\upsilon_{\theta}(\rho,\zeta,\tau)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\Biggr{\\{}\frac{2(\mathcal{R}+1)}{k_{n}^{2}\mathcal{H}}\left[\frac{\mathcal{H}}{k_{n}\rho}-A_{n}\rmI_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)-B_{n}\rmK_{1}\left(\frac{k_{n}}{\mathcal{H}}\rho\right)\right]$
(65)
$\displaystyle\quad\quad\mbox{}+\sum_{m=1}^{\infty}C_{mn}R_{m}(\rho)\exp\left({-\frac{\mathcal{H}^{2}\mu_{m}^{2}+k_{n}^{2}}{\Rey}\tau}\right)\Biggr{\\}}\sin\left(k_{n}\zeta\right).$

Setting $\mathcal{H}^{2}\mu_{m}^{2}+k_{n}^{2}\approx k_{1}$ in the exponential
term above reduces (65) to the approximate solution (40). The solution
procedure outlined above and in Section 2.2.1 is reproduced in a Mathematica
notebook available at https://github.com/cysdavid/magnetoStokes.

## References

  * Anderson et al. (1959) Anderson, Oscar, Baker, William R., Bratenahl, Alexander, Furth, Harold P. & Kunkel, Wulf B. 1959 Hydromagnetic Capacitor. Journal of Applied Physics 30 (2), 188–196.
  * Aris (1989) Aris, Rutherford 1989 Vectors, tensors, and the basic equations of fluid mechanics. Dover publications.
  * Barker & Latter (2015) Barker, A. J. & Latter, H. N. 2015 On the vertical-shear instability in astrophysical discs. Monthly Notices of the Royal Astronomical Society 450 (1), 21–37.
  * Baylis & Hunt (1971) Baylis, J. A. & Hunt, J. C. R. 1971 MHD flow in an annular channel; theory and experiment. Journal of Fluid Mechanics 48 (3), 423–428.
  * Bradski (2000) Bradski, G. 2000 The OpenCV Library. Dr. Dobb’s Journal of Software Tools 25 (11), 120–125.
  * Braginsky (1959) Braginsky, Stanislav I. 1959 Magnetohydrodynamics of Weakly Conducting Fluids. Soviet Journal of Experimental and Theoretical Physics 37 (5), 1417–1430.
  * Burns et al. (2020) Burns, Keaton J., Vasil, Geoffrey M., Oishi, Jeffrey S., Lecoanet, Daniel & Brown, Benjamin P. 2020 Dedalus: A flexible framework for numerical simulations with spectral methods. Physical Review Research 2 (2), 023068.
  * Canny (1986) Canny, John 1986 A Computational Approach to Edge Detection. IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-8 (6), 679–698.
  * Davidson (2016) Davidson, P. A. 2016 Introduction to Magnetohydrodynamics, 2nd edn. Cambridge University Press.
  * Digilov (2007) Digilov, Rafael M. 2007 Making a fluid rotate: Circular flow of a weakly conducting fluid induced by a Lorentz body force. American Journal of Physics 75 (4), 361–367.
  * Early & Dow (1950) Early, H. C. & Dow, W. G. 1950 Supersonic Wind at Low Pressures Produced by Arc in Magnetic Field. Physical Review 79 (1), 186–186.
  * Favier et al. (2011) Favier, B., Godeferd, F. S., Cambon, C., Delache, A. & Bos, W. J. T. 2011 Quasi-static magnetohydrodynamic turbulence at high Reynolds number. Journal of Fluid Mechanics 681, 434–461.
  * Font & Dudley (2004) Font, Gabriel I. & Dudley, Scott C. 2004 Magnetohydrodynamic Propulsion for the Classroom. The Physics Teacher 42 (7), 410–415.
  * Fromang & Lesur (2019) Fromang, S. & Lesur, G. 2019 Angular momentum transport in accretion disks: a hydrodynamical perspective. EAS Publications Series 82, 391–413.
  * Gleeson et al. (2004) Gleeson, James P., Roche, Olivia M., West, Jonathan & Gelb, Anne 2004 Modelling Annular Micromixers. SIAM Journal on Applied Mathematics 64 (4), 1294–1310.
  * Goldburg et al. (1997) Goldburg, W.I., Rutgers, M.A. & Wu, X.-l. 1997 Experiments on turbulence in soap films. Physica A: Statistical Mechanics and its Applications 239 (1-3), 340–349.
  * Greenspan & Howard (1963) Greenspan, H. P. & Howard, L. N. 1963 On a time-dependent motion of a rotating fluid. Journal of Fluid Mechanics 17 (03), 385.
  * Heiser & Shercliff (1965) Heiser, W. H. & Shercliff, J. A. 1965 A simple demonstration of the Hartmann layer. Journal of Fluid Mechanics 22 (04), 701.
  * Hele-Shaw (1898) Hele-Shaw, H. S. 1898 The Flow of Water. Nature 58 (1489), 34–36.
  * Hester et al. (2021) Hester, Eric W., Vasil, Geoffrey M. & Burns, Keaton J. 2021 Improving accuracy of volume penalised fluid-solid interactions. Journal of Computational Physics 430, 110043.
  * Hunt & Stewartson (1965) Hunt, J. C. R. & Stewartson, K. 1965 Magnetohydrodynamic flow in rectangular ducts. II. Journal of Fluid Mechanics 23 (03), 563.
  * Isdale & Morris (1972) Isdale, J.D. & Morris, R. 1972 Physical properties of sea water solutions: density. Desalination 10 (4), 329–339.
  * Isdale et al. (1972) Isdale, J.D., Spence, C.M. & Tudhope, J.S. 1972 Physical properties of sea water solutions: viscosity. Desalination 10 (4), 319–328.
  * Ji & Goodman (2023) Ji, H. & Goodman, J. 2023 Taylor–Couette flow for astrophysical purposes. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 381 (2246), 20220119.
  * Khal’zov & Smolyakov (2006) Khal’zov, I. V. & Smolyakov, A. I. 2006 On the calculation of steady-state magnetohydrodynamic flows of liquid metals in circular ducts of a rectangular cross section. Technical Physics 51 (1), 26–33.
  * Knaepen & Moreau (2008) Knaepen, Bernard & Moreau, René 2008 Magnetohydrodynamic Turbulence at Low Magnetic Reynolds Number. Annual Review of Fluid Mechanics 40 (1), 25–45.
  * Lamb (1906) Lamb, Horace 1906 Hydrodynamics, 3rd edn. Cambridge University Press.
  * Lee et al. (2011) Lee, Chia-Yen, Chang, Chin-Lung, Wang, Yao-Nan & Fu, Lung-Ming 2011 Microfluidic Mixing: A Review. International Journal of Molecular Sciences 12 (5), 3263–3287.
  * LeMaster & Shepley (1970) LeMaster, E. W. & Shepley, L. C. 1970 A Simple MHD Pump. American Journal of Physics 38 (3), 389–390.
  * Melcher & Warren (1965) Melcher, J. R. & Warren, E. P. 1965 Demonstration of Magnetic Flux Constraints and a Lumped Parameter Alfven Wave. IEEE Transactions on Education 8 (2), 41–47.
  * Milne-Thomson (1938) Milne-Thomson, L. M. 1938 Theoretical Hydrodynamics, 1st edn. Macmillan & Co. Ltd.
  * Moresco & Alboussire (2004) Moresco, P. & Alboussire, T. 2004 Experimental study of the instability of the Hartmann layer. Journal of Fluid Mechanics 504, 167–181.
  * Ortiz-Pérez et al. (2017) Ortiz-Pérez, A. S., García-Ángel, V., Acuña-Ramírez, A., Vargas-Osuna, L. E., Pérez-Barrera, J. & Cuevas, S. 2017 Magnetohydrodynamic flow with slippage in an annular duct for microfluidic applications. Microfluidics and Nanofluidics 21 (8), 138.
  * Pamme (2006) Pamme, Nicole 2006 Magnetism and microfluidics. Lab on a Chip 6 (1), 24–38.
  * Park (1964) Park, Kilho 1964 Partial equivalent conductance of electrolytes in sea water. Deep Sea Research and Oceanographic Abstracts 11 (5), 729–736.
  * Poyé et al. (2020) Poyé, A., Agullo, O., Plihon, N., Bos, W. J. T., Desangles, V. & Bousselin, G. 2020 Scaling laws in axisymmetric magnetohydrodynamic duct flows. Physical Review Fluids 5 (4), 043701.
  * Pérez-Barrera et al. (2016) Pérez-Barrera, J., Ortiz, A. & Cuevas, S. 2016 Analysis of an Annular MHD Stirrer for Microfluidic Applications. In Recent Advances in Fluid Dynamics with Environmental Applications (ed. Jaime Klapp, Leonardo Di G. Sigalotti, Abraham Medina, Abel López & Gerardo Ruiz-Chavarría), pp. 275–288. Springer International Publishing, series Title: Environmental Science and Engineering.
  * Pérez-Barrera et al. (2015) Pérez-Barrera, J., Pérez-Espinoza, J. E., Ortiz, A., Ramos, E. & Cuevas, S. 2015 Instability of electrolyte flow driven by an azimuthal Lorentz force. Magnetohydrodynamics 51 (2), 203–214.
  * Stelzer et al. (2015a) Stelzer, Zacharias, Cébron, David, Miralles, Sophie, Vantieghem, Stijn, Noir, Jérôme, Scarfe, Peter & Jackson, Andrew 2015a Experimental and numerical study of electrically driven magnetohydrodynamic flow in a modified cylindrical annulus. I. Base flow. Physics of Fluids 27 (7), 077101.
  * Stelzer et al. (2015b) Stelzer, Zacharias, Miralles, Sophie, Cébron, David, Noir, Jérôme, Vantieghem, Stijn & Jackson, Andrew 2015b Experimental and numerical study of electrically driven magnetohydrodynamic flow in a modified cylindrical annulus. II. Instabilities. Physics of Fluids 27 (8), 084108.
  * Taylor (1923) Taylor, Geoffrey I. 1923 VIII. Stability of a viscous liquid contained between two rotating cylinders. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 223 (605-615), 289–343.
  * Taylor (1967) Taylor, Geoffrey I. 1967 Film Notes for Low Reynolds-Number Flows. National Committee for Fluid Mechanics Film.
  * Valenzuela-Delgado et al. (2018a) Valenzuela-Delgado, Mónica, Flores-Fuentes, Wendy, Rivas-López, Moisés, Sergiyenko, Oleg, Lindner, Lars, Hernández-Balbuena, Daniel & Rodríguez-Quiñonez, Julio 2018a Electrolyte Magnetohydrondyamics Flow Sensing in an Open Annular Channel—A Vision System for Validation of the Mathematical Model. Sensors 18 (6), 1683.
  * Valenzuela-Delgado et al. (2018b) Valenzuela-Delgado, M., Ortiz-Pérez, A.S., Flores-Fuentes, W., Bravo-Zanoguera, M.E., Acuña-Ramírez, A., Ocampo-Díaz, J.D., Hernández-Balbuena, D., Rivas-López, M. & Sergiyenko, O. 2018b Theoretical and experimental study of low conducting fluid MHD flow in an open annular channel. International Journal of Heat and Mass Transfer 127, 322–331.
  * Verma (2017) Verma, Mahendra K 2017 Anisotropy in Quasi-Static Magnetohydrodynamic Turbulence. Reports on Progress in Physics 80 (8), 087001.
  * Wang & Ruuth (2008) Wang, Dong & Ruuth, Steven J. 2008 Variable step-size implicit-explicit linear multistep methods for time-dependent partial differential equations. Journal of Computational Mathematics 26 (6), 838–855, publisher: Institute of Computational Mathematics and Scientific/Engineering Computing.
  * West et al. (2003) West, Jonathan, Gleeson, James P., Alderman, John, Collins, John K. & Berney, Helen 2003 Structuring laminar flows using annular magnetohydrodynamic actuation. Sensors and Actuators B: Chemical 96 (1-2), 190–199.
  * Xu et al. (2022) Xu, Yufan, Horn, Susanne & Aurnou, Jonathan M. 2022 Thermoelectric precession in turbulent magnetoconvection. Journal of Fluid Mechanics 930, A8.
  * Yuen et al. (1990) Yuen, Hk, Princen, J, Illingworth, J & Kittler, J 1990 Comparative study of Hough Transform methods for circle finding. Image and Vision Computing 8 (1), 71–77.