2 Experimental set-up and procedure

The experimental set-up is shown in figure 2. It consisted of an open cylindrical cavity with a diameter of 85.8 mm and a depth of 25.4 mm carved in an acrylic block. A thin copper ring of 10 mm height which served as the outer electrode was fitted into the container while a copper cylinder of 25.4 mm in diameter placed at the centre of the container acted as the inner electrode. The width of the gap between the electrodes was 30.2 mm. The cavity was filled with a weak electrolytic solution of sodium bicarbonate (NaHCO₃; 8.6 % by weight). Experimentally, the thickness of the undisturbed fluid layer was controlled by pouring a prescribed volume of electrolyte into the cavity. The container was placed on top of a vertically polarised permanent disk magnet with a diameter of 80 mm. An electric potential difference was applied between the electrodes, and, therefore, a radial electric direct current of magnitude $I_{0}$ traversed the layer of conducting fluid. The interaction of the radial current and the predominantly vertical magnetic field led to the appearance of an azimuthal Lorentz force which caused fluid motion.

Figure 2. Experimental set-up: 1, disk magnet; 2, acrylic glass block; 3, outer copper electrode; 4, inner copper electrode; 5, electric current source; 6, wires; 7, non-magnetic insert; 8, electrolyte layer.

When viewed from above, the fluid rotated clockwise (anticlockwise) when the applied current flowed radially outwards (inwards). The physical properties of the electrolyte and electrodes, the geometrical parameters of the experimental device as well as the characteristic values of the applied magnetic field and electric current are given in table 1. Experiments were performed for various values of the applied current and the electrolyte layer thickness.

The experimental flow was analysed with two different methods.

First, the flow patterns were visualised using a dye tracer. The experiments started from rest with a small blob of liquid dye dripped from the top and initially remaining near the free surface of the quiescent fluid. To demonstrate the three-dimensionality of the flow, the dye was also dripped from further above the layer surface so that it would penetrate deeper on impact. Once the applied current was fixed, the development of the flow pattern produced by the Lorentz force was revealed by smearing of the dye blob and, eventually, the appearance of vortices if the applied current was sufficiently strong. The flows were recorded using a Sony Handycam videorecorder (model HDR-XR100, with a resolution of $1080\times 720$ pixels and a capture frequency of 30 f.p.s.) placed above the experimental set-up. The recorded flow patterns were used to identify the onset of instability and the number of vortices that appeared for given conditions.

Second, detailed measurements of the free surface velocity field were carried out using particle image velocimetry (PIV). The experiments with PIV also started from rest, and a uniform illumination was provided by a set of three LED lamps placed on top of the set-up around the liquid container. Due to the relatively low velocities (of the order of centimetres per second) of the analysed flows, the same videorecorder was also used for the PIV analysis. Silver-covered glass spheres with an average diameter of 10  $\unicode[STIX]{x03BC}$ m were used as tracers. A centrifuge was used to select only the neutrally buoyant spheres for PIV measurements. The subsequent analysis was performed using the PIVlab toolbox in MATLAB. A more detailed discussion of the experimental set-up and procedure can be found in Pérez-Barrera et al. (Reference Pérez-Barrera, Pérez-Espinoza, Ortiz, Ramos and Cuevas2015).

3 Mathematical problem formulation and governing equations

We consider an electrolyte layer of thickness $h$ over a non-conducting solid bottom with a free top surface. The layer is contained between two vertical coaxial cylindrical electrodes with radii $R_{1}$ and $R_{2}$ , as shown in figure 3. The system is placed above a permanent disk magnet which creates a predominantly vertical magnetic field. Since this field is generally non-uniform, we characterise its strength by the field magnitude $B_{0}$ in the corner between the solid bottom and the inner electrode.

Figure 3. Sketch of the problem geometry and non-dimensional coordinate system.

When an electric potential difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}$ is applied between the electrodes, the total current

flows predominantly radially through the layer. Expression (3.1) is obtained upon adopting two assumptions, namely that the surfaces of the electrodes are essentially equipotential and that the influence of the electrolyte flow on the distribution of electric potential in a layer is negligible. The former is justified by the fact that the conductivity of the electrodes (copper) is seven orders of magnitude higher than that of the electrolyte (see table 1), while the validity of the latter is due to the negligibility of magnetic induction effects which will be demonstrated numerically in the subsequent sections.

The interaction of the current and the applied magnetic field creates a Lorentz force $\boldsymbol{F}_{L}=\boldsymbol{j}\times \boldsymbol{B}$ , where $\boldsymbol{j}$ and $\boldsymbol{B}$ are the applied current density and the magnetic induction field respectively. This predominantly unidirectional azimuthal force drives the flow which structure we investigate here.

Under the small-magnetic-Reynolds-number approximation (Davidson Reference Davidson2001; Pérez-Barrera, Ortiz & Cuevas Reference Pérez-Barrera, Ortiz, Cuevas, Klapp, Sigalotti, Medina Ovando, López Villa and Ruíz Chavarría2016), the non-dimensional equations governing the steady axisymmetric fluid motion are the Poisson equation for the electric potential, the momentum equations and the continuity equation for an incompressible fluid. They are written in cylindrical coordinates as

where $\unicode[STIX]{x1D719}$ is the electric potential, $(u_{r},u_{\unicode[STIX]{x1D703}},u_{z})$ and $(B_{r},0,B_{z})$ are the velocity and magnetic field components in the radial ( $r$ ), azimuthal ( $\unicode[STIX]{x1D703}$ ) and vertical ( $z$ ) directions respectively and $p$ is the pressure, including the hydrostatic component. From Ohm’s law, the electric current density components are

