Exact and Approximate Solutions for Magnetohydrodynamic Flow Control in Hele-Shaw Cells

Kyle I. McKee\aff1,2 \corresp [email protected] \aff1Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA \aff2Gulliver Laboratory, ESPCI Paris, PSL University, 75231 Paris cedex 05, France

Abstract

Consider the motion of a thin layer of electrically conducting fluid, between two closely spaced parallel plates, in a classical Hele-Shaw geometry. Furthermore, let the system be immersed in a uniform external magnetic field (normal to the plates) and let electrical current be driven between conducting probes immersed in the fluid layer. In the present paper, we analyse the ensuing fluid flow at low Hartmann numbers. We first elucidate the mechanism of flow generation both physically and mathematically. We proceed by presenting mathematical solutions for a class of canonical multiply-connected geometries, in terms of the prime function developed by crowdy2020solving. Notably, those solutions can be written explicitly as series, and are thus exact, in doubly-connected geometries. Note that in higher connectivities, the prime function must be evaluated numerically. We then demonstrate how recently developed fast numerical methods may be applied to accurately determine the flow-field in arbitrary geometries when exact solutions are inaccessible.

1 Introduction

Magnetohydrodynamic flows find applications in a variety of fields from geophysics and astrophysics to metallurgy, in situations where the fluid under study interacts significantly with external or self-induced electromagnetic fields (moffatt1978field; davidson2002introduction). Although magnetic effects tend to transform the Navier-Stokes equations into an even more formidable form, they also generate a host of unique phenomena that have been studied over the past century including Alfvén waves (alfven1942existence) and geodynamos (moffatt1978field; moffatt2019self). In the present paper, we analyse mathematically the manner in which electrical and magnetic effects modify arguably the simplest possible base flow — a two-dimensional potential flow as is physically realized in a Hele-Shaw cell.

First, consider a Hele-Shaw flow in the absence of electromagnetic effects. The flow is obtained through the solution of a two-dimensional boundary-value problem for a harmonic pressure-field whose negative gradient gives the velocity field. Since the pressure field must be single-valued, the fluid flow must possess zero circulation, thus removing the usual paradox of undetermined circulations present in aerofoil theory (see gonzalez2022variational for a recent development and survey of the aerofoil problem). Since pressure-driven Hele-Shaw flows possess exactly zero circulation around any closed contour, circulation cannot be used to mix or control such flows. Since the Hele-Shaw cell constitutes a model for flows through porous media and microfluidic devices, it is of technological interest to overcome the no-circulation limitation. We now outline some relevant literature regarding circulation generation via electromagnetic effects.

moffatt1991electromagnetic discussed mechanisms for stirring using time-dependent electromagnetic fields, in general three-dimensional contexts. The mechanisms presented therein rely mainly on truly three-dimensional aspects of flow, and the Hele-Shaw limit was not considered. Recently, mirzadeh2020vortices showed that circulation can be generated in electro-osmotic Hele-Shaw flows if the gap thickness is made inhomogeneous. Henceforth, when discussing Hele-Shaw cells, we restrict our attention to the case of uniform thickness.

bau2001minute and zhong2002magneto fabricated Hele-Shaw type devices containing thin layers of electrolytes, through which electrical current was driven between electrodes immersed in the fluid. The former work placed electrodes on the base of the device whereas the latter placed electrodes on the walls of a concentric annlulus. When the devices were placed in a uniform magnetic field, Lorentz forces induced a fluid flow in the bulk. yi2002magnetohydrodynamic showed that under certain conditions, the application of periodic AC electrical currents in similar devices may lead to chaotic advection. Later, homsy2005high manufactured a high current density DC microfluidic pump which found applications in nuclear magnetic resonance (homsy2007magnetohydrodynamic); notably, the pump avoided significant Joule heating. A recent comprehensive review of applications of magnetohydrodynamics in the context of microfluidics is given by bau2022applications.

The papers discussed in the paragraph above were mainly experimental in nature. In the cases where fluid flows were analysed analytically, such as in bau2001minute and zhong2002magneto, solutions were obtained as infinite series in highly specialized coordinate systems, using methods only applicable to special geometries. One aim of the present paper is to demonstrate how conformal maps obviate these complicated analyses and yield a class of flow geometries with closed-form solutions; in doubly-connected problems, the prime function has a series representation and solutions are exact (baddoo2020exact).

More recently, david2023magnetostokes established a mathematical analogy, in magneto-hydrodynamics, to reversible Stokes flow in a Taylor-Couette geometry. Experiments presented therein beautifully reproduced the classic kinematic reversibility experiments of kinrev without the need for moving boundaries. At the end of their paper, the authors noted that in the limit of a shallow fluid layer, the flow resembles potential flow. Note that their experiments were conducted with a free surface and no lid. Although the authors correctly pointed out the potential flow limit, no analytical solutions or further theoretical developments were explored in their work.

In the present paper, we analyse generally the flow of a conducting fluid in a Hele-Shaw cell, immersed in an external uniform and constant magnetic field that is directed normal to the cell walls, $\boldsymbol{B}=B_{0}\hat{\boldsymbol{z}}$ (see figure 1). We consider scenarios where flow is driven by applying voltages to conducting probes that are immersed in the fluid, the probes being impermiable to fluid flow. We assume an ohmic fluid to model the electrical current flow, $\boldsymbol{J}=\sigma\boldsymbol{E}$ , where $\sigma$ is the electrical conductivity of the fluid and $\boldsymbol{E}$ is the electric field experienced by the fluid. The presence of electrical current leads to Lorentz forces acting on the fluid bulk, $\boldsymbol{F}=\boldsymbol{J}\times\boldsymbol{B}$ , which ultimately induces a fluid flow. We also consider the effect of adding obstacles to the flow that are either electrical insulators or conductors. We present mathematical solutions for a large class of multiply-connected geometries in terms of the prime function given by crowdy2020solving. In general geometries where exact solutions are not possible, we show how highly accurate approximate solutions can be obtained using series solution methods as described by trefethen2018series.

