Magnetohydrodynamic flow control in Hele-Shaw cells

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. Physically, the system is particularly interesting because it allows for circulation in the flow, which is not possible in the standard pressure-driven Hele-Shaw cell. We first elucidate the mechanism of flow generation both physically and mathematically. After formulating the problem using complex variables, we present mathematical solutions for a class of canonical multiply connected geometries in terms of the prime function framework developed by Crowdy (Solving Problems in Multiply Connected Domains, SIAM, 2020). We then demonstrate how recently developed fast numerical methods may be applied to accurately determine the flow field in arbitrary geometries.

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 (Moffatt 1978; Davidson 2001). 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 (Alfvén 1942) and geodynamos (Moffatt 1978; Moffatt & Dormy 2019). 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 Gonzalez & Taha (2022) for a recent development and survey of the aerofoil problem). Because 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.

Moffatt (1991) 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, Mirzadeh et al. (2020) 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.

Bau, Zhong & Yi (2001) and Zhong, Yi & Bau (2002) 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 annulus. When the devices were placed in a uniform magnetic field, Lorentz forces induced a fluid flow in the bulk. Yi, Qian & Bau (2002) showed that, under certain conditions, the application of periodic AC electrical currents in similar devices may lead to chaotic advection. Later, Homsy et al. (2005) manufactured a high current density DC microfluidic pump which found applications in nuclear magnetic resonance (Homsy et al. 2007); notably, the pump avoided significant Joule heating. A recent comprehensive review of applications of magnetohydrodynamics in the context of microfluidics is given by Bau (2022).

The papers discussed in the paragraph above were mainly experimental in nature. In the cases where fluid flows were analysed analytically, such as in Bau et al. (2001) and Zhong et al. (2002), 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.

More recently, David et al. (2023) investigated magnetically driven flow in a Taylor–Couette geometry (concentric annulus) with a free surface. Experiments presented therein beautifully reproduced the classic kinematic reversibility experiments of Taylor (1967) without the need for moving boundaries. After investigating mixing in this flow geometry, the authors pointed out that, under certain assumptions, and in the limit of a shallow fluid layer, the flow is approximated by potential flow.

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 impermeable 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 Crowdy (2020). In general geometries where the relevant conformal mappings needed to apply Crowdy's prime function machinery are difficult to find, we show how highly accurate approximate solutions can be obtained using series solution methods as described by Trefethen (2018). The principal aims of the present paper are to elucidate the physical mechanism of circulation generation, to develop a mathematical formulation of magnetically driven Hele-Shaw flow, to provide a framework for solving the mathematical problem and to demonstrate how such flows can be used to augment known pressure-driven flows.

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, $\boldsymbol {\nabla } V \boldsymbol {\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.

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 Crowdy (2020). One such solution gives the flow in a geometry explored experimentally by David et al. (2023). In § 3.5, another solution is compared with 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.

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 electrolyte 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 (Zaltzman & Rubinstein 2007). 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 laboratory frame is given by the gradient of a harmonic potential, $\boldsymbol {E}_{{lab}}=-\boldsymbol {\nabla }V$. The electric field drives an electrical current density 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 laboratory frame by the relation $\boldsymbol {E}=\boldsymbol {E}_{{lab}}+\boldsymbol {u}\times \boldsymbol {B}$, owing to the fact that the electric field is not invariant under Galilean transformations (see Moffatt 1978, p. 34). 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 $1\ {\rm V}$; voltages higher than roughly $1.23\ {\rm V}$ create bubbles due to water electrolysis. With typical separation distances between electrical probes being of the order of centimetres, the typical electric field magnitude is $E_0\approx 1\ {\rm V}\ \mathrm {cm}^{-1}$. Typical velocities in the cell are $U_0\approx 1\ \mathrm {mm}\ \mathrm {s}^{-1}$, 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 with $\boldsymbol {E}_{{lab}}$ scales as $U_0B_0/E_0=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

(2.1a)\begin{gather} \rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}\right) = \mu \nabla^2 \boldsymbol{u}- \boldsymbol{\nabla} P+\sigma \boldsymbol{E}\times \boldsymbol{B}, \end{gather}

(2.1b)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}

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 Batchelor 1967, p. 222). Hence, the left-hand side of (2.1a) 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 (Davidson 2001, figure 5.10 on p. 153); see also the work of Rossow (1958) 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 negligible 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

(2.2a)\begin{gather} \boldsymbol{u} = \frac{1}{2\mu}\left( -\boldsymbol{\nabla} P+\sigma B_0\boldsymbol{E}\times \hat{\boldsymbol{z}}\right)z(h-z), \end{gather}

(2.2b)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}