Here, the aspect ratio of the electrolyte layer $\unicode[STIX]{x1D716}$ , the square of the Hartmann number $Ha^{2}$ characterising electromagnetic effects, the square of the Froude number $Fr^{2}$ describing the influence of gravity on a free surface flow and the Reynolds number $Re\,$ quantifying viscous effects are

where the velocity scale is defined as $U_{0}=\unicode[STIX]{x1D70E}_{e}\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{0}B_{0}h^{2}/2\unicode[STIX]{x1D707}(R_{2}-R_{1})$ . The typical experimental values of these parameters are $\unicode[STIX]{x1D716}\sim 10^{-1}$ , $Ha^{2}\sim 10^{-5}$ , $Fr^{2}\sim 10^{-1}$ and $Re\,\sim 10^{3}$ . Equations (3.2)–(3.9) are made non-dimensional using the following scalings:

where the star symbol denotes dimensional quantities.

The no-slip/no-penetration velocity boundary conditions at the cylindrical electrodes and non-conducting bottom are

where $\unicode[STIX]{x1D6FC}=(R_{2}+R_{1})/(R_{2}-R_{1})$ ( $\unicode[STIX]{x1D6FC}\approx 1.84$ in our experiments and computations). As will be justified in § 8, we neglect the deformation of the free surface (see also Satijn et al. Reference Satijn, Cense, Verzicco, Clercx and van Heijst2001) and apply continuity of the tangential stress condition. This leads to

The boundary conditions for the electric potential at the electrode surfaces are

The condition $j_{z}=0$ at the non-conducting bottom and at the free surface leads to

4 Magnetic field

The reconstruction of the realistic applied magnetic field is of fundamental importance for proper modelling of the flow. To calculate a non-uniform magnetic field created by a permanent vertically polarised magnet, we assume that it is homogeneously magnetised (i.e. its bulk magnetisation $M_{z}$ is assumed to be constant), so that the characteristics of the field depend only on the shape of the magnet and the distance from its surface. Experiments conducted with rectangular magnets revealed that the observed flow patterns remain largely the same as those arising when a smaller disk magnet is used as long as the horizontal dimensions of the rectangular magnet are larger than the diameter of the disk magnet, and the diameter of the magnet is not smaller than the outer diameter of the experimental container, so that no magnetic field reversal occurs within the electrolyte layer. Therefore, for definiteness, all computations reported here are performed for a vertically polarised disk magnet of height $2c$ with the radius equal to the radius $R_{2}$ of the outer electrode. Due to axial symmetry, the azimuthal component of the magnetic field of such a magnet is $B_{\unicode[STIX]{x1D703}}^{\ast }=0$ , as has already been taken into account in the governing equations. Experimentally, the magnitude of the magnetic field driving the flow is controlled by changing the width $d$ of the gap between the upper surface of the magnet and the bottom of the electrolyte layer ( $d$ is the sum of the bottom container wall thickness and the air gap between the container and the magnet); see figure 3 and table 1 for details. In particular, variation of the gap width $d$ from 6 to 27 mm leads to a decrease of the maximum magnetic field intensity by a factor of approximately 3 (from approximately 0.06 to 0.02 T in experiments). The distribution of the magnetic field also changes, as demonstrated in figure 4. In particular, comparison of $(a)$ and $(b)$ indicates that the uniform field assumption commonly used in similar studies is quickly invalidated away from the magnet surface, so that, in general, the gap width $d$ needs to be considered as an independent experimental control parameter, the influence of which cannot be accounted for in full by a simple rescaling of the characteristic magnetic field strength $B_{0}$ .

Figure 4. Scaled vertical (solid lines) and radial (dashed lines) components of the magnetic field along the bottom (thick lines) and the free surface (thin lines) of the electrolyte layer of depth $h=7.5$ mm placed $(a)$ $d=6$ mm and $(b)$ $d=27$ mm above the disk magnet. Only the magnetic field in the flow region is shown. The outer edge of the disk magnet coincides with the right edges of the graphs.

The magnetic field induced by the weak electric current flowing between the electrodes is several orders of magnitude smaller than the magnetic field $\boldsymbol{B}^{\ast }$ of the permanent magnet. Thus, once the induced magnetic field is neglected, $\boldsymbol{B}^{\ast }$ has to satisfy the following two conditions:

To compute the magnetic field, we choose to introduce a scalar magnetic potential $\unicode[STIX]{x1D719}_{m}^{\ast }$ , so that $\boldsymbol{B}^{\ast }=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{m}$ , the second of conditions (4.1a,b ) is satisfied identically and the first one results in the Laplace equation $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{m}^{\ast }=0$ . Subsequently, $\unicode[STIX]{x1D719}_{m}^{\ast }$ is computed (up to a multiplicative constant) at discrete collocation points required by the numerical method (described in Suslov & Cuevas (Reference Suslov and Cuevas2017)) by evaluating numerically the double integrals in