The remainder of this paper is arranged as follows. In §2, we discuss the simplifying assumptions in our mathematical formulation, and subsequently outline the physical mechanism driving the fluid flow. In §3, we present a complex variables formulation of the model described in §2. We proceed by deriving mathematical solutions for the fluid flow in a variety of multiply-connected geometries, using the framework of crowdy2020solving. One such solution gives the flow in a geometry explored experimentally by david2023magnetostokes. In §3.5, another solution is compared to a new experiment which serves as further motivation for the present theoretical study. In §4, we show how more general geometries can be solved to high accuracy using series solutions.

Refer to caption — Figure 1: (a) Schematic of the system under consideration in the present paper. Voltages are applied to conducting bodies immersed in a thin layer of conducting fluid, of thickness $h$ , that is bounded above and below by parallel walls. The entire system is immersed in a uniform magnetic field oriented along the $z$ -axis and normal to the bounding walls. (b) Top-view of the flow geometry. Each conducting probe serves as an impermeable obstacle in the fluid flow. Each probe is held at a fixed voltage. We also investigate the effect of insulating obstacles in the flow, which obey the zero Neumann condition, $\nabla V\cdot\boldsymbol{n}=0$ , in place of the Dirichlet condition satisfied by conductors, where $\boldsymbol{n}$ is the unit normal vector to the surface of the obstacle.

2 Assumptions and Physical Picture

Consider a thin layer of conducting fluid occupying the region between two rigid walls separated by a distance $h$ , as in figure 1(a). The conducting fluid may be taken, for example, to be saltwater or liquid mercury. Perfectly conducting probes, held at fixed voltages, penetrate the entire cell thickness $h$ . Meanwhile, the system is immersed in an external magnetic field perpendicular to the flow, $\boldsymbol{B}=B_{0}\hat{\boldsymbol{z}}$ . The application of a voltage difference between two conductors gives rise to an electric field in the conducting fluid, which induces an electrical current flow according to Ohm’s law. The external magnetic field then induces a Lorentz force on fluid parcels which gives rise to a fluid flow whose direction can be predicted using the right-hand rule.

Restricting our attention to steady and quasi-steady problems, the voltage $V(\boldsymbol{x})$ is determined generally as the solution to a Poisson equation. In an uncharged ohmic material, the Poisson equation reduces to the Laplace equation so that $V(\boldsymbol{x})$ must be a harmonic function. We now justify the proposed physical picture.

In the studies mentioned in §1, as well as our experiment in §3.5, the conducting fluid is a binary electolyte and so the “uncharged” assumption requires extra justification. A binary electrolyte comprises two oppositely charged species, whose concentrations we denote $c^{+}$ and $c^{-}$ . If at any point of space $c^{+}\neq c^{-}$ holds, then the Poisson equation does not reduce to the Laplace equation, since the charge density (source term) is not identically zero in all space. In such a case, $V(\boldsymbol{x})$ fails to be harmonic across the entire domain. To ensure that $V(\boldsymbol{x})$ is indeed harmonic, we must invoke the so-called local electroneutrality condition, $c^{+}=c^{-}$ , which is a standard assumption in the modelling of electrolytes (zaltzman2007electro). This assumption is valid in the bulk of the fluid but is violated in the electric double layer, within a Debye length of conducting boundaries. In the present study, we make no effort to model double layer effects. Instead we treat boundaries as simplified ideal conductors and insulators, and the bulk fluid as a locally electroneutral ohmic conductor. We now proceed under the assumption that the electrical potential can be modelled as a harmonic function in the fluid domain.

Under the stated assumptions, the electric field as measured in the lab frame is given by the gradient of a harmonic potential, $\boldsymbol{E}_{\mathrm{lab}}=-\nabla V$ . The electric field drives an electrical current according to Ohm’s law, $\boldsymbol{J}=\sigma\boldsymbol{E}$ , where $\boldsymbol{E}$ is the electric field felt in the frame of a fluid parcel. The electric field that a fluid parcel experiences is related to the electric field in the lab frame by the relation $\boldsymbol{E}=\boldsymbol{E}_{\mathrm{lab}}+\boldsymbol{u}\times\boldsymbol{B}$ , owing to the fact that the electric field is not invariant under Gallilean transformations (see moffatt1978field). However, we will now show that the second term is negligible in the context of the Hele-Shaw flow under consideration here. In our system, typical voltage differences between probes are $1V$ ; voltages higher than roughly $1.23V$ create bubbles due to water electrolysis. With typical separation distances between electrical probes being on the order of centimeters, the typical electric field magnitude is $E_{0}\approx 1V/\mathrm{cm}$ . Typical velocities in the cell are $U_{0}\approx 1\mathrm{mm}/\mathrm{s}$ , while the maximum magnetic field is roughly $B_{0}\approx 2340\mathrm{G}$ . Hence, the relative importance of $\boldsymbol{u}\times\boldsymbol{B}$ as compared to $\boldsymbol{E}_{\mathrm{lab}}$ scales as $U_{0}B_{0}/E_{0}=\mathcal{O}(10^{-6})$ and we thus safely neglect the cross product term. Henceforth, we use $\boldsymbol{E}$ to denote the electric field and ignore any distinction between the field experienced in different frames.