To make further progress, we note that the Lorentz force, $\boldsymbol {F}=-\sigma B_0 \boldsymbol {\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}=\boldsymbol {\nabla }\phi$, where $\nabla ^2\phi =0$. In general $\phi$ need not be single valued. Equation (2.2a) then reduces to

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

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, (2.2b) and (2.3) 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 \boldsymbol {\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}=\boldsymbol {\nabla } \phi$, may be obtained. Finally, a single-valued pressure field $P$, as appears in (2.3), 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 are prescribed so that $V=V_j$ for $j\in \{1,2,\ldots, N_C\}$, where $N_C$ is the total number of conductors. On insulating boundaries, $\partial D_k$, the normal derivative must vanish so that $\boldsymbol {\nabla } V \boldsymbol {\cdot } \boldsymbol {n}=0$, where $k\in \{1,2,\ldots, 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+{\rm i}y$. 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

(3.1a)\begin{gather} \mathrm{Re}\left\{W_E(z)\right\} = V_j,\quad z \in \partial B_j,\ j\in\{1,2,\ldots, N_C\}, \end{gather}

(3.1b)\begin{gather}\mathrm{Im}\left\{ W_E(z)\right\} = C_j,\quad z \in \partial D_j,\ j\in\{1,2,\ldots, N_I\}, \end{gather}

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 {{\rm d}W_E/{\rm d}z}$, from which the current may be obtained via Ohm's law, $J=-\sigma \overline {{\rm d}W_E/{\rm d}z}$. Here, the complex current density is defined as $J=J_x+\mathrm {i}J_y$, where $\boldsymbol {J}=(J_x,J_y)$. 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 $-{\rm \pi} /2$ in the plane, giving $F=\mathrm {i}\sigma B_0\overline {{\rm d}W_E/{\rm d}z}$. The Lorentz potential and force are then written as

(3.2a)\begin{gather} W_L ={-}\mathrm{i}\sigma B_0 W_E \end{gather}

(3.2b)\begin{gather}F = \overline{{\rm d}W_L/{\rm d}z}. \end{gather}

In complex notation, the two-dimensional velocity $\tilde {\boldsymbol {u}}$ may thus be written as

(3.3)\begin{equation} u= \overline{\frac{{\rm d}}{{\rm d}z}\left(W_L-W_P\right) }\equiv \overline{\frac{{\rm d}}{{\rm d}z}\left(W_{{flow}}\right) }, \end{equation}

where $W_P$ is an analytic function satisfying $\mathrm {Re}\{W_P\}=P(x,y)$. In (3.3), we have defined the flow potential function $W_{{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. We reiterate that while the mathematical problem (3.1) describes an electrostatic problem, this is just a mathematical step leveraged to attain the steady electrical current flow, wherein charges do in fact flow. 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. We now solve the electrical problem with the electrical permittivity equal to unity, $\epsilon =1$, since the final solution will clearly be independent of $\epsilon$. Each surface charge manifests as a source term in the complex potential $W_E$ so that the potential satisfies the relation

(3.4)\begin{equation} W_E ={-}\sum_{j=1}^{N_C}\frac{Q_j}{2{\rm \pi}} \log{\left(z-z_j\right)}+\mathrm{single\ valued\ function}, \end{equation}

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.2a), which converts the source term present in (3.4) into a circulation term through a pre-multiplication by $\mathrm {i}$. Since the pressure is a single-valued function, we immediately deduce that

(3.5)\begin{equation} W_{{flow}}={-}\mathrm{i}\frac{ B_0}{2{\rm \pi}}\sum_{j=1}^{N_C} I_j \log{\left(z-z_j\right)}+W_{{flow}}^{{SV}}, \end{equation}

where $W_{{flow}}^{{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, $\varGamma _j=B_0 I_j$.

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

(3.6)\begin{equation} \mathrm{Im}\left\{W_{{flow}}^{{SV}}(z)\right\} = \mathrm{Im}\left\{\mathrm{i} \frac{ B_0}{2{\rm \pi}}\sum_{j=1}^{N_C} I_j \log{\left(z-z_j\right)}\right\}+\psi_k, \end{equation}

where $z \in \partial B_k$ for $k\in \{1,2,\ldots, N_C\}$, and $\psi _k$ are real constants that are unknown a priori. 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-hand side of (3.5) with the term $\bar {U}z$ and the right-hand side of (3.6) with the term $-\mathrm {Im}\{\bar {U}z\}$.

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_{{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.2a) 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.2a) 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 (3.5)–(3.6) 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 (Gonzalez & Taha 2022). 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 (Ahlfors 2010, p. 65). A significant literature exists regarding bounds and comparison theorems for the capacitance of different geometries (Ahlfors 2010, p. 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 that 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.2a) and (3.3) 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).

Figure 2. (a) Conformally mapped domain of the physical geometry given in panel (b). The circles $|\zeta |=1$ and $|\zeta |=\rho$ map to the two conductor boundaries in the physical plane. (b) Schematic for electrostatic problem between two perfectly conducting cylinders held at fixed voltages in the physical $z$-domain.

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

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

where the inner annular radius is given by $\rho =( ( d-\sqrt {d^2-s^2} )/s )^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 to 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

(3.8)\begin{equation} W_{E,\zeta}(\zeta)=V_0\frac{\log{\left(\zeta/\rho\right)}}{\log{\left(1/\rho\right)}}. \end{equation}

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

(3.9)\begin{equation} W_E(z)=W_{E,\zeta}(\zeta(z))=V_0 \frac{\log{\left(\dfrac{1}{\sqrt{\rho}}\dfrac{\sqrt{d^2-s^2}-z}{\sqrt{d^2-s^2}+z}\right)}}{\log{\left(1/\rho\right)}}. \end{equation}

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

(3.10)\begin{equation} W_{{flow}}(z)={-}\mathrm{i} \sigma B_0V_0 \frac{\log{\left(\dfrac{1}{\sqrt{\rho}}\dfrac{\sqrt{d^2-s^2}-z}{\sqrt{d^2-s^2}+z}\right)}}{\log{\left(1/\rho\right)}}. \end{equation}

It is clear from (3.10) that the applied voltage difference manifests as a fluid flow circulation of magnitude $\varGamma = 2{\rm \pi} \sigma B_0 V_0$ around each of the cylinder conductors. The induced fluid flow is plotted in figure 3(a).

Figure 3. (a) Visualization of the flow solution given in (3.10). Grey lines represent electric field lines and the direction of current flow, $\boldsymbol {J}$. Black lines are streamlines of the fluid flow. A circulation of magnitude $\varGamma =|\sigma B_0 V_0|$ is induced around each body. (b) Black lines represent streamlines of a uniform flow past two conducting bodies. The given flow may be driven by either pressure or an external electric field since both solutions are identical. (c) Streamlines of a uniform flow past two cylinders when a potential difference $\Delta V =1$ is applied between the two cylinders, which drives electrical current and hence a circulation around each body. (d) Same plot as panel (c) with a larger applied voltage differential, and hence larger electrical current, between the cylinders. The larger current induces more circulation around each body as compared with panel (c).

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. Note that the conformal maps being used to solve the Laplace problems in the present section follow directly from theoretical developments laid out in earlier work summarized in Crowdy (2020).

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.

The flow solution can be obtained directly by the methods of Crowdy (2020, chap. 15), who presents a calculus for potential theory. With vanishing circulations around each body, and the behaviour at infinity specified ($W_{{{flow}}}\sim U z$ 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 representation in the doubly connected case (Baddoo et al. 2020). 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 Crowdy (2020, p. 294) and Crowdy 2006b,a).

Consider two unit cylinders a distance $2d$ apart, with each centred on the real axis, as in figure 2(b). Take the left-most circle to be centred on the origin. The Möbius transformation $\zeta =1/z$ then maps the origin-centred circle to itself. However, the second circle is mapped to a circle of radius $\delta <1$ inside the unit circle and centred 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 (Crowdy 2020).

Let us return to the two-cylinder problem of current interest. In this case (figure 2b), 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 (3.7). By direct manipulation, it may be shown that the map can be re-expressed as $z=(1-|a|^2)/(\zeta -a)-\bar {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_{{flow}}\sim U z$ as $|z|\rightarrow \infty$.

The free-stream condition in the $\zeta$-plane becomes $W_{{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}\{ W_{{flow},\zeta }(\zeta ) \}=\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:

(3.11a)\begin{gather} W_{{flow},\zeta}(\zeta) = U\left(1-|a|^2\right)\phi_0(\zeta,a), \end{gather}

(3.11b)\begin{gather}\phi_0(\zeta,a) ={-}\frac{1}{a}\mathcal{K}\left(\zeta,a\right)+\frac{1}{a}\mathcal{K}\left(\zeta,\frac{1}{\bar{a}}\right), \end{gather}

where $\mathcal {K}(\zeta,a)=a\partial \log {(\omega (z,a))}/\partial a$ and $\omega$ is the prime function as developed in Crowdy (2020). 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 Crowdy 2020, p. 83). Hence, $\phi _0$ satisfies the impermeability condition on each body. The $\mathcal {K}$ functions in (3.11) possess a simple series representation in the doubly connected case (Crowdy 2020, p. 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_{{{flow}},\zeta }{(\zeta (z))}$, which is computed by substituting (3.7) into (3.11).

We note that if instead there were $N>2$ cylinders in the flow, the same solution (3.11) 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 Crowdy & Marshall (2007) and Crowdy (2020, chap. 14).

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.11). 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 (3.8), the electrostatic capacitance of the two-cylinder system is $C=2{\rm \pi} /\log {(1/\rho )}$. 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 Levi (2023) 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 Crowdy (2020). 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, although several general constructions exist including a generalized Schwarz–Christoffel formula for mapping to multiply connected polygonal regions (Crowdy 2007).

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. Crowdy (2020) 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 Crowdy (2020) 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.1a) 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.

Figure 4. (a) Experimental image from experiments performed by David et al. (2023) in a shallow free-surface annulus. (b) Exact solution, as obtained through the procedure of § 3.4.3, for the experimental geometry from panel (a). Flow streamlines are depicted in grey. The voltage distribution is indicated in colour. (c) Another mathematical solution, as obtained through the procedure of § 3.4.3. The central circle is held at $1\ \mathrm {V}$ while the outer circle is held at $0\ \mathrm {V}$. Other circles are taken to be floating conductors.

Crowdy (2020) 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. Crowdy & Marshall (2007) and Crowdy (2020, chap. 14) present efficient methods for computing the prime function. Codes for computing the prime function are readily available (Crowdy, Kropf & Green 2016).

There exists one such function, $v_j(z)$, for each of the excised circles so that $j\in \{1,2,\ldots,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.1a); 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 A.2. 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 A.2), are presented in figures 4(b) and 4(c). Figure 5(a) shows a geometry that was explored experimentally by David et al. (2023), wherein the authors did not seek an exact 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. Note that, in the case when insulating bodies are also present, a new generalized Schwarz integral formula due to Miyoshi & Crowdy (2023, § 2) can give explicit solutions to the electrostatic problem in terms of the prime function in some cases; however, we do not elaborate further on that method in the present paper. Secondary prime functions (Vasconcelos, Marshall & Crowdy 2015) can also be useful for dealing mixed boundary-value problems.

Figure 5. (a) Top-down view of the Hele-Shaw cell filled with saltwater. The cell sits atop a permanent neodymium magnet. Two aluminium cylinders, of diameters $7\ \mathrm {mm}$ and $10\ \mathrm{mm}$, completely penetrate the Hele-Shaw and are held at a voltage difference of $1\ \mathrm{V}$. (b) Same as panel (a) with mathematical solutions, as presented in § 3.4.3, plotted atop the experimental streamlines. The bubbles in the top left of the Hele-Shaw cell are modelled as impermeable obstacles in the exact solution.

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 5a), 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 with 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, 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 (Crowdy et al. 2016).

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 (Trefethen 2018). It is worth noting that such solvers are by no means the only method for solving such boundary-value problems as, for example, boundary integral methods are also applicable. However, we focus on these series methods due to their flexibility and simplicity in implementation. The solution procedure is described as follows.

The complex potential described by (3.1) is first expressed as a sum of Laurent series centred 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. Brubeck, Nakatsukasa & Trefethen (2021) 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 (Brubeck et al. 2021, 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 Baddoo (2020), 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 (Gopal & Trefethen 2019b). 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 (Gopal & Trefethen 2019b; Baddoo 2020).

4.1. Numerical solution of electrostatic problem

We seek the electrostatic potential in the unbounded domain 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 $\boldsymbol {\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,\ldots,N_C\}$ and the interior of the $(j-N_C){\mathrm {th}}$ insulator for $j\in \{N_C+1,\ldots,N_C+N_I\}$.

The electrical potential then takes the form (3.4), where $Q_i$ are undetermined constants. The solution can be written generally as a sum of a Laurent series and its logarithm terms as

(4.1a)\begin{gather} W_E(z) = a_0 + \sum_{j=1}^{N_C} Q_j \log{\left(z-z_j\right)}+\sum_{j=1}^N\sum_{k=1}^{\infty} \frac{a_k^{(j)}}{\left(z-z_j\right)^k}, \end{gather}

(4.1b)\begin{gather}\sum_{j=1}^{N_C} Q_j = 0, \end{gather}

(4.1c)\begin{gather}\mathrm{Re}\left\{W_E\right\} = V_j,\quad z\in \partial B_j,\ j\in \{1,2,\ldots,N_C\}, \end{gather}

(4.1d)\begin{gather}\mathrm{Im}\left\{W_E\right\} = C_j,\quad z\in \partial D_j,\ j\in \{1,2,\ldots,N_I\}, \end{gather}

where $C_j$ and $Q_j$ are a priori unknown real constants. Note that no logarithm terms are centred inside of electrical insulators. The condition (4.1b) is required to ensure that the potential is regular at infinity.

To convert (4.1) into a least-squares fitting problem, we must truncate each Laurent series. Let the expansion centred around the $j{\mathrm {th}}$ body be truncated at $N_L^{j}$ terms. Combining (4.1a) and (4.1b), the form of our approximation becomes

(4.2)\begin{equation} W_E(z) \approx a_0 + \sum_{j=2}^{N_C} Q_j \log{\left(\frac{z-z_j}{z-z_1}\right)}+\sum_{j=1}^N\sum_{k=1}^{N_L^{j}} \frac{a_k^{(j)}}{\left(z-z_j\right)^k}, \end{equation}

where we have implemented (4.1b) by setting $Q_1=-\sum _{j=2}^{N_C}Q_j$. Through this implementation, we guarantee that (4.1b) is satisfied exactly so that $W_E(z)$ is finite as $|z|\rightarrow \infty$. If (4.1b) were instead implemented as a constraint in the least-squares problem, then (4.1b) would only hold approximately and $W_E(z)$ would not be guaranteed to be finite as $|z|\rightarrow \infty$. We next choose $N_b^{j}$ sample points on the boundary of the $j{\mathrm {th}}$ body and fit $W_E$ to the appropriate boundary conditions, (4.1c)–(4.1d).

We now formulate the least-squares problem. Let $z_j^{(k)}$ denote the $j{\mathrm {th}}$ sample point on the surface of the $k{\mathrm {th}}$ body, where $k\in \{1,\ldots,N_I+N_C\}$ and $j\in \{1,\ldots,N_b^k\}$. In order to construct the Vandermonde matrix, we first define the vector

(4.3) \begin{align} \boldsymbol{V}_j^{(k)} &= \left[\left. (z_j^{(1)}-z_1)^{{-}1} \quad \cdots \quad (z_j^{(1)}-z_1)^{{-}N_L^{1}}\right| \cdots \left| (z_j^{(1)}-z_N)^{{-}1} \right.\right.\nonumber\\ &\left.\quad \cdots \quad (z_j^{(1)}-z_N)^{{-}N_L^{N}} \right] . \end{align}

We next define a vector of logarithm term evaluations corresponding to the $j{\mathrm {th}}$ sample point on the $k{\mathrm {th}}$ boundary as follows:

(4.4)\begin{equation} \boldsymbol{L}_j^{(k)} = \left[ \begin{array}{ccc} \log{\left(\dfrac{z_j-z_2}{z_j-z_1} \right)} & \ldots & \log{\left(\dfrac{z_j-z_{N_C}}{z_j-z_1} \right)} \end{array} \right] . \end{equation}

We also define the Laurent series coefficient vector as

(4.5)\begin{equation} \boldsymbol{a} = \left[\left. a_1^{(1)} \quad \cdots \quad a_{N_L^1}^{(1)}\right| \cdots \left| a_1^{(N)} \quad \cdots\quad a_{N_L^N}^{(N)} \right.\right] . \end{equation}

Now let $\boldsymbol{\mathsf{V}}^{(k)}$ be the matrix of containing all sample points on the $k{\mathrm {th}}$ body defined by

(4.6)\begin{equation} \boldsymbol{\mathsf{V}}^{(k)} = \left[ \begin{array}{c} \boldsymbol{\mathsf{V}}_1^{(k)}\\ \boldsymbol{\mathsf{V}}_2^{(k)}\\ \cdots\\ \boldsymbol{\mathsf{V}}_{N_{b}^{k}}^{(k)} \end{array} \right]. \end{equation}

Similarly, we let $\boldsymbol{\mathsf{L}}^{(k)}$ define the matrix of logarithm evaluations on the $k{\mathrm {th}}$ body

(4.7)\begin{equation} \boldsymbol{\mathsf{L}}^{(k)} = \left[ \begin{array}{c} \boldsymbol{\mathsf{L}}_1^{(k)}\\ \boldsymbol{\mathsf{L}}_2^{(k)}\\ \cdots\\ \boldsymbol{\mathsf{L}}_{N_{b}^{k}}^{(k)} \end{array} \right]. \end{equation}

Then, the least-squares problem for the undetermined coefficients can be expressed in terms of the matrix $\boldsymbol{\mathsf{A}}$ defined by

(4.8) \begin{equation} \boldsymbol{\mathsf{A}} = \left[ \begin{array}{cccc|c|cccc} \boldsymbol{0} & \boldsymbol{1} & -\mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(1)}\right\} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(1)}\right\} & \mathrm{Re}\left\{\boldsymbol{\mathsf{L}}^{(1)}\right\} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{1} & \vdots & \vdots & \vdots & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{1} & -\mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N_C)}\right\} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N_C)}\right\} & \mathrm{Re}\left\{\boldsymbol{\mathsf{L}}^{(N_C)}\right\} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0}\\ \boldsymbol{1} & \boldsymbol{0} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N_C+1)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N_C+1)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{L}}^{(N_C+1)}\right\} & -\boldsymbol{1} & \boldsymbol{0} & \cdots & \boldsymbol{0}\\ \boldsymbol{1} & \boldsymbol{0} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N_C+2)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N_C+2)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{L}}^{(N_C+2)}\right\} & \boldsymbol{0} & -\boldsymbol{1} & \boldsymbol{0} & \cdots\\ \boldsymbol{1} & \boldsymbol{0} & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ \boldsymbol{1} & \boldsymbol{0} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{L}}^{(N)}\right\} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & -\boldsymbol{1}\\ \end{array} \right], \end{equation}