using the standard MATLAB two-dimensional quadrature routine quad2d. However, due to the finite accuracy of the numerical integration, the thus-computed scalar magnetic potential satisfies the Laplace equation only approximately. Therefore, in further computations, we use the corrected magnetic potential $\unicode[STIX]{x1D6F7}_{m}^{\ast }=\unicode[STIX]{x1D719}_{m}^{\ast }-\unicode[STIX]{x1D719}_{m}^{\ast \prime }$ , where the correction $\unicode[STIX]{x1D719}_{m}^{\ast \prime }$ is obtained by solving the Poisson equation $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{m}^{\ast \prime }=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}_{m}^{\ast }$ , with the homogeneous condition $\unicode[STIX]{x1D719}_{m}^{\ast \prime }=0$ along the boundary of the computational domain. The numerical values of $\unicode[STIX]{x1D6F7}_{m}^{\ast }$ are scaled so that $\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}_{m}^{\ast }/\unicode[STIX]{x2202}z\ast =-B_{0}$ at the reference location shown by a black circle in figure 3. The typical distributions of the vertical and radial components of the magnetic field along the bottom of the electrolyte layer and the free surface are shown in figure 4.

5 Quasi-two-dimensional and thin-layer approximations

The symmetry of the problem suggests that an approximate steady two-dimensional azimuthal flow solution $\boldsymbol{u}=(0,u_{\unicode[STIX]{x1D703}}(r,z),0)$ , $\boldsymbol{j}=(j_{r}(r,z),0,j_{z}(r,z))$ , $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}(r,z)$ may exist which satisfies the following reduced equations:

It should be noted that while the $Ha\rightarrow 0$ limit is appropriate in the context of our experimental studies of electrolyte flows, we choose not to drop terms multiplied by $Ha$ from the above equations. As will be seen from the subsequent discussion, this does not overly complicate the algebraic manipulations but will enable us to derive a more general analytical solution which remains valid for liquid metal flows characterised by small but finite Hartmann numbers. Given that $u_{\unicode[STIX]{x1D703}}=u_{\unicode[STIX]{x1D703}}(r,z)$ , as is enforced by the no-slip boundary condition at the bottom of the layer, (5.3) and (5.4) contradict one another. This means that the problem solution is inherently three-dimensional, with all velocity components being non-zero. The quasi-two-dimensional solution with $p\approx \text{const.}$ can only exist in the limit

that is, in thin layers where $\unicode[STIX]{x1D716}\rightarrow 0$ . Upon assuming that $B_{z}\approx 1$ and $B_{r}\approx 0$ (these are reasonable assumptions if $d\rightarrow 0$ , i.e. if the bottom of the fluid layer is located close to the surface of a magnet with large lateral dimensions), (5.1) and (5.2) become

Away from the cylindrical electrodes, these equations have the exact solution

Figure 5. (a,b) The quasi-two-dimensional azimuthal velocity profiles and (c) the Rayleigh criterion for inviscid centrifugal instability for an electrolyte layer placed in a uniform vertical magnetic field at $\unicode[STIX]{x1D716}=0.083$ , $Ha=1.69\times 10^{-3}$ , $Re\,=1.1\times 10^{2}$ and $Fr^{2}=4.4\times 10^{-3}$ (corresponding to $B_{0}=B_{z}=0.02$  T,  $I_{0}=0.05$  A and $h=2.5$  mm). The profiles are shown for $z=1.0$ , $0.5$ and $-0.5$ in panels (a) (top to bottom) and (c), and for $r=1.30$ , $1.98$ and $2.60$ (right to left) in panel (b). The solid and dashed lines represent the numerical and exact solutions (5.9) respectively. The dash-dotted and dotted lines show the depth-averaged profiles obtained from (5.7) and (5.19) respectively.

It is worthwhile to point out that the distribution of electric potential given by (5.8) is independent of $Ha$ , even though it enters the expression (5.9) for the azimuthal velocity and the governing equation (5.6). However, this solution does not satisfy the no-slip boundary conditions at the surface of the electrodes, where boundary layers are formed. Thus, the problem is naturally rendered to a singular perturbation analysis similar to that reported, for example, in Delacroix & Davoust (Reference Delacroix and Davoust2014). In the current paper, however, we do not pursue such an analysis any further and resort to the high-accuracy numerical spectral solution shown in figure 5. As seen from this figure, away from the cylindrical walls, the analytical and numerical solutions coincide within the plotting resolution, but deviate from one another significantly in the boundary layers (consider also approximate solutions reported in Digilov (Reference Digilov2007) for an infinitely deep layer).

It is important to note that in thin electrolyte layers, the maximum of the azimuthal velocity is necessarily located near the inner cylinder. As a result, the largest flow shear is found at small values of $r$ . This is where the shear instability would be expected to set first, which, however, is not seen in our experiments. The variation of the geometry and strength of the magnetic field created by a realistic magnet and shown in figure 4 is expected to only enhance this effect, as it weakens the driving Lorentz force near the outer edge of the layer.

In figure 5(c), the rate of change of the angular momentum of the fluid is shown. As can be derived from the Rayleigh criterion for inviscid centrifugal instability (Rayleigh Reference Rayleigh1916), such an instability can only occur if $\text{d}/\text{d}r(ru_{\unicode[STIX]{x1D703}})<0$ . Since $u_{\unicode[STIX]{x1D703}}\sim 1/r$ , see (5.9), the bulk of the flow is marginally stable with respect to this mechanism. However, the boundary layer near the outer electrode satisfies this necessary condition and thus centrifugal instability should be considered as a potential physical mechanism that can lead to the appearance of the vortices.