The equations of conservation of fluid momentum and mass then become, {subeqnarray} ρ(∂u∂t+u⋅∇u) & = μ∇^2 u-∇P+σE×B,
∇⋅ u = 0.

Since the cell thickness, $h$ , is much smaller than the lateral extent of the system, the Hele-Shaw approximation is justified, implying that viscous forces dominate inertial forces (see batchelor1967introduction). Hence, the left side of (2a) can be neglected, and the flow is determined by a balance of viscous forces and both pressure and magnetic driving forces. In order to invoke the usual Hele-Shaw approximation, it is important to ensure that the flow across the entire gap thickness, $h$ , is fully developed. The presence of a magnetic field has the ability to shrink the boundary layer so that the fully-developed parabolic velocity profile may not be reached (davidson2002introduction, figure 5.10 on pp. 153); see also the work of rossow1958flow for an interesting application of this concept. In our system, the Hartmann number is small, $Ha=\sqrt{B_{0}^{2}h^{2}\sigma/\mu}\approx 0.01$ , indicating that magnetic dissipation is neglibible and the flow does indeed adopt the fully-developed viscous parabolic profile; our calculations use the typical gap thickness, $h\approx 0.7\mathrm{mm}$ (see §3.5). For comparison, liquid mercury in a Hele-Shaw cell of the same thickness, and subject to the same magnetic field, has a larger Hartman number $Ha\approx 4$ ; thus, a thinner gap is necessary to attain the Hele-Shaw limit ( $Ha\ll 1$ ) for the liquid metal.

Under the Hele-Shaw approximation, the top-down velocity becomes two-dimensional and we can immediately write {subeqnarray} u & = 12μ( -∇P+σB_0E×^z)z(h-z),
∇⋅ u = 0.

To make further progress, we note that the Lorentz force, $\boldsymbol{F}=-\sigma B_{0}\nabla V\times\hat{\boldsymbol{z}}$ , is clearly solenoidal and irrotational so that $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{F}=0$ and $\boldsymbol{\nabla}\times\boldsymbol{F}=0$ . It is therefore possible to represent it as the gradient of a harmonic function, $\boldsymbol{F}=\nabla\phi$ , where $\nabla^{2}\phi=0$ . In general $\phi$ need not be single-valued. Equation (2a) then reduces to

\boldsymbol{u}=\frac{z(h-z)}{2\mu}\nabla\left(\phi(x,y)-P(x,y)\right),

(1)

where the gradient is two-dimensional. Henceforth, we shall suppress the vertical structure of the flow, $z(h-z)$ , for convenience since we are only concerned with a top-down view of the flow. Note that since $\phi$ is a harmonic function, (2b) and (1) together imply that $P$ must be harmonic too.

While $P$ must be a single-valued function, there is no such restriction on $\phi$ . In fact, the main manner in which the Lorentz force generates flow, as will be shown, is by inducing flow circulation which actually requires $\phi$ to be multi-valued. Physically, the circulation is induced as follows. Since electric field lines must exit perpendicular to conducting surfaces, so too must the current density, $\boldsymbol{J}$ , according to Ohm’s law. The right-hand rule then reveals that the Lorentz force induces a circulatory flow around each conductor.

3 Complex Variables Formulation and Solution Procedure

We proceed by formulating the boundary value problem for the two-dimensional velocity field, $\tilde{\boldsymbol{u}}\equiv\nabla(\phi-P)$ , using a complex variables approach. The solution procedure for obtaining $\tilde{\boldsymbol{u}}$ is as follows. First, one must solve the electrostatic problem for $V(x,y)$ in the fluid domain. Once $V$ is known, the potential $\phi$ , which describes the Lorentz force according to $\boldsymbol{F}=\nabla\phi$ , may be obtained. Finally, a single-valued pressure field $P$ , as appears in (1), must be obtained to enforce the impermeability condition on the surface of each flow obstacle. The procedure is outlined in detail in the remainder of §3.

3.1 Complex Variables Formulation

Consider the electrostatic boundary value problem in the domain illustrated in figure 1(b). We seek a real harmonic function $V(x,y)$ in the fluid domain satisfying constant Dirichlet conditions on each conducting boundary. More generally, we consider also the addition of perfectly insulating obstacles, on the surface of which $V$ satisfies a zero Neumann boundary condition.

On each conducting surface, $\partial B_{i}$ , constant Dirichlet boundary data is prescribed so that $V=V_{j}$ for $j\in\{1,2,...,N_{C}\}$ , where $N_{C}$ is the total number of conductors. On insulating boundaries, $\partial D_{k}$ , the normal derivative is required to vanish so that $\nabla V\cdot\boldsymbol{n}=0$ where $k\in\{1,2,...,N_{I}\}$ and $N_{I}$ is the total number of electrical insulators. To make progress, we now express the boundary value problem in the language of complex variables.

We may take $V$ to be the real part of an analytic function $W_{E}(z)$ , where $z=x+iy$ . We require $W_{E}(z)$ to be analytic over the entire fluid domain to ensure the harmonicity of $V(x,y)$ . We thus seek an analytic function $W_{E}$ with the properties {subeqnarray} Re{W_E(z)} & = V_j , z ∈∂B_j , j∈{1,2,…, N_C},
Im{ W_E(z)} = C_j , z ∈∂D_j , j∈{1,2,…, N_I}, where $C_{j}$ are real constants that are unknown a priori. In (3.1b), the zero Neumann condition on $V$ has been re-expressed as a constant condition on its harmonic conjugate.