along with the accompanying coefficient vector $\boldsymbol {c}$ defined by

(4.9)\begin{equation} \boldsymbol{c} = \left[ \mathrm{Im}\left\{a_0\right\} \quad \mathrm{Re}\left\{a_0\right\} \quad\mathrm{Im}\left\{\boldsymbol{a}\right\} \quad \mathrm{Re}\left\{\boldsymbol{a}\right\} \quad Q_2 \quad \dots \quad Q_{N_C} \quad C_1 \quad \dots \quad C_{N_I} \right].\end{equation}

Note that the vectors $\boldsymbol {1}$ and $\boldsymbol {0}$ in the definition of $\boldsymbol{\mathsf{A}}$ represent vectors of ones and zeros of the appropriate length; for example, the vector $\boldsymbol {1}$ in the top row has a length of $N_b^1$ while the $\boldsymbol {1}$ in the bottom row has a length of $N_b^N$. The load vector, having a length equal to the total number of boundary sample points, $\sum _{j=1}^N N_b^j$, then takes the form

(4.10)\begin{equation} \boldsymbol{f} = \left[ \left. V_1 \quad \cdots \quad V_1\right| \cdots \left| V_{N_C} \quad \cdots \quad V_{N_C}\right| 0 \quad \cdots \quad 0 \right], \end{equation}

