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 Reference Gonzalez and Taha2022). 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 Reference Ahlfors2010, p. 65). A significant literature exists regarding bounds and comparison theorems for the capacitance of different geometries (Ahlfors Reference Ahlfors2010, 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 (Reference Crowdy2020).

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 (Reference Crowdy2020, 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. Reference Baddoo, Kurt, Ayton and Moored2020). 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 (Reference Crowdy2020, p. 294) and Crowdy Reference Crowdy2006b,Reference Crowdya).

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 Reference Crowdy2020).

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 (Reference Crowdy2020). 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 Reference Crowdy2020, 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 Reference Crowdy2020, 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 (Reference Crowdy and Marshall2007) and Crowdy (Reference Crowdy2020, chap. 14).

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 (Reference Crowdy2020). 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 Reference Crowdy2007).

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 (Reference Crowdy2020) 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 (Reference Crowdy2020) 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. (Reference David, Hester, Xu and Aurnou2023) 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 (Reference Crowdy2020) 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 (Reference Crowdy and Marshall2007) and Crowdy (Reference Crowdy2020, chap. 14) present efficient methods for computing the prime function. Codes for computing the prime function are readily available (Crowdy, Kropf & Green Reference Crowdy, Kropf, Green and Nasser2016).

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. (Reference David, Hester, Xu and Aurnou2023), 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 (Reference Miyoshi and Crowdy2023, § 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 Reference Vasconcelos, Marshall and Crowdy2015) 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.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. (Reference Brubeck, Nakatsukasa and Trefethen2021) and an application to potential flow is presented by Baddoo (Reference Baddoo2020). Arnoldi orthogonalization may be applied to the least-squares problem (4.11) through simple modifications to code given by Brubeck et al. (Reference Brubeck, Nakatsukasa and Trefethen2021).

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 Reference Gopal and Trefethen2019a). 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 (Reference Gopal and Trefethen2019b), Costa & Trefethen (Reference Costa and Trefethen2023) and Xue, Waters & Trefethen (Reference Xue, Waters and Trefethen2024). Note that Baddoo (Reference Baddoo2020) 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.