Suppose now that $W_{E}$ , satisfying (3.1), is known. Then the complex electric field is given by $E\equiv E_{x}+\mathrm{i}E_{y}=-\overline{dW_{E}/dz}$ , from which the current may be obtained via Ohm’s law, $J=-\sigma\overline{dW_{E}/dz}$ . The cross product giving the Lorentz force is performed through a pre-multiplication by $-\mathrm{i}$ , which corresponds to a rotation through an angle of $-\pi/2$ in the plane, giving $F=\mathrm{i}\sigma B_{0}\overline{dW_{E}/dz}$ . The Lorentz potential and force are then written as {subeqnarray} W_L & = -iσB_0 W_E
F = ¯dW_L/dz. In complex notation, the two-dimensional velocity $\tilde{\boldsymbol{u}}$ may thus be written as

u=\overline{\frac{d}{dz}\left(W_{L}-W_{P}\right)}\equiv\overline{\frac{d}{dz}\left(W_{\mathrm{flow}}\right)},

(2)

where $W_{P}$ is an analytic function satisfying $\mathrm{Re}\left\{W_{P}\right\}=P(x,y)$ . In (2), we have defined the flow potential function $W_{\mathrm{flow}}=W_{L}-W_{P}$ .

The precise mechanism for the generation of circulation is now evident and can be described as follows. First, the electrostatic problem (3.1) is solved in the fluid domain exterior to the collection of insulators and conductors held at different electrical potentials. Depending on the geometry and the applied voltages, some amount of surface charge, $Q_{j}$ , accumulates on each conductor surface $\partial B_{j}$ in the electrostatic problem. Each such surface charge manifests as a source term in the complex potential $W_{E}$ so that the potential satisfies the relation,

W_{E}=-\sum_{j=1}^{N_{C}}\frac{Q_{j}}{2\pi}\log{\left(z-z_{j}\right)}+\mathrm{single\;valued\;function},

(3)

where $z_{j}\in B_{j}$ . The amount of current, $I_{j}$ , leaving each conductor, $B_{j}$ , is then trivially related to the induced charge in the electrostatic problem through Ohm’s law so that $I_{j}=\sigma Q_{j}$ .

The Lorentz force potential is obtained through (3.1a), which converts the source term present in (3) into a circulation term through a pre-multiplication by $\mathrm{i}$ . Since the pressure is a single-valued function, we immediately deduce that

W_{\mathrm{flow}}=-\mathrm{i}\frac{B_{0}}{2\pi}\sum_{j=1}^{N_{C}}I_{j}\log{\left(z-z_{j}\right)}+W_{\mathrm{flow}}^{\mathrm{SV}},

(4)

where $W_{\mathrm{flow}}^{\mathrm{SV}}$ is a single-valued function. Thus, the amount of electrical current leaving a particular conductor alone dictates the induced circulation around that same conductor in the induced fluid flow, $\Gamma_{j}=B_{0}I_{j}$ .

After obtaining $Q_{i}$ through the solution to the electrostatic problem (3.1), the boundary value problem for $W_{\mathrm{flow}}$ becomes fully specified by (4) supplemented by impermeability conditions on each conducting body which may be expressed as follows,

\mathrm{Im}\left\{W_{\mathrm{flow}}^{\mathrm{SV}}(z)\right\}=\mathrm{Im}\left\{\mathrm{i}\frac{B_{0}}{2\pi}\sum_{j=1}^{N_{C}}I_{j}\log{\left(z-z_{j}\right)}\right\}+\psi_{k},

(5)

where $z\in\partial B_{k}$ for $k\in\{1,2,...,N_{C}\}$ , and $\psi_{k}$ are real constants that are unknown apriori. We have assumed the absence of additional driving forces. However, we note that a background free stream of complex velocity $U$ , for example, may be included simply by augmenting the right side of (4) with the term $\overline{U}z$ and the right side of (5) with the term $-\mathrm{Im}\left\{\overline{U}z\right\}$ .

Note that even in the case of a strictly magnetically driven flow, a non-zero pressure field ( $W_{P}\neq 0$ ) may be required to enforce the impermeability boundary conditions, to be explained as follows. In the case that all obstacles in the flow are perfect conductors ( $N_{I}=0$ ), it is clear that $W_{\mathrm{flow}}=W_{L}=-\mathrm{i}\sigma B_{0}W_{E}$ alone satisfies the fluid flow boundary conditions and $W_{P}=0$ . Since the electric field lines of $W_{E}$ are necessarily perpendicular to the boundary of a perfect conductor, the multiplication by $\mathrm{i}$ in (3.1a) ensures that the fluid flow streamlines are tangent to each conducting surface and thus $W_{L}$ is the valid fluid flow potential .

However, in cases where some obstacles are electrical insulators, $W_{L}$ alone does not satisfy the impermeability boundary conditions. On the surface of each insulator, the electric field lines of $W_{E}$ are necessarily tangent to the insulating surface. Hence, the multiplication by $\mathrm{i}$ in (3.1a) produces fluid flow streamlines that are normal to each insulating surface, so that $W_{L}$ alone violates the impermeability condition on insulators. Hence, when insulators exist in the flow, a non-zero pressure field develops to enforce the impermeability condition on insulators. When electrically insulating obstacles are present, the fluid flow must thus be obtained in two steps. First, one must solve the electrostatic problem and thus obtain a set of charges $Q_{j}$ accumulated on the surface of each conducting body. Second, those charges are used to specify the circulation in the fluid flow problem so that the flow problem described by (4)-(5) is fully posed and can be solved. We note again that since the circulation around each body is specified, the problem is well-posed and we do not encounter the issue of undetermined circulations that plagues high Reynolds number aerofoil theory (gonzalez2022variational). We proceed in §3.2 with a brief discussion relating the electrostatic problem (3.1) to the induced circulations. We subsequently derive mathematical solutions for the flow in various multiply-connected geometries.