with which the approximate form of (4.1) is expressed as the following matrix equation:

(4.11)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{c}^T \approx \boldsymbol{f}^T. \end{equation}

Equation (4.11) can be solved with the backslash operator in MATLAB for the coefficient vector $\boldsymbol {c}$, from which the complex potential can be readily evaluated using (4.2). Generally, the solution accuracy can be greatly improved by using Arnoldi orthogonalization on the Vandermonde part of the matrix $\boldsymbol{\mathsf{A}}_{{flow}}$. We omit the details of the Vandermonde orthogonalization procedure for brevity. Details of the procedure may be found in Brubeck et al. (2021) and an application to potential flow is presented by Baddoo (2020). Arnoldi orthogonalization may be applied to the least-squares problem (4.11) through simple modifications to code given by Brubeck et al. (2021).

4.2. Finding the fluid flow

If all of the bodies in the flow are perfect electrical conductors, then the complex potential describing the fluid flow is obtained by pre-multiplying the electrical potential from § 4.1 by $-\mathrm {i}\sigma B_0$, which follows from (3.2a) and (3.3), since in this case $W_P=0$ as was discussed in § 3.1.

When one or more of the obstacles in the flow is an electrical insulator, the potential $-\mathrm {i}\sigma B_0 W_E$ violates the impermeability condition on electrical insulators. This is illustrated in figure 6(a), where the flow streamlines of $-\mathrm {i}\sigma B_0 W_E$ (shown in grey) clearly penetrate the insulating body. In cases where insulating bodies are present, the fluid flow must thus be obtained in two steps. First, the electrostatic problem must be solved to specify the induced circulation in the fluid flow around each body as in (3.5). Then, the fluid flow boundary-value problem, subject to impermeable boundary conditions on each obstacle, must be solved. When electrical insulators are present in the flow, a non-zero pressure arises whose role is to enforce the impermeability condition.