To find the approximate radial profile of the azimuthal velocity component in the complete flow domain, Pérez-Barrera et al. (Reference Pérez-Barrera, Ortiz, Cuevas, Klapp, Sigalotti, Medina Ovando, López Villa and Ruíz Chavarría2016) suggested that one should proceed with a depth-averaged model that is commonly used in studies of thin fluid layers (see, for example, Sommeria (Reference Sommeria1988), Andreev, Heberstroh & Thess (Reference Andreev, Heberstroh and Thess2001) for discussion in an MHD context). Such a model is obtained by integrating (5.6) and (5.7) over the layer thickness and taking into account the boundary conditions (3.12), (3.13) and (3.15). This results in

where the overbar denotes depth averaging,

Integrating (5.10) once and eliminating $\bar{\unicode[STIX]{x1D719}}$ from (5.11) results in

Since the solution (5.8) for the electric potential is regular in $r$ , it can be used to eliminate the integration constant $C$ from the above equations,

It is tempting to use (5.9) to obtain an explicit expression for the bottom friction. However, because of the boundary layer singularity, this cannot be justified near the electrodes, where the explicit expression for the bottom friction is not known. Thus, instead, we adopt the quasi-two-dimensionality assumption (see, for example, the discussion in Dolzhanskii, Krymov & Manin (Reference Dolzhanskii, Krymov and Manin1992)), and parameterise it as

where the value of the proportionality constant can be determined using (5.9),

Expression (5.17) can be written in dimensional form as

In the context of a quasi-two-dimensional approximation (Dolzhanskii et al. Reference Dolzhanskii, Krymov and Manin1992), the first term in (5.18) is recognised as Rayleigh friction. It accounts for the viscous friction at the bottom of the layer. The second term represents Hartmann friction which characterises the flow resistance due to the electromagnetic effect. It only becomes essential in fluids with high conductivity, such as liquid metals, and plays a negligible role in the electrolytes considered here. Finally, neglecting the higher-order terms in $Ha$ , we obtain the forced Bessel equation for the depth-averaged azimuthal velocity profile,

with homogeneous boundary conditions

This equation is structurally similar to that discussed in Pérez-Barrera et al. (Reference Pérez-Barrera, Ortiz, Cuevas, Klapp, Sigalotti, Medina Ovando, López Villa and Ruíz Chavarría2016) but contains rigorously derived coefficients. Its solution, shown by the dotted lines, is compared with the depth-averaged solution of (5.7), shown by the dash-dotted lines, in figure 5. The two solutions are virtually indistinguishable except for a small difference in the boundary layers, where the one-dimensional solution slightly overestimates the depth-averaged quasi-two-dimensional solution. This is traced back to the inaccuracy of assumption (5.16) illustrated in figure 6, where the approximation error of (5.16),

is presented. In the vicinity of cylindrical walls, approximation (5.16) underestimates the bottom friction and thus predicts a slightly faster depth-averaged velocity than the full set of equations. It should be noted, however, that even though the overall performance of the model (5.19) is quite satisfactory, the depth-averaged velocity and thus the average flow shear are significantly smaller than the maximum local values that define the stability properties of the flow. Therefore, it is unlikely that the depth-averaged model could be used to produce accurate information on the stability of the considered flows. Thus, here, we will not conduct any comprehensive computations using the depth-averaged equations, and resort to a numerical solution of the quasi-two-dimensional and axisymmetric three-dimensional equations.

Figure 6. The difference between the bottom friction values calculated from (5.7) and using the approximation (5.16).

6 Numerical quasi-two-dimensional solutions for realistic magnetic fields

When a realistic magnetic field is considered, as shown in figure 4(a,b), no analytical solutions are readily available. Thus, (5.1) and (5.2) are solved numerically using the Chebyshev pseudo-spectral collocation method proposed by Ku & Hatziavramidis (Reference Ku and Hatziavramidis1984) and Hatziavramidis & Ku (Reference Hatziavramidis and Ku1985) and adapted for large-scale flow stability calculations by Suslov & Paolucci (Reference Suslov and Paolucci1995a ,Reference Suslov and Paolucci b ). The generalisation of the method for the multi-dimensional spatial approximation used here and the details of implementation are discussed in Suslov & Cuevas (Reference Suslov and Cuevas2017). The numerical results reported here were obtained using $45\times 31$ spectral modes, which ensures convergence to at least $10^{-6}$ . The numerical solutions were also validated in the limiting cases against the analytical results; see figure 5.

Figure 7. The difference between the numerical solution $\unicode[STIX]{x1D719}$ of (5.1) with a realistic magnetic field and the approximate analytical solution for the electric potential $\unicode[STIX]{x1D719}_{0}$ given by (5.8) for $z=1$ , $0.5$ , $-0.5$ and $-1$ (bottom to top) and the parameter values specified for figure 5.

When the magnetic field is given and fixed, the governing equations are linear, and their solutions are obtained by a direct matrix inversion. A comparison of the thus-obtained fully numerical and approximate analytical solutions for the electric potential is given in figure 7. For the small Hartmann numbers that are of interest here, the difference between the two is negligible. Therefore, with a very good accuracy, the expression (5.8) can be used, and (3.2) and (5.1) can be eliminated from the respective systems of governing equations.

Figure 8. The same as figure 5 but for a realistic magnetic field created by a disk magnet (see figure 4 a) placed 6 mm below the bottom of the electrolyte layer.