3.2 Connection Between Induced Circulation Electrostatic Capacitance

It is worth noting that in the special case where only two distinct voltage values are prescribed on the boundary, (3.1) is an electrostatic capacitance problem. The electrostatic capacitance is the measure of the magnitude of induced charge on each conducting surface per unit applied voltage difference, and it is determined by the geometry. It is a conformal invariant and is the reciprocal of the extremal distance (ahlfors2010conformal, pp. 65). A significant literature exists regarding bounds and comparison theorems for the capacitance of different geometries (ahlfors2010conformal, pp. 65). Recall that in §3.1, it was shown that the circulation of the induced fluid flow is solely determined by the charges present in the electrostatic problem (3.1). Thus, theorems which indicate how one might geometrically tune the capacitance of the electrostatic problem also indicate how one might tune the circulation of the induced fluid flow. For example, the introduction of an insulator anywhere into a given probe configuration necessarily decreases the capacitance and hence also the induced circulation. Meanwhile, the introduction of a floating conductor necessarily increases the capacitance and circulation.

3.3 Magnetically-Driven Flows Around a Collection of Conductors ( $W_{P}=0$ )

Through example, we now outline a procedure for obtaining solutions when there are no electrically insulating obstacles in the flow ( $N_{I}=0$ ). We consider the geometry in figure 2(b), which is a special case of the general geometry of figure 1(b) where $N_{I}=0$ and $N_{C}=2$ .

Two conducting obstacles lie in a two-dimensional conducting fluid. The first is held at the voltage $V_{1}=1$ while the second is held at $V_{2}=0$ . Otherwise the flow is infinite in extent and there are no external pressures driving fluid flow. In order to determine the induced fluid flow, we first must solve (3.1) for the electrical potential, and then the resulting flow is given by (3.1a) and (2) after taking $W_{P}=0$ , as was outlined in §3.1, since all the immersed bodies are perfect conductors. We solve for $W_{E}$ by first conformally mapping the physical domain in figure 2(b) onto the concentric annulus of figure 2(a).

The conformal map from the physical geometry in figure 2(b) to the an annulus of figure 2(a) is given by,

\zeta(z)=\sqrt{\rho}\frac{\sqrt{d^{2}-s^{2}}-z}{\sqrt{d^{2}-s^{2}}+z},

(6)

where the inner annular radius is given by $\rho=\left(\left(d-\sqrt{d^{2}-s^{2}}\right)/s\right)^{2}$ . Note that for unit circles, $s=1$ . The left conductor in the physical domain is mapped to the inner circle $|\zeta|=\rho$ while the right conductor is mapped the the outer circle of the annulus, $|\zeta|=1$ . It is simple to solve the Dirichlet problem in the conformally mapped domain. In order to achieve $V=1$ on $|\zeta|=\rho$ and $V=0$ on $|\zeta|=1$ , the potential in the annulus must be

W_{E,\zeta}(\zeta)=V_{0}\frac{\log{\left(\zeta/\rho\right)}}{\log{\left(1/\rho\right)}}.

(7)

The complex potential in the physical domain is then given simply by

W_{E}(z)=W_{E,\zeta}(\zeta(z))=V_{0}\frac{\log{\left(\frac{1}{\sqrt{\rho}}\frac{\sqrt{d^{2}-s^{2}}-z}{\sqrt{d^{2}-s^{2}}+z}\right)}}{\log{\left(1/\rho\right)}}.

(8)

The complex potential for the flow is then obtained simply through multiplication by $-\mathrm{i}\sigma B_{0}$ , as follows

W_{\mathrm{flow}}(z)=-\mathrm{i}\sigma B_{0}V_{0}\frac{\log{\left(\frac{1}{\sqrt{\rho}}\frac{\sqrt{d^{2}-s^{2}}-z}{\sqrt{d^{2}-s^{2}}+z}\right)}}{\log{\left(1/\rho\right)}},

(9)

It is clear from (9) that the applied voltage difference manifests as a fluid flow circulation of magnitude $\Gamma=2\pi\sigma B_{0}V_{0}$ around each of the cylinder conductors. The induced fluid flow is plotted in figure 3(a).

3.4 Uniform Stream Past Perfect Conductors

We now derive the solution for a strictly pressure-driven flow — which possesses zero circulation — past the same two cylinders. We proceed by demonstrating, through direct calculations, how magnetic fields may induce circulation in the pressure-driven flow.

3.4.1 The Pressure-Driven Uniform Stream (No Circulation)

Consider again the circular obstacle geometry given in figure 2(b), except now a uniform steam $U\in\mathbb{R}$ is directed along the real axis past the two cylinders with zero circulation around each cylinder. Such a flow might be generated, for example, by an applied pressure gradient along the $\hat{\boldsymbol{x}}$ direction.

In either case, the flow solution can be obtained directly by the methods of crowdy2020solving, who presents a calculus for potential theory. With vanishing circulations around each body, and the behaviour at infinity specified ( $W_{{}_{\mathrm{flow}}}\sim Uz$ as $|z|\rightarrow\infty$ ), a class of mathematical solutions for the magnetohydrodynamically driven flow in the Hele-Shaw cell become expressible in terms of the prime function. Notably, the prime function has a closed-form series reprentation in the doubly-connected case baddoo2020exact. In what follows, we present the solution procedure for the two-cylinder problem. We later outline how N-body solutions follow by an identical procedure (for more details, see crowdy2020solving).