Figure 6. (a) Numerical solution of the electrostatic problem involving two perfect conductors held at voltages $V=1$ and $V=0$ along with a perfect insulator, obtained though the approach outlined in § 4.1. The charge magnitude on each conducting surface, $Q$, is obtained through the solution of the least-squares problem.(b) The fluid flow induced by the voltage distribution of (a), in the presence of an out-of-plane magnetic field. The flow is computed through the procedure outlined in § 4.2. The magnitude of the circulation around each conductor is specified through the solution presented in (a), $\varGamma =\sigma B_0 Q$. (c) A free-stream fluid flow past the same obstacles in a Hele-Shaw cell in the absence of magnetic effects (no circulation). (d) The free-stream flow from (c) combined with magnetically induced circulation due to the electrical configuration in part (a). The fluid flow depicted in (d) is essentially a superposition of the flows depicted in (b,c).

The complex potential for the fluid flow problem can be written as follows:

(4.12a)\begin{gather} W_{{flow}}(z) = b_0 + \sum_{j\in I_C} \mathrm{i}\sigma B_0 Q_j \log{\left(z-z_i\right)}+\sum_{j=1}^N\sum_{k=1}^{\infty} \frac{b_k^{(j)}}{\left(z-z_j\right)^k} +\bar{U}z, \end{gather}