In figure 8, the azimuthal velocity profiles obtained for a realistic magnetic field (shown in figure 4) with the same reference magnitude $B_{0}$ as that of the uniform vertical magnetic field ( $B_{z}=1$ , $B_{r}=0$ ) used in figure 5 are presented. Comparison of the two figures demonstrates a remarkable robustness of the velocity profiles. Even though the details of the magnetic fields used to obtain the two solutions are drastically different, the velocity profiles remain qualitatively unchanged. However, quantitatively, the velocity profiles become fuller in the middle part of the flow region away from the electrodes: the profiles lie somewhat higher (to the right) in figure 8(a) (figure 8 b) than the dashed lines representing the reference analytical solutions. This trend continues as the depth of the layer and the electric current increase, and the shape of the velocity profile eventually deviates from that given by the depth-averaged model; see figure 9. The curvature of the velocity profiles becomes uniformly negative, but the largest velocity gradient is still found near the inner cylinder. This suggests that a potential shear instability would still set in preferentially near the inner cylinder. However, this is not what is observed experimentally.

As is seen from figure 8(c), in a realistic magnetic field, the region near the outer wall, where the Rayleigh criterion for inviscid centrifugal instability is satisfied, widens at the expense of the central core flow. However, the magnitude of the radial derivative of the angular momentum is reduced compared with the idealised uniform magnetic field situation of figure 5. Therefore, it also remains unclear whether the inviscid centrifugal mechanism of instability can be directly responsible for the formation of the observed vortices in thin layers of viscous electrolytes. In turn, this raises doubts about whether the quasi-two-dimensional approximation given by (5.1) and (5.2) can be sufficiently accurate to describe the flow at hand in principle, even though the value of the Hartmann number characterising the Lorentz force driving the flow remains small. Indeed, it should be recalled that the validity of the quasi-two-dimensional approximation is dictated by the value of the square of the Froude number defined in (5.5), which increases linearly with the depth of the layer and quadratically with the magnitude of the current flowing between the electrodes. As discussed in § 5, when $Fr^{2}$ is not negligible, an essentially three-dimensional flow with non-zero radial and vertical velocity components $u_{r}$ and $u_{z}$ develops, which we investigate next.

Figure 9. Numerical solutions of (5.2) for the quasi-two-dimensional azimuthal profiles for an electrolyte layer placed in a magnetic field created by a disk magnet at $\unicode[STIX]{x1D716}=0.248$ , $Ha=5.08\times 10^{-3}$ , $Re\,=1.34\times 10^{3}$ and $Fr^{2}=0.213$ (corresponding to $B_{0}=0.02$  T,  $I_{0}=0.2$  A,  $h=7.5$  mm and $d=6$  mm). The profiles are shown for $z=1.0$ , $0.5$ and $-0.5$ in panels (a) (top to bottom) and (c), and for $r=1.30$ , $1.98$ and $2.60$ (right to left) in panel (b).

7 Steady toroidal flows

Figure 10. (a) Radial profiles of the azimuthal (solid lines) and vertical (dashed lines) velocity components $u_{\unicode[STIX]{x1D703}}$ and $u_{z}$ , (b) vertical $u_{\unicode[STIX]{x1D703}}$ profiles and (c) the Rayleigh criterion for inviscid centrifugal instability for the type 1 steady three-dimensional solution of (3.2)–(3.9) for the same values of the parameters as in figure 9. The dotted lines in (a) show the Stewartson boundary layer solutions (Stewartson Reference Stewartson1957).

Figure 11. The experimental (symbols) and computed (dashed line) azimuthal velocity profiles at the free surface for $h=5$ mm and $I_{0}=100$ mA. The circles show the averaged PIV data. The solid line extrapolates the experimental velocity profile towards the walls, where accurate PIV measurements are impossible, which explains the discrepancy between the computed and measured profiles in the near-wall regions.

Figure 12. (a) Horizontal velocity hodograph (Ekman-type spiral) and (b) the corresponding Ekman angle for the type 1 steady three-dimensional solution of (3.2)–(3.9) for the same values of the parameters as in figure 9. The dashed line in (b) shows the values for a bottom Ekman layer of infinite lateral extent.

In this section, we discuss steady axisymmetric solutions of the system of governing equations (3.2)–(3.9). Given that the three-dimensional governing equations are nonlinear, we solve them using Newton-type iterations. An example of the velocity profiles obtained by solving a three-dimensional axisymmetric problem is shown in figure 10. Despite the fact that the fields were computed for the same values of the governing parameters as used for figure 9, the azimuthal velocity profiles are strikingly different. First, the maximum value of the azimuthal velocity component of the three-dimensional axisymmetric solution is more than three times smaller than that obtained from a quasi-two-dimensional approximation. Second, the location of the maximum velocity shifts from the inner to the outer cylinder. This is in good agreement with the experimentally measured surface velocity profile shown in figure 11. Third, the shape of the three-dimensional axisymmetric velocity profile becomes a sensitive function of the distance from the bottom of the layer: near the free surface, the fastest flow occurs close to the outer cylinder (see lines 1 and 2 in figure 10 a), while closer to the bottom, the azimuthal velocity profile becomes more of a ‘plug-flow’ type (line 3 in figure 10 $a$ ). Fourth, figure 10 $(b)$ demonstrates that, somewhat counterintuitively, the fastest azimuthal flow occurs in the bulk of the layer rather than at the free surface, as predicted by the quasi-two-dimensional approximation. Fifth, figure 10(c) shows that one of the major effects of the flow three-dimensionality is a significant redistribution of the fluid angular momentum in such a way that the core of the flow not only retains its marginal stability with respect to centrifugal instability, but actually becomes stable, so that $\text{d}/\text{d}r(ru_{\unicode[STIX]{x1D703}})\approx \text{const.}>0$ there. This is mostly due to the development of the radial depth-dependent cross-flow depicted in figure 12(a), where the hodograph of the horizontal velocity is shown. It has the characteristic shape of an Ekman spiral (Ekman Reference Ekman1905), with a positive (outward) radial velocity component $u_{r}$ near the free surface and negative (inward) $u_{r}$ near the bottom of the layer. The direction of the horizontal flow is shown in figure 12(b) by the angle $\unicode[STIX]{x1D6FE}=\tan ^{-1}(u_{r}/u_{\unicode[STIX]{x1D703}})$ as a function of the vertical coordinate $z$ . Similarly to the bottom Ekman layer near the solid surface and close to the inner cylindrical boundary, this angle tends to $-45^{\circ }$ . However, comparison with the analytical results for an unbounded Ekman layer, shown by the dashed line in figure 12(b), indicates that this similarity is rather superficial. It breaks both near the free surface and in the vicinity of the outer cylindrical wall.