Consider two unit cylinders a distance $2d$ apart, with each centered on the real axis, as in figure 2(b). Take the left-most circle to be centered on the origin. The Möbius transformation $\zeta=1/z$ then maps the origin-centered circle to itself. However, the second circle is mapped to a circle of radius $\delta<1$ inside the unit circle and centered at the point $q$ . More generally, if additional cylinders were added to the physical domain, the same Möbius transformation would take each additional cylinder to a distinct excised circle contained inside the unit disk. The unit circle, with a number of excised interior circles, represents a canonical domain, where the so-called prime function $\omega(z_{1},z_{2})$ becomes useful for obtaining mathematical solutions to certain boundary value problems (crowdy2020solving). Let us return to the two-cylinder problem of current interest. In this case (figure 2(b)), it is convenient to choose the canonical domain to be a concentric annulus ( $q=0$ ) as drawn in figure 2(a), since the solution possesses an explicit sum representation there.

The particular map to the concentric annulus was already presented in (6). By direct manipulation, it may be shown that the map can be re-expressed as $z=(1-|a|^{2})/(\zeta-a)-\overline{a}-d$ , where $a=-d+\sqrt{d^{2}-1}$ . Recasting the map in this form makes it clear that $\zeta=a$ is the pre-image of infinity. We now seek a complex potential which, in the physical domain, possesses zero circulation around each body and which tends to $W_{\mathrm{flow}}\sim Uz$ as $|z|\rightarrow\infty$ .

The free-stream condition in the $\zeta$ -plane becomes $W_{\mathrm{flow},\zeta}\sim U(1-|a|^{2})/(\zeta-a)$ as $\zeta\rightarrow a$ . The impermeability condition on rigid boundaries may be expressed as $\mathrm{Im}\left\{W_{\mathrm{flow},\zeta}(\zeta)\right\}=\tilde{C}_{j}$ for $z\in\partial B_{j}$ for some set of constants $\{\tilde{C}_{j}\}$ . A function which possesses all of the necessary properties, in the canonical domain, can be written concisely as follows, {subeqnarray} W_flow,ζ(ζ) & = U(1-—a—^2)ϕ_0(ζ,a),
ϕ_0(ζ,a) = -1aK(ζ,a)+1aK(ζ,1¯a). where $\mathcal{K}\left(\zeta,a\right)=a\partial\log{\left(\omega(z,a)\right)}/\partial a$ and $\omega$ is the prime function as developed in crowdy2020solving. The function $\phi_{0}$ has the two important properties: it has a simple pole with unit residue at $\zeta=a$ , and it maps the circles of the canonical domain to slits parallel to the real axis, which possess a constant imaginary part (see crowdy2020solving). Hence, $\phi_{0}$ satisfies the impermeability condition on each body. The $\mathcal{K}$ functions in (3.4.1) possess a simple series representation in the doubly-connected case (crowdy2020solving, pp. 280), which was used to generate the plot of flow streamlines given in figure 3(b). Note that the solution in the physical domain is simply given by $W_{{}_{\mathrm{flow}},\zeta}{\left(\zeta(z)\right)}$ , which is computed by substituting (6) into (3.4.1).

We note that if instead there were $N>2$ cylinders in the flow, the same solution (3.4.1) applies. However the definition of $\mathcal{K}$ changes. In such a case, a Möbius map converts the physical domain into a canonical domain comprising a unit disk with $N-1$ excised disks. The definition of $\mathcal{K}$ depends on the form of this canonical domain through $\omega$ and, in higher connectivities, it must be computed numerically by methods described in crowdy2007computing and crowdy2020solving.

3.4.2 Magnetic Driving Induces Circulation

In this section, we demonstrate how magnetic driving may be used to induce circulation into an otherwise circulation-free pressure driven flows such as that described in §3.4.1 and illustrated in figure 3(b). The solution for a flow comprising a pressure-driven free stream in addition to magnetic driving can be obtained simply via superposition of the solutions from §3.3 and §3.4.1.

In figure 3(b), streamlines for the flow of a uniform stream in a Hele-Shaw cell are plotted using (3.4.1). Meanwhile, in figure 3(a), the flow streamlines of the magnetic-driven flow from §3.3 are plotted. Figure 3(c) then plots flow streamlines when a pressure-driven uniform stream is combined with a circulatory magnetic flow, of the type described in §3.3. Figure 3(d) shows the streamlines when the intensity of magnetic flow is increased relative to the situation in 3(c); therein, the relative voltage between the conducting probes is increased by a factor of five which increases the magnitude of circulation around each body. The resulting flow streamlines in figure 3(d) are successfully diverted between the cylinders.

Increasing the relative voltage between probes increases the magnitude of induced circulations, leading to a stronger diversion of the fluid flow from the uniform stream profile. Alternatively, the circulation can be enhanced for a fixed voltage differential by adjusting the geometry. As seen in (7), the electrostatic capacitance of the two-cylinder system is $C=2\pi/\log{\left(1/\rho\right)}$ . Since the capacitance is a conformal invariant, any configuration of conductors that can be conformally mapped to an annulus with inner radius $\rho$ and outer radius unity will induce the same circulation around each of the two conductors in the flow (for a fixed voltage difference). Thus, by changing the geometry in a conformally inequivalent manner, the induced circulation may thus be changed without modifying the applied voltage difference. See Leviconds for a related discussion in the electrostatic context.

3.4.3 General $N$ -body Solutions Using Framework of Crowdy

Up to this point, we have examined the two-body problem to illustrate the physical mechanism for circulation generation.