(4.12b)\begin{gather}\mathrm{Im}\left\{W_E\right\} = \psi_j,\quad z\in \partial B_j,\ j\in\{1,2,\ldots,N\}, \end{gather}

where $\{b_k^{(j)}\}$ and $\{\psi _j\}$ are unknown coefficients, $U$ is a background pressure-driven flow and $Q_j$ are the known charges from the electrostatic solution of § 4.1. Taking $U=0$ gives the situation of a flow driven exclusively by magnetic effects. Note that we can immediately take $\mathrm {Re}\{b_0\}=0$, since the boundary-value problem does not place any restriction on the real part of this constant.

The potential can be approximated by once again truncating the Laurent series centred in each body

(4.13a)\begin{equation} W_{{flow}}(z) \approx b_0 + \sum_{j\in I_C} \mathrm{i}\sigma B_0 Q_j \log{\left(z-z_i\right)}+\sum_{j=1}^N\sum_{k=1}^{N_L^{(j)}} \frac{b_k^{(j)}}{\left(z-z_j\right)^k} +\bar{U}z, \end{equation}

which is still subject to (4.13b).

The sets of coefficients, $b_k^{j}$ and $\psi _k$, can be obtained in a manner similar to § 4.1. In this case, there are no undetermined coefficients associated with logarithm terms. Compared with the electrostatic problem of § 4.1, the matrix $\boldsymbol{\mathsf{A}}$ is changed in two ways in the fluid flow problem. Firstly, logarithm terms appear in $\boldsymbol {f}$ instead of in $\boldsymbol{\mathsf{A}}_{{flow}}$. Second, since all bodies possess an undetermined stream function value, the right partition of $\boldsymbol{\mathsf{A}}$ is changed accordingly. In the flow problem, the flow matrix $\boldsymbol{\mathsf{A}}$ becomes