It might appear that the considered flow configuration is closer to that of the Bödewadt boundary layer existing in an infinite layer of fluid rotating over a stationary solid plane (Bödewadt Reference Bödewadt1940). However, computational results presented by MacKerrell (Reference MacKerrell2005) indicate that the velocity vector in such a layer would form an angle of $\unicode[STIX]{x1D6FE}\approx -50.6^{\circ }$ with the azimuthal direction near the solid bottom. Davidson & Pothérat (Reference Davidson and Pothérat2002) considered the Bödewadt–Hartmann layer that occurs when conducting fluid placed in a constant vertical magnetic field rotates over a stationary solid bottom. The approximate solutions obtained by these authors were given in terms of the Elsasser number $A$ , which in our notation can be written as

In this limit, the solution provided by Davidson & Pothérat (Reference Davidson and Pothérat2002) predicts $\unicode[STIX]{x1D6FE}\approx -59^{\circ }$ at the bottom of the layer. Although both of these set-ups appear to be similar to the one considered here, our computations do not confirm the existence of such large angles in our flows. Therefore, while offering certain guidance, the numerous stability results available in the literature for Ekman and Bödewadt layers (see, for example, Faller (Reference Faller1963, Reference Faller1991), Lilly (Reference Lilly1966), Savaş (Reference Savaş1987), MacKerrell (Reference MacKerrell2005) and references therein, to name a few) are not likely to be fully relevant to stability studies of the flows of interest here.

It is also instructive to compare the obtained radial $u_{\unicode[STIX]{x1D703}}$ profiles for the axisymmetric three-dimensional solutions shown in figure 10(a) with those existing on a stationary solid cylindrical surface surrounding a fluid rotating almost as a solid body. In our notation, the leading term of the boundary layer solution derived by Stewartson (Reference Stewartson1957) reads as

where $U$ is the effective azimuthal velocity at the edge of the boundary layer. In the bulk of the fluid, this solution predicts a shear-free flow, and thus it is very different from the velocity profiles computed for our problem. However, agreement could reasonably be expected near the outer cylindrical wall. Yet this does not appear to be the case either: figure 10(a) shows that Stewartson boundary layer is somewhat thicker than that existing in the problem considered here.

Figure 13. Meridional velocity (a) and azimuthal vorticity (b) fields for the type 1 three-dimensional solution of (3.2)–(3.9) for the same values of the parameters as in figure 9. The light (dark) regions in (a) correspond to large (near-zero) values of the azimuthal velocity component $u_{\unicode[STIX]{x1D703}}$ . The light (dark) regions in (b) correspond to positive (negative) values of the azimuthal vorticity $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}$ . (c) Streamlines projected onto the meridional plane.

The dashed lines in figure 10(a) show that relatively weak vertical flows (it should be recalled that their magnitude is characterised by $\unicode[STIX]{x1D716}u_{z}$ , while $u_{z}$ is shown in the figure) exist near the cylindrical walls. However, they have a rather complicated non-monotonic structure. To make the visualisation of a complete flow field easy, this is shown in figure 13, where the non-dimensional transverse flow velocity $(u_{r},\unicode[STIX]{x1D716}u_{z})$ , the azimuthal component $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}$ of the vorticity vector

and the azimuthal projection of the streamlines onto the meridional plane are shown (it should be recalled that $u_{r}=u_{z}=\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}\equiv 0$ if the quasi-two-dimensional approximation is used). Even though both the maximum radial and maximum vertical velocity components remain substantially smaller than the maximum azimuthal velocity (in the example shown, by factors of 3.2 and 5.5 respectively) they have a strongly pronounced effect on the overall flow and, ultimately, on its stability. Since, in the considered geometry, fluid moves circumferentially, its inertia attempts to drive it towards the outer cylinder in the meridional plane shown in figure 13(a). However, the top–bottom symmetry is broken by the different dynamic conditions (zero tangential stress at the top boundary and no-slip at the bottom boundary). As a result, a flow occurs in the radial direction, with fluid moving outward along the free surface, returning to the centre near the bottom and rising (sinking) near the inner (outer) cylindrical boundary. The corresponding vortex is illustrated in figure 13(b), with meridional streamlines shown in (c).