Any two-body problem can be conformally mapped onto a concentric annulus having some inner radius $\rho$ and an outer radius $1$ , where the value of $\rho$ depends on the geometry. In contrast to the Riemann mapping theorem — applicable for simply-connected domains — not all doubly-connected domains are necessarily conformally equivalent (they are equivalent only if $\rho$ happens to be the same). More generally, any $N$ -body problems can be mapped to the interior of a unit circle with $N-1$ excised circles whose positions and radii are determined by the physical geometry. The concentric annulus is the special case of exactly one excised circle. To obtain the fluid flow around a number of physical obstacles, such as $N$ cylinders, it suffices to solve a problem in a conformally equivalent canonical domain, and then to map back to the physical domain of interest. Exact solutions for problems in the canonical domain are accessible through the framework of crowdy2020solving. Note that to evaluate the solution in the physical domain, one must possess a conformal map to the physical domain of interest. In arbitrary geometries, the map between the physical and canonical domains may be difficult to attain analytically.

The procedure for solving the $N$ -body problem is as follows. First, one must find a mapping from the physical domain of interest to a canonical domain comprising the unit circle with $N-1$ excised disks. crowdy2020solving showed that this map can be written exactly in terms of the so-called prime function $\omega(z,a)$ , for a number of physically relevant domains. For the case of the unbounded domain exterior to finitely many cylinders of arbitrary position and radius, a simple Möbius transformation brings the domain onto a canonical domain, as was described in §3.4.1. Once the mapping to the canonical domain is found, the boundary value problem can be solved there instead. We now focus on solutions in canonical domains.

The solution in the canonical domain is amenable to the techniques developed by crowdy2020solving and is expressible exactly in terms of the prime function. In relation to §3.3, we seek a complex analytic function $W_{E}$ defined on a canonical domain which satisfies 3.1(a) on the unit circle and each excised disk. The geometry in figure 4(c) represents $N_{C}=6$ , so that there are five excised circles. In what follows we only consider only situations where all boundaries are perfectly conducting, so that $N_{I}=0$ , and where external voltages are applied to each conductor. We also allow the possibility that some of the conductors are floating and not explicitly connected to a voltage source. We now describe how one may obtain solutions for the induced fluid flow in such geometries.

crowdy2020solving introduced a special set of functions called integrals of the first kind, $\{v_{j}(z)\}$ , defined in the canonical domain in terms of the prime function $\omega$ , which will serve as the basis for the solutions obtained in the present section. crowdy2007computing and crowdy2020solving present efficient methods for computing the prime function. Codes for computing the prime function are readily available (crowdy2016schottky).

There exists one such function, $v_{j}(z)$ , for each of the excised circles so that $j\in\{1,2,...,N_{C}-1\}$ , each possessing the following two important properties. First, $v_{j}(z)$ introduces a unit circulation around the $j^{\mathrm{th}}$ excised circle and exactly zero circulation around each of the other excised circles. Second, the imaginary part of $v_{j}(z)$ is constant on all other boundaries. It is thus clear that a linear superposition of the functions $\mathrm{i}v_{j}(z)$ is capable of meeting the boundary conditions in 3.1(a); the coefficients of the superposition are determined through the solution of a linear system of $N_{C}$ equations which is easily solved with the backslash operator in MATLAB. Details are presented in Appendix LABEL:vjsecapp. Once the electrostatic potential has been obtained, pre-multiplication by $-\mathrm{i}\sigma B_{0}$ gives the complex potential of the fluid flow as was described in §3.1.

Solutions in two canonical domains, as obtained through the solution of a linear system of equations (see Appendix LABEL:vjsecapp), are presented in figures 4(b) and 4(c). Figure 5(a) shows a geometry that was explored experimentally by david2023magnetostokes, wherein the authors were unable to attain an analytic solution. Figure 4(b) shows the streamlines of the exact solution for the fluid flow in the same geometry as obtained through the method of the present section. The central cylinder is held at $V=1$ and the outer cylinder at $V=0$ . The off-centre cylinder is taken to as a floating conductor, since it was not connected to voltage source in the experiments. As a consequence, there is no net circulation around the floating conductor. In figure 5(c), we plot flow streamlines in a domain of higher connectivity where $N_{C}=6$ . The central cylinder is held at $V=1$ and the outer cylinder at $V=0$ . The other conductors are chosen to be floating. Note that any combination of the conductors can be taken to be floating or possessing a prescribed voltage.

3.5 Experiment

Here we outline a simple experiment which realizes one of the exact solutions from §3.4.3. This experiment is meant to motivate and supplement the present theoretical paper but we note that it does not constitute a completely general experimental investigation of the magnetohydrodynamic Hele-Shaw cell. Future studies may focus on a more detailed experimental investigation of the system.

A Hele-Shaw cell was constructed from two thin circular sheets of transparent acrylic in a configuration similar to figure 1(a). Each transparent sheet had a diameter of $8\mathrm{cm}$ and a thickness of $1.5\mathrm{mm}$ . Two circles, with diameters of $d_{1}=7\mathrm{mm}$ and $d_{2}=10\mathrm{mm}$ , were cut from each disk. Then, six small holes of diameter $1\mathrm{mm}$ were cut in the top acrylic sheet along the line passing through the two circles of radii $d_{1}$ and $d_{2}$ , for the purposes of streamline visualization. A small drop of blue dye was place atop each such hole just before the experiment began. All cuts were made using a laser cutter.

Note that four additional circular holes can be seen in seen in figure 5(a) near the boundary of each acrylic sheet; these holes serve no function in the present experiment. However, near the top left such hole, an unintentional bubble appeared when filling the cell, which will be addressed below.

The two sheets were then separated by spacers of thickness $h=0.7\mathrm{mm}$ . The sheets were glued in place at each spacer, at several points along the boundary. Two metal cylinders, with diameters $d_{1}$ and $d_{2}$ , were then fitted into the cell. The cylinders fit snugly into the bottom acrylic sheet in such a way that no glue was necessary to prevent leakage. In the upper acrylic sheet, the holes for the cylinders were made slightly larger than the cylinder diameters to allow gas to escape in the event of electrolysis.