(4.14) \begin{equation} \boldsymbol{\mathsf{A}}_{{flow}} = \left[ \begin{array}{c|cc|cccc} \boldsymbol{1} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(1)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(1)}\right\} & -\boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{1} & \vdots & \vdots & \boldsymbol{0} & -\boldsymbol{1} & \boldsymbol{0} & \vdots\\ \boldsymbol{1} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N_C)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N_C)}\right\} & \vdots & \boldsymbol{0} & -\boldsymbol{1} & \ddots\\ \boldsymbol{1} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N_C+1)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N_C+1)}\right\} & \vdots & \vdots & \boldsymbol{0} & \ddots\\ \boldsymbol{1} & \cdots & \cdots & \vdots & \vdots & \vdots & \ddots\\ \boldsymbol{1} & \mathrm{Re}\left\{\boldsymbol{\mathsf{V}}^{(N)}\right\} & \mathrm{Im}\left\{\boldsymbol{\mathsf{V}}^{(N)}\right\} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \cdots\\ \end{array} \right], \end{equation}

where $\boldsymbol {1}$ and $\boldsymbol {0}$ represent vectors of ones and zeros of the appropriate length; for example, the vector $\boldsymbol {1}$ in the top row has a length of $N_b^1$ while the $\boldsymbol {1}$ in the bottom row has a length of $N_b^N$. The corresponding coefficient vector becomes

(4.15)\begin{equation} \boldsymbol{c}_{\mathrm{flow}} = \left[ \mathrm{Im}\left\{a_0\right\} \quad\mathrm{Im}\left\{\boldsymbol{a}\right\} \quad \mathrm{Re}\left\{\boldsymbol{a}\right\} \quad \psi_1 \quad \dots \quad \psi_{N}\right]. \end{equation}

The length of the new load vector $\boldsymbol {f}_{{flow}}$ is the total number of sample points on all boundaries, $\sum _{k=1}^N N_b^{k}$ components. For the boundary points on the $k{\mathrm {th}}$ body, the corresponding section of $\boldsymbol {f}_{{flow}}$ is given by a row vector $\boldsymbol {f}^{(k)}$ whose $n{\mathrm {th}}$ component is given by

(4.16)\begin{equation} f^{(k)}_n=\mathrm{Im}\left\{-\sum_{j=1}^{N_C} \mathrm{i}\sigma B_0 Q_j\log(z^{(k)}_n-z_j)-\bar{U}z^{(k)}_n\right\}. \end{equation}

The new load vector is then given by concatenating the vectors corresponding to each boundary so that $\boldsymbol{f}_{\!{flow}}=[\boldsymbol {f}^{(1)}\ \dots \ \boldsymbol {f}^{(N)}]$. The boundary-value problem is then reduced to the following least-squares problem,

(4.17) \begin{equation} \boldsymbol{\mathsf{A}}_{{flow}}\boldsymbol{c}_{{flow}}^T \approx \boldsymbol{f}_{{flow}}^T. \end{equation}

Once again, we note that the solution accuracy can be improved by using Arnoldi orthogonalization on the Vandermonde part of the matrix $\boldsymbol{\mathsf{A}}_{\mathrm {flow}}$ when solving (4.17). Figure 6(b) shows the fluid flow induced by the voltage configuration illustrated in figure 6(a).