It is important to note here that, in contrast to the two-dimensional approximation, the three-dimensional azimuthal velocity profile possesses the largest shear near the outer cylinder, as can be seen by comparing figures 9(a) and 10(a). This indicates that the likelihood of a hydrodynamic instability is expected to be higher in the outer flow region, which is indeed observed in experiments. Thus, it is reasonable to conclude that the experimentally observed vortical structures appearing near the outer cylinder are caused by the three-dimensionality of the base axisymmetric flow.

Figure 14. The same as figure 10 but for the type 2 three-dimensional solution of (3.2)–(3.9) for the same values of the parameters as in figure 9.

Figure 15. The same as figure 12 but for the type 2 three-dimensional solution of (3.2)–(3.9) for the same values of the parameters as in figure 9.

While solving the axisymmetric three-dimensional problem, we found that the convergence of iterations is very sensitive to the choice of the initial guess. It turned out that one of the reasons for this is the apparent existence of multiple solutions for the same set of governing parameters. Figures 14 and 16 show the example of a second, type 2, solution that was found for the same set of governing parameters as those used to obtain the type 1 solutions shown in figures 10 and 13. The type 1 and 2 solutions appear to be quite similar, as can be seen by comparing figures 10, 12 and 14, 15, but with a few notable distinctions.

First, in the bulk of the layer, the type 2 azimuthal velocity profile has a somewhat larger gradient near the outer cylinder (compare curves 3 in figures 10 a and 14 a), which is expected to promote the hydrodynamic instability detected experimentally in this region.

Figure 16. The same as figure 13 but for the type 2 three-dimensional solution of (3.2)–(3.9).

Figure 17. Toroidal vortex visualisation in an electrolyte layer of thickness $h=5$ mm driven by a radial current of $I_{0}=100$  mA. The stagnation line of the radial flow visible on the free surface as a dark ring was located approximately 3 mm away from the outer electrode.

Second, as seen in figure 16(b), the type 2 solution is characterised by the existence of a secondary counter-rotating vortex in the top corner near the outer cylinder. Comparison of the streamlines azimuthally projected onto the meridional plane, shown for the type 1 and 2 solutions in figures 13(c) and 16(c), indicates that this vortex reflects very different flow topologies in the two cases. While the type 1 solution is characterised by a single-torus flow, the most prominent feature of the type 2 solution is the presence of a second smaller counter-rotating toroidal flow in the upper outer corner of the meridional cross-section. The structure of this axisymmetric flow near the outer cylinder may appear to be similar to that computed for the flow between a pair of shrouded co-rotating disks by Abrahamson, Eaton & Koga (Reference Abrahamson, Eaton and Koga1989). However, there is an important difference: while a pair of symmetric counter-rotating vortices was reported in Abrahamson et al. (Reference Abrahamson, Eaton and Koga1989), the symmetry between the dynamic boundary conditions at the solid bottom and at the free surface is broken and only one toroidal vortex is found in our computations. Such a vortex creates a noticeable radial counterflow along the free surface, with a radial stagnation point some distance away from the outer cylinder; see figure 16(a). The existence of such a toroidal structure was also indicated by experimental observations: the dye tracer released at the free surface away from the side boundaries drifted towards the outer cylinder, confirming the three-dimensional structure of the flow, with a non-zero radial component, but never reached it and remained at a distance of approximately 12 % of the total annular width from the wall; see figure 17. This corresponds to a non-dimensional radial position of $r_{s}\approx 2.62$ , which is very close to the value of $2.60$ for the computed toroidal vortex seen in figure 16. A careful inspection of the experimental photograph in figure 17 also indicates that this toroidal vortex is the locus of instability giving rise to the formation of the anticyclonic vortices seen in figure 1. Thus, these observations suggest that such vortices could be caused by the interplay between the inviscid Rayleigh-type centrifugal instability and the shear instability of a Stewartson-type shear layer (Stewartson Reference Stewartson1957) near the outer cylinder.

Following up the hypothesis that the anticyclonic vortices observed in our experiments could be related to the existence of a secondary single-torus structure, it is of interest to see whether a similarity to other experimental configurations, where such a structure has been seen before, can be established. Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004) detected it in the flow between two counter-rotating disks. In contrast to the symmetric two-torus situation discussed in Abrahamson et al. (Reference Abrahamson, Eaton and Koga1989), in the experiments of Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004), the top–bottom symmetry was broken by a difference in rotation speeds between the two disks (even though, for a different physical reason, such a symmetry breaking also occurs in our flows). Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004) argued that the appearance and location of the instability structures that they observed were related to the existence of a stagnation circle at the slower rotating disk, where the azimuthal vorticity component $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}$ changes sign. Further, Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004) also observed that the location of this stagnation circle was close to the location of the maximum of the vertical vorticity component $\unicode[STIX]{x1D714}_{z}$ . These observations enabled the authors to conclude that, in their experiment, the instability was caused by a classical Kelvin–Helmholtz instability in the shear layer separating the faster rotating outer region from the slower rotating inner one and containing a stagnation circle. Somewhat earlier, Lopez et al. (Reference Lopez, Hart, Marques, Kittelman and Shen2002) considered a similar experimental and computational set-up, but with a larger angular velocity differential. They also found that the instability, somewhat resembling in shape the one observed in our experiments, was that of a free shear layer ‘cutting’ the flow domain into two subregions topologically similar to two tori. Despite the topological similarity to flows discussed by Lopez et al. (Reference Lopez, Hart, Marques, Kittelman and Shen2002) and Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004), our two-torus flows differ from them in a fundamental way: the azimuthal velocities within the two regions have the same sign. More precisely, in our configuration, the location of the stagnation circle on the free surface, where $u_{r}=0$ (it should be noted that, at the free surface next to which the secondary torus is formed, $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}\equiv 0$ , so that it cannot be used to define a boundary between two flow regions, as was done in Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004)), does not coincide with the location of the maximum vorticity component $\unicode[STIX]{x1D714}_{z}$ (which is achieved at the outer cylindrical wall in our case). Therefore, the Kelvin–Helmholtz-type instability responsible for the bulk flow patterns observed by Lopez et al. (Reference Lopez, Hart, Marques, Kittelman and Shen2002) and Moisy et al. (Reference Moisy, Doaré, Pasutto, Daube and Rabaud2004) is not an active physical mechanism in the flows discussed here.