The entire Hele-Shaw cell was then placed atop a DZ08-N52 Neodymium Disc Magnet, as ordered from K&J Magnetics, with nominal surface field strength of $2340\mathrm{G}$ . The cell was then completely filled with saltwater using a needle. The saltwater was produced by simply adding salt packets to tap water. A small drop of blue dye was then placed atop each of the six $1\mathrm{mm}$ diameter holes just before the experiment began. A voltage difference of $1\mathrm{V}$ was then applied between the two cylinders. As the flow developed, the dye traced out streamlines as can be seen in figure 5(a).

Streamlines of exact solutions, as obtained though the procedure of §3.4.3, are plotted in figure 5(b). Only the six thoreotical streamlines which intersect the dye release points are plotted. The unintended bubble (top left) is treated as an impermeable boundary in the theoretical solution as is drawn in figure 5(b). There is a reasonable agreement between the experimental and theoretical streamlines.

Streamlines around the right-most cylinder are diffuse because the first few droplets of dye were dropped into place from too high and consequently spread beyond the small $1\mathrm{mm}$ opening. At the dye release points near the left cylinder, drops were added more carefully by gently touching a droplet onto the small opening, giving rise to more defined streamlines.

Note that the solution in figure 5(b) treats the bubbles as floating electrical conductors since the framework of §3.4.3 only allows for conducting obstacles. In reality, a bubble is better approximated by an insulator. However, since there are only two electrodes with fixed applied voltages in the experimental geometry (figure 5(a)), the flow streamlines are unaffected by this assumption. That is to say, the streamlines in figure 5(b) would be identical to those found if bubbles were treated as perfect insulators. However, the magnitude of the induced flow velocities does in fact depend on the type of boundary condition imposed: the induced circulation is higher in the presence of floating conductors compared to insulators. Generally, when there are multiple applied voltages, one must apply the appropriate insulating boundary condition on insulating obstacles in order to obtain accurate flow streamlines. Even in the case of figure 5(a), the proper insulating boundary condition needs to be applied in order to obtain accurate values for the velocity magnitudes at each point in the flow. A numerical scheme which allows as well for the presence of insulating obstacles, which thus overcomes the limitations of §3.4.3, is presented in §4.

4 Series Solutions

When no conformal mapping is known between the physical domain of interest and a canonical domain, or if insulating bodies exist in the flow, solutions are not attainable by the method presented in §3.4.3. Moreover, even when solutions can be written explicitly in terms of the prime function, as in §3.4.3, the numerical computation of the prime function can be expensive and actually involves the numerical solution of a different boundary value problem (crowdy2016schottky).

All of these facts lead us to seek out accurate and efficient numerical solutions for the fluid flow, in general settings. In this section, we demonstrate how series solutions can be adapted to our problem to solve for the fluid flow with many digits of accuracy in just a few seconds of computing time on a standard laptop (trefethen2018series). The procedure is described as follows.

The complex potential described by (3.1) is first expressed as a sum of Laurent series centered inside each body (exterior to the flow domain). The Laurent series are then truncated and their coefficients are determined through a least squares problem which enforces the conditions (3.1) at a finite number of points on each boundary.

The Vandermonde matrix in the least squares problem becomes exponentially ill-conditioned in the degree of approximation. As a result, the numerical solution eventually saturates with an increasing degree of approximation. brubeck2021vandermonde showed that by performing Arnoldi orthogonalization, one can improve the condition number and thus achieve more digits of accuracy in the numerical solution. In some least squares problems, the difference in accuracy between solutions obtained with and without Arnoldi can be quite dramatic, even reaching ten digits (brubeck2021vandermonde, see figure 3.1). Note that for the examples considered in the present section, the Arnoldi orthogonalization is not essential for attaining quite accurate solutions. In more general geometries, the Arnoldi procedure may be necessary to achieve high accuracy solutions.

The numerical method described herein is similar to that presented by baddoo2020lightning, but with modified boundary conditions and an extension to multiply-connected domains. Note that when sharp corners are present in the solution domain, strong singularities in the solutions emerge, and rapidly converging lightning solvers can be applied (gopal2019solving). Lightning solvers supplement the Laurent series representation with a set of poles clustered near sharp corners to approximate singularities. In the present work, we focus on smooth bodies which do not require the placement of such poles. However, the procedure for including such poles is rather straightforward (baddoo2020lightning; gopal2019solving).

4.1 Numerical Solution of Electrostatic Problem

We seek the electrostatic potential in the unbounded domian exterior to $N=N_{I}+N_{C}$ distinct bodies, with $N_{C}$ conducting surfaces $\{\partial B_{j}\}$ held generally at different voltages $V_{j}$ , and $N_{I}$ insulating surfaces $\{\partial D_{j}\}$ obeying zero Neumann conditions $\nabla V_{j}\boldsymbol{\cdot}\boldsymbol{n}=0$ . For notational convenience, let $z_{j}$ be a point interior to the $j^{\mathrm{th}}$ conductor for $j\in\{1,...,N_{C}\}$ and the interior of the $\left(j-N_{C}\right)^{\mathrm{th}}$ insulator for $j\in\{N_{C}+1,...,N_{C}+N_{I}\}$ .

The electrical potential then takes the form (3), where $Q_{i}$ are undetermined constants. The solution can be written generally as a sum of a Laurent series and its logarithm terms as {subeqnarray} W_E(z) & = a_0 + ∑_j=1^N_C Q_j log