By checking the deviation of the imaginary part of the $W_{{flow}}$ from a constant, we can attain an estimate for the error in $\mathrm {Im}\{W_{{flow}}\}$ over the entire domain. Checking the boundary convergence on a set that is 16 times as dense as the sample points, we find an accuracy of ten digits on the circular boundaries when using 40 Laurent terms in each body and 200 uniformly sample points per body. The accuracy on the trefoil shaped boundary in figure 6 is limited to a more modest seven digits owing to the finger-like geometry. In more extreme finger-like geometries, one should consider the placement of poles in addition to the Laurent series used in the present paper, to account for the well-known crowding phenomenon (Gopal & Trefethen 2019a). For work exploiting simple pole placement to achieve accurate solutions in the case of sharp corners and highly curved geometries, see the works of Gopal & Trefethen (2019b), Costa & Trefethen (2023) and Xue, Waters & Trefethen (2024). Note that Baddoo (2020) applied such methods to the potential flow problem past a simply connected body with a corner and he demonstrated rapid convergence when simple poles are exponentially clustered near corners.

5. Discussion and conclusion

We have analysed the flow of an uncharged ohmic fluid inside a magnetically driven Hele-Shaw cell, at low Hartmann numbers. The problem was first cast the into a complex variables framework. Within this framework, we elucidated – both physically and mathematically – the mechanism by which an external magnetic field may convert electric current into a tuneable fluid flow circulation. Whereas pressure-driven Hele-Shaw flows exhibit identically zero circulation around any closed contour, the circulation in magnetically driven Hele-Shaw flows need not vanish. Moreover, we have demonstrated that by changing the voltage difference applied between conducting probes in the fluid, as well as the probe geometries, one can control the induced fluid flow circulation.

Within our complex variables framework, we presented mathematical solutions, in terms of the prime function, for a class of geometries involving circular conductors; one such solution describes an experiment of David et al. (2023) whose geometry is illustrated in figure 3. Notably, the prime function has a series representation in doubly connected geometries (Baddoo et al. 2020). We note that more exact solutions in periodic domains are also available through the theoretical developments of Baddoo & Ayton (2021); however, such results were not explored in the present paper.

We subsequently noted two limitations of the aforementioned mathematical solutions. First, they are only applicable in canonical circular domains. To obtain a solution in non-circular geometries, a conformal map between a canonical domain and the physical domain of interest must be known. As such, the class of solutions that may be written explicitly in terms of the prime function is limited. Second, the exact solutions do not apply when any of the obstacles in the flow is electrically insulating. Moreover, it should be noted that the prime function must also be computed numerically in domains with a connectivity greater than two.

Because of these limitations, we proceeded to present a numerical scheme capable of finding the fluid flow in more general geometries and in situations where obstacles in the flow may be either insulating or conducting. We outlined a procedure for obtaining accurate Laurent series solutions for the flow in a manner similar to that described by Trefethen (2018).

Approximate solutions must be obtained in two steps. First, one must obtain a series solution to the relevant electrostatic problem. The electrostatic field $\boldsymbol {E}$ is then converted to an electric current via Ohm's law, $\boldsymbol {J}=\sigma \boldsymbol {E}$. The induced flow circulation is related to the product of the net current $I_j$ leaving each conducting body and the external magnetic field strength, $B_0$. Both geometry and applied voltage magnitudes thus affect the magnitude of induced circulations. The key result of the first step is that the electrostatic problem specifies the fluid flow circulation around each conducting body. With the circulation around each body fixed, the potential flow problem is fully posed. The second step of the solution is to find a series solution for the potential flow subject to impermeable boundary conditions on each obstacle and with the circulations around each obstacle imposed by the electrostatic problem. A uniform stream, or other background flows, can be added at this stage by superposition.

Note that in the magnetohydrodynamic Hele-Shaw cell, the circulation around each body in the flow is specified. In the absence of magnetic fields, the velocity potential in a Hele-Shaw cell is proportional to the pressure and is therefore single valued, corresponding to a flow with identically zero circulation. Meanwhile, when a magnetic field is introduced, the Lorentz force generates a well-defined and calculable amount of circulation around each body. This situation is in contrast with potential aerofoil theory, where the circulation around each body is unknown. Since the velocity potential in aerofoil theory is not related to a physical quantity (such as pressure), no physical constraint enforces the circulation value around each body. In some special cases, such as for aerofoils with a single sharp corner, Kutta's condition specifies a unique circulation to de-singularize the velocity field. Recent developments by Gonzalez & Taha (2022) conjecture a new criterion to replace the Kutta condition in more general aerofoil geometries.

Magnetohydrodynamic flow control, as explored in the present study, has many potentially interesting applications in, for example, microfluidic devices. The voltage difference applied between probes, in the magnetically driven flow, may be actively controlled during an experiment to divert flow along desired paths. Probe geometries may also be designed to achieve desired flows. More engineering applications that might be explored in future studies are outlined by Bau (2022). Future theoretical work might examine the possibility of bubble manipulation using magnetically driven flows in a Hele-Shaw cell (Booth, Griffiths & Howell 2023). In realizing such applications, the principles and solution methods of the present paper may be useful for quick and iterative design. Probe geometries and voltages can be adjusted very simply within the framework of § 4 making it and attractive tool for iterative design.