Given the existence of multiple solutions, a continuation procedure, where a converged steady axisymmetric solution obtained for one set of parameters is used as an initial guess for a set of slightly varied parameters, was implemented to trace the parametric evolution of the computed flow. However, this continuation procedure was found to fail in certain regions of the parameter space. This indicates another likely feature of the problem: the existence of bifurcation points where the type of solution changes abruptly. Our computations also indicated that steady solutions might not always exist. At some parametric values, the iterations neither converged nor diverged, with the equation residuals oscillating around some small values. Such a behaviour of steady-state solvers typically occurs when steady solutions bifurcate to time-periodic solutions. For example, by using the continuation procedure, we were able to trace the existence of the type 1 and type 2 solutions depicted in figures 10–16 up to a current value of $I_{0}\approx 0.235$ A (keeping the rest of the physical parameters fixed), but no steady solutions were detected for larger currents. A similar behaviour of the numerical solutions was obtained when the magnitude of the magnetic field was varied instead of the current. These observations provide even stronger motivation for further flow stability studies which will be reported in a separate dedicated publication.

8 Estimation of the free surface deformation

To conclude the discussion of methodological aspects of the study of electrically driven flows in thin cylindrical layers, we have to consider one more important aspect. In all of the computations reported so far, it was assumed that the free surface of the electrolyte layer remains flat. While this is a commonly used assumption (Satijn et al. Reference Satijn, Cense, Verzicco, Clercx and van Heijst2001; Delacroix & Davoust Reference Delacroix and Davoust2014, e.g.) which greatly simplifies modelling of slow flows in thin cylindrical layers, the physical reality is that the free surface starts to deform as soon as fluid is set in motion. A depression develops near the centre and the fluid level rises near the outer boundary (stirred flow in a tea cup is a common everyday analogy). In the experiments that motivated the present study (Pérez-Barrera et al. Reference Pérez-Barrera, Pérez-Espinoza, Ortiz, Ramos and Cuevas2015), this deformation remained small (not visible with the naked eye) in all tested regimes, including the onset of instability and beyond, so that the problem formulation adopted here ignored it completely. However, the obtained numerical results can be used to approximately estimate the shape of the deformed free surface and then obtain a quantitative measure of the error introduced by ignoring the curvature of the free surface.

Taking the divergence of the momentum equations written in a vector form and using the continuity equations for incompressible fluid and magnetic induction, we obtain the following Poisson equation for pressure:

where $\unicode[STIX]{x1D714}^{2}=\unicode[STIX]{x1D714}_{r}^{2}+\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D703}}^{2}+\unicode[STIX]{x1D714}_{z}^{2}$ and $u^{2}=u_{r}^{2}+u_{\unicode[STIX]{x1D703}}^{2}+u_{z}^{2}$ . Given the magnetic field and velocities, the above equation is solved (up to an arbitrary additive constant, which is set to make the pressure at the left lower corner $(r,z)=(\unicode[STIX]{x1D6FC}-1,-1)$ zero) subject to the following boundary conditions:

These boundary conditions are consistent with the momentum equations in the vicinity of the physical boundaries, and take into account the velocity boundary conditions (3.12) and (3.13).

Figure 18. The cross-layer pressure isolines with the subtracted hydrostatic component (a) and the estimated free surface shape (b) for the same values of the parameters as in figure 9. The pressure increases towards the outer cylinder (from left to right in a).

The resulting solution for the pressure distribution is illustrated in figure 18(a). Apart from thin layers adjacent to the vertical boundaries, the pressure increases almost linearly in the radial direction from the inner to the outer cylinder. If the free surface is to remain flat, this distribution is physically impossible, since, in this case, the pressure (and thus the stress normal to the free surface) would have to be discontinuous, given that the pressure just above the fluid layer has to be constant and equal to atmospheric pressure. In reality, this mismatch does not occur because the free surface deforms and the variable hydrostatic pressure component appears and balances the dynamic pressure within the moving fluid with its atmospheric counterpart. To obtain an estimation of the free surface shape $1+h(r)$ , we can then write the non-dimensional hydrostatic balance equation

The integration constant $C$ is determined from the total mass conservation condition

This estimation predicts the qualitatively expected shape of the free surface illustrated in figure 18(b): the fluid surface develops a trough near the inner cylinder and is elevated near the outer vertical boundary, and its slope remains close to zero at both radial extremes. The plot in figure 18 confirms that the vertical deviation of the free surface in the considered regimes does not exceed a few tenths of a per cent. Therefore, neglect of it is not expected to introduce a significant error, let alone change the obtained solutions qualitatively, if the values of $Fr$ remain small, as in the considered examples. However, our numerical experiments show that the issue of properly resolving the free surface deformation quickly becomes pressing and deteriorates the convergence of iterations significantly as the Froude number increases.