1 Introduction
Fluid flows in two spatial dimensions have been the subject of substantial research efforts in recent decades. For the greater part of the twentieth century, it was generally considered that two-dimensional (2D) flows were merely a theoretical idealization with limited practical relevance. This conception has changed drastically since the 1980s, when experiments in thin electrolyte layers (Bondarenko, Gak & Dolzhanskiy Reference Bondarenko, Gak and Dolzhanskiy1979), soap films (Couder Reference Couder1984) and liquid metals (Sommeria & Moreau Reference Sommeria and Moreau1982) demonstrated that nearly 2D flows can indeed be realized in the laboratory. Today, experimental approximations of 2D flows are widely employed as models of atmospheric and oceanic flows (Boffetta & Ecke Reference Boffetta and Ecke2012; Dolzhansky Reference Dolzhansky2013). Being theoretically and experimentally more amenable than their three-dimensional (3D) counterparts, 2D flows have also served as platforms for studying new phenomena such as turbulent cascades (Sommeria Reference Sommeria1986; Tabeling et al. Reference Tabeling, Burkhart, Cardoso and Willaime1991), coherent structures (Sommeria, Meyers & Swinney Reference Sommeria, Meyers and Swinney1988) and mixing (Haller & Yuan Reference Haller and Yuan2000).
Perhaps one of the best-known examples of 2D flows is the one introduced by Andrey Kolmogorov in 1959 as a mathematical problem for studying hydrodynamic stability (Arnold & Meshalkin Reference Arnold and Meshalkin1960). The Kolmogorov flow represents the motion of a viscous fluid in two dimensions (we will refer to these as
$x$
and
$y$
) driven by a forcing that points along the
$x$
-direction and varies sinusoidally in the
$y$
-direction. The fluid flow is considered incompressible,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$
, and is governed by the 2D Navier–Stokes equation,
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}.\end{eqnarray}$$
Here,
$\boldsymbol{u}=(u_{x},u_{y})$
is the velocity field,
$p$
is the 2D pressure field, and
$\boldsymbol{f}=A\sin (\unicode[STIX]{x1D705}y)\hat{\boldsymbol{x}}$
represents the driving force with amplitude
$A$
and wavenumber
$\unicode[STIX]{x1D705}$
. The parameters
$\unicode[STIX]{x1D70C}$
and
$\unicode[STIX]{x1D708}$
are the density and the kinematic viscosity of the fluid, respectively. Kolmogorov flow has served as a convenient model for understanding a wide variety of hydrodynamic phenomena in two dimensions, such as fluid instabilities (Meshalkin & Sinai Reference Meshalkin and Sinai1961; Iudovich Reference Iudovich1965; Kliatskin Reference Kliatskin1972; Nepomniashchii Reference Nepomniashchii1976), 2D turbulence (Green Reference Green1974) and coherent structures (Armbruster et al.
Reference Armbruster, Heiland, Kostelich and Nicolaenko1992; Smaoui Reference Smaoui2001; Chandler & Kerswell Reference Chandler and Kerswell2013).
Practically realizable flows, however, are never strictly 2D. Experimental approximations of Kolmogorov flow have often been carried out either in shallow layers of electrolytes (Bondarenko et al.
Reference Bondarenko, Gak and Dolzhanskiy1979) or in soap films (Burgess et al.
Reference Burgess, Bizon, McCormick, Swift and Swinney1999), wherein geometric confinement suppresses the component of velocity along one of the spatial directions (
$z$
). The remaining two velocity components, however, generally depend on both extended and confined coordinates, making the flow ‘quasi-two-dimensional’ (Q2D). To account for the dependence on the confined coordinate, Q2D flows in shallow layers have often been modelled by adding a linear friction term to the 2D Navier–Stokes equation (1.1):
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D6FC}\boldsymbol{u}+\boldsymbol{f},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}$
is a constant. Here,
$\boldsymbol{u}$
corresponds to the velocity field at the electrolyte–air interface. The addition of this term was first suggested by Bondarenko et al. (Reference Bondarenko, Gak and Dolzhanskiy1979) to model a Q2D flow generated in a homogeneous shallow electrolyte layer. In such a flow, the bottom of the fluid layer is constrained to be at rest because it is in contact with the solid surface of the container holding the fluid. This no-slip constraint at the bottom of the fluid layer causes a gradient in the magnitude of the horizontal velocity along the confined direction
$z$
. Bondarenko et al. (Reference Bondarenko, Gak and Dolzhanskiy1979) rationalized that the dissipation due to this shear, for sufficiently shallow fluid layers, is captured by the linear friction term. In the context of Q2D flows in electrolyte layers, this term has come to be known as ‘Rayleigh friction’. Experimental flows in thin layers, and their 2D approximations employing (1.2), are now commonly referred to as ‘Kolmogorov-like’ when the forcing profile is nearly sinusoidal. Note that linear friction models, similar to that in (1.2), have also been employed to describe Q2D flows in liquid metals (Sommeria Reference Sommeria1986) and soap films (Couder, Chomaz & Rabaud Reference Couder, Chomaz and Rabaud1989; Burgess et al.
Reference Burgess, Bizon, McCormick, Swift and Swinney1999). The motivation behind the addition of the friction term in these models is different from that in (1.2). In this article we are only concerned with flows in shallow electrolyte layers.
Experimental realizations of Q2D flows in recent years have employed two-fluid-layer set-ups: either a set-up with miscible layers made up of a heavy electrolyte fluid (salt water) beneath a lighter non-conducting fluid (pure water) (Marteau, Cardoso & Tabeling Reference Marteau, Cardoso and Tabeling1995; Paret & Tabeling Reference Paret and Tabeling1997; Kelley & Ouellette Reference Kelley and Ouellette2011), or a set-up with immiscible layers made up of a heavy dielectric fluid beneath a lighter electrolyte (Rivera & Ecke Reference Rivera and Ecke2005; Akkermans et al.
Reference Akkermans, Kamp, Clercx and Van Heijst2008, Reference Akkermans, Kamp, Clercx and van Heijst2010). The rationale behind these modifications was that, in addition to confinement, density stratification and immiscibility should enhance two-dimensionality in the top layer. Theoretical models of Q2D experimental flows realized in stratified layers of fluids, however, have not accurately modelled the effect of inhomogeneity in fluid properties as well as the gradient in the magnitude of horizontal velocity
$\boldsymbol{u}(x,y)$
along the confined direction
$z$
. Consequently, experiments were compared with simulations based on the 2D model (1.2) with empirically estimated parameters (Jüttner et al.
Reference Jüttner, Marteau, Tabeling and Thess1997; Boffetta & Ecke Reference Boffetta and Ecke2012).
To address this deficiency, Suri et al. (Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014) have investigated the variation in the horizontal velocity
$\boldsymbol{v}(x,y,z,t)$
along the confined direction
$z$
for a stratified two-immiscible-layer set-up. Following Dovzhenko, Obukhov & Ponomarev (Reference Dovzhenko, Obukhov and Ponomarev1981), the Q2D velocity was approximated as
$$\begin{eqnarray}\boldsymbol{v}(x,y,z,t)=P(z)\boldsymbol{u}(x,y,t)=P(z)[u_{x}(x,y,t)\hat{\boldsymbol{x}}+u_{y}(x,y,t)\hat{\boldsymbol{y}}],\end{eqnarray}$$
where
$\boldsymbol{u}(x,y,t)$
corresponds to the 2D velocity field at the electrolyte–air interface and
$P(z)$
models the variation of the horizontal velocity along
$z$
. By substituting the form of velocity in (1.3) into the 3D Navier–Stokes equation and integrating along the
$z$
-direction, the following modified version of (1.2) was derived:
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\unicode[STIX]{x1D6FD}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\frac{1}{\bar{\unicode[STIX]{x1D70C}}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D6FC}\boldsymbol{u}+\boldsymbol{f}.\end{eqnarray}$$
In the above equation,
$\boldsymbol{f}$
is the depth-averaged force density and the parameters
$\unicode[STIX]{x1D6FD}$
,
$\bar{\unicode[STIX]{x1D70C}}$
,
$\unicode[STIX]{x1D708}$
, and
$\unicode[STIX]{x1D6FC}$
are given by:
$$\begin{eqnarray}\unicode[STIX]{x1D6FD}=\frac{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D70C}P^{2}\,\text{d}z}{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D70C}P\,\text{d}z},\quad \bar{\unicode[STIX]{x1D70C}}=\frac{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D70C}P\,\text{d}z}{h},\quad \unicode[STIX]{x1D708}=\frac{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D707}P\,\text{d}z}{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D70C}P\,\text{d}z},\quad \unicode[STIX]{x1D6FC}=\frac{\unicode[STIX]{x1D707}\displaystyle \left(\frac{\text{d}P}{\text{d}z}\right)_{z=0}}{\displaystyle \int _{0}^{h}\unicode[STIX]{x1D70C}P\,\text{d}z},\end{eqnarray}$$
where
$h$
is the total thickness of the two fluid layers. The vertical profile
$P(z)$
is very weakly dependent on the horizontal flow profile
$\boldsymbol{u}(x,y,t)$
. The profile
$P(z)$
that corresponds to a sinusoidal horizontal flow (described as the straight flow below) was computed and validated against experimental measurements in Suri et al. (Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014).
The prefactor
$\unicode[STIX]{x1D6FD}$
reflects the change in the mean inertia of the fluid layer due to the variation
$P(z)$
of the horizontal velocity in the vertical direction. Since
$P(z)\neq 1$
in experiments,
$\unicode[STIX]{x1D6FD}\neq 1$
, which distinguishes (1.4) from all previous 2D models of flows in shallow electrolyte layers. For multilayer set-ups, the coefficients
$\unicode[STIX]{x1D6FD}$
,
$\bar{\unicode[STIX]{x1D70C}}$
,
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D6FC}$
account for both the inhomogeneity in fluid properties as well the vertical profile
$P(z)$
, as suggested by (1.5). Equations (1.1) and (1.2) can be treated as special cases of (1.4) with suitable choices of the parameters
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
. Furthermore, (1.2) can also be obtained from (1.4) by rescaling the variables (cf. appendix A).
In this article we study instabilities of a Q2D Kolmogorov-like flow realized in a set-up with two immiscible fluid layers and compare experimental results with direct numerical simulations (DNS) of (1.4). Most previous studies of Kolmogorov-like flow that compared experiments with theoretical predictions assumed a perfectly sinusoidal shear flow on an unbounded or periodic domain (Bondarenko et al.
Reference Bondarenko, Gak and Dolzhanskiy1979; Dovzhenko, Krymov & Ponomarev Reference Dovzhenko, Krymov and Ponomarev1984; Batchaev & Ponomarev Reference Batchaev and Ponomarev1989; Krymov Reference Krymov1989; Dolzhanskii, Krymov & Manin Reference Dolzhanskii, Krymov and Manin1992; Thess Reference Thess1992). While some of these studies reported quantitative agreement between theory and experiment with regard to the primary instability, none were able to match simultaneously both the critical Reynolds number and the critical wavenumber. Even matching one of these required treating the Rayleigh friction coefficient
$\unicode[STIX]{x1D6FC}$
as an adjustable parameter. To address these shortcomings, we have performed a systematic investigation of the effects of lateral confinement using numerical simulations with three different sets of boundary conditions. Furthermore, we investigate how the observed flow patterns and their stability are affected by the deviations in the forcing profile from perfect periodicity in the extended directions and by the variation of the forcing profile in the confined direction. Finally, we compare the results of numerical simulations with the experimental observations for the secondary instability, which introduces time dependence into the flow.
Given that none of the previous models of Q2D flows were quantitatively accurate, the availability of an experimental set-up and a matching 2D model that are in quantitative agreement is quite important for a number of reasons. In particular, this allows us to make substantial progress (Suri et al. Reference Suri, Tithof, Grigoriev and Schatz2017) in understanding the role of coherent structures in turbulent flows (Hussain Reference Hussain1986; Kawahara, Uhlmann & van Veen Reference Kawahara, Uhlmann and van Veen2012; Chandler & Kerswell Reference Chandler and Kerswell2013; Gallet & Young Reference Gallet and Young2013; Haller Reference Haller2015). Recent advances in transitional flows and weak turbulence rely on a deterministic, geometrical description where the evolution of the flow is guided by non-chaotic, unstable solutions of the Navier–Stokes equation, often referred to as exact coherent structures (ECS) (Nagata Reference Nagata1997; Waleffe Reference Waleffe1998; Kerswell Reference Kerswell2005; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007; Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2009). The bulk of numerical studies have explored the role of ECS in 3D flows simulated on periodic domains with simple geometries, such as pipe flow, plane Couette flow and plane Poiseuille flow.
However, experimental evidence for the role of ECS in 3D flows has been scarce (Hof et al. Reference Hof, van Doorne, Westerweel, Nieuwstadt, Faisst, Eckhardt, Wedin, Kerswell and Waleffe2004; de Lozar et al. Reference de Lozar, Mellibovsky, Avila and Hof2012; Dennis & Sogaro Reference Dennis and Sogaro2014), in part due to technical limitations in obtaining spatially and temporally resolved 3D velocity fields. Q2D flows, on the other hand, can be quantified using 2D planar velocity fields, which are relatively easy to measure. Recently, Chandler & Kerswell (Reference Chandler and Kerswell2013) and Lucas & Kerswell (Reference Lucas and Kerswell2014, Reference Lucas and Kerswell2015) have identified dozens of ECS in numerical simulations of a weakly turbulent 2D Kolmogorov flow, governed by (1.1) with periodic boundary conditions, which, however, do not describe flows that can be realized in experiments (Bondarenko et al. Reference Bondarenko, Gak and Dolzhanskiy1979; Dolzhanskii et al. Reference Dolzhanskii, Krymov and Manin1992; Suri et al. Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). Hence, the analysis presented herein should provide the much needed foundation for further studies of 2D turbulence that focus on experimental validation of theoretical predictions, building on the results of Suri et al. (Reference Suri, Tithof, Grigoriev and Schatz2017).
This article is organized as follows. In § 2, we describe the experimental set-up employed to generate a Q2D Kolmogorov-like flow. In § 3, we introduce a realistic model of the forcing in the experiment and discuss different types of lateral boundary conditions, which are used to study the effects of confinement theoretically. In § 4, we compare the flow fields obtained from experimental measurements with those from the numerical simulations for different flow regimes and characterize the bifurcations associated with increasing the forcing strength. In § 5, we discuss how the nature of the primary instability depends on the lateral boundary conditions. Conclusions are presented in § 6.
2 Experimental set-up
We generate a Q2D Kolmogorov-like flow in the experiment using a stratified set-up with two immiscible fluid layers, first introduced by Rivera & Ecke (Reference Rivera and Ecke2005). In this configuration, a lighter electrolyte is suspended on top of a denser dielectric, which serves as a lubricant between the electrolyte layer and the solid surface at the bottom of the container that holds the fluids. The fluid layers are set in motion using Lorentz forces resulting from the interaction of a direct current passing through the electrolyte and a spatially varying magnetic field.
We use a magnet array consisting of 14 NdFeB magnets (Grade N42) to generate a magnetic field that varies roughly sinusoidally along one direction, approximating the forcing in the Kolmogorov flow. Each magnet in the array is 15.24 cm long and 1.27 cm wide, with a thickness of
$0.32\pm 0.01~\text{cm}$
. The magnetization is parallel to the thickness dimension, with a surface field strength of approximately 0.2 T. The magnets are positioned side-by-side along their width to form a
$15.24~\text{cm}\times (14\times 1.27~\text{cm})\times 0.32~\text{cm}$
array such that the adjacent magnets have fields pointing in opposite directions, normal to the plane of the array. This magnet array is placed on a flat aluminium plate of dimensions
$30.5~\text{cm}\times 30.5~\text{cm}\times 1.0~\text{cm}$
, and rectangular pieces of aluminium with the same thickness as the magnets
$(0.32\pm 0.02~\text{cm})$
are placed beside the magnet array to create a level surface. Manufacturing imperfections in the individual magnets and the aluminium siding result in a surface that is not adequately smooth. Hence, a thin glass plate measuring
$25.4~\text{cm}\times 25.4~\text{cm}$
in area with a thickness of
$0.079\pm 0.005~\text{cm}$
is placed atop the magnets and siding to provide a uniform surface. A thin layer of black, adhesive contact paper (with approximate thickness 0.005 cm) is placed on top of the glass plate to serve as a dark background for imaging. The surface of the contact paper serves as the bottom boundary for the fluids. We place the origin of our coordinate system at this height and the lateral centre of the magnet array, with the
$x$
-coordinate aligned with the magnets’ longest side, the
$y$
-coordinate pointing in the direction of the magnet array periodicity, and the
$z$
-coordinate in the vertical direction. A schematic diagram is shown in figure 1.
Figure 1. A schematic diagram of the two-immiscible-layer experimental set-up for generating Kolmogorov-like flow viewed (a) from above and (b) from the side. The vectors
$\boldsymbol{J}$
,
$\boldsymbol{B}$
and
$\boldsymbol{F}$
denote, respectively, the directions of the electric current, the magnetic field and the resulting Lorentz force. The flow is bounded by two endwalls, two sidewalls (electrodes) and a no-slip bottom surface, while the top surface is a free electrolyte–air interface. This container is mounted on an aluminium plate that is levelled and submerged in a water bath that is temperature-regulated such that the electrolyte is maintained at
$23.0\pm 0.2\,^{\circ }\text{C}$
.
Rectangular bars of acrylic are affixed directly onto the contact paper to create the lateral boundaries of the container that will hold the fluids. Parallel to the
$y$
-direction, two bars are placed at a distance of 17.8 cm apart, centred about the origin. These solid boundaries for the fluid are henceforth referred to as the ‘endwalls’. Similarly, running parallel to the
$x$
-direction, two electrodes mounted on rectangular bars of acrylic are placed at a distance of 22.9 cm, symmetrically relative to the origin. These boundaries are henceforth referred to as the ‘sidewalls’ and are used to drive the current through the electrolyte. The placement of the endwalls and sidewalls leaves a buffer region of
$d_{x}=1.3~\text{cm}$
and
$d_{y}=2.5~\text{cm}$
, respectively, between the edge of the magnet array and these solid boundaries.
The aluminium plate upon which the magnets are mounted is supported by three screws, which are adjusted to level the system. The interior of the container is filled with
$122\pm 4~\text{ml}$
of a dielectric fluid and
$122\pm 2~\text{ml}$
of an electrolyte to form two immiscible layers that are
$0.30\pm 0.01~\text{cm}$
and
$0.30\pm 0.005~\text{cm}$
thick, respectively. The dielectric fluid used is perfluorooctane, which has a viscosity of
$\unicode[STIX]{x1D707}_{d}=1.30~\text{mPa}~\text{s}$
and a density of
$\unicode[STIX]{x1D70C}_{d}=1769~\text{kg}~\text{m}^{-3}$
at
$23.0\,^{\circ }\text{C}$
. The electrolyte fluid is a solution consisting of 60 % 1 M copper sulfate solution and 40 % glycerol by weight. The electrolyte’s viscosity is
$\unicode[STIX]{x1D707}_{c}=5.85~\text{mPa}~\text{s}$
and the density is
$\unicode[STIX]{x1D70C}_{c}=1192~\text{kg}~\text{m}^{-3}$
at
$23.0\,^{\circ }\text{C}$
. Note that a large viscosity ratio
$\unicode[STIX]{x1D707}_{c}/\unicode[STIX]{x1D707}_{d}=4.5$
has been chosen to enhance the two-dimensionality of the electrolyte, as described by Suri et al. (Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). A small amount of viscosity-matched surfactant is added to the electrolyte to lower the surface tension, and a glass lid is placed on top of the container to limit evaporation.
A direct current, which serves as the control parameter, is then passed through the electrolyte; the resulting current density
$J$
ranges from approximately 10 to
$40~\text{A}~\text{m}^{-2}$
across the different runs. The interaction of this current with the spatially alternating magnetic field
$\boldsymbol{B}$
results in a spatially alternating Lorentz force
$\boldsymbol{F}$
, which drives the electrolyte (cf. figure 1
a). The viscous coupling between the electrolyte and the dielectric fluids sets the dielectric fluid in motion as well. Since passing a current through a resistive conductor (the electrolyte) results in Joule heating, a calibrated thermistor is placed in the corner of the fluid domain to monitor the fluid temperature, and the aluminium plate is immersed in a temperature-controlled water bath. The water bath is regulated such that the temperature of the electrolyte is maintained to
$23.0\pm 0.2\,^{\circ }\text{C}$
. By limiting the temperature fluctuations, the associated change in viscosity of the fluids is kept to a minimum.
For flow visualization, we add hollow glass microspheres (Glass Bubbles K15, manufactured by 3M), sieved to obtain particles with mean radius
$24.5\pm 2~\unicode[STIX]{x03BC}\text{m}$
and mean density
$150~\text{kg}~\text{m}^{-3}$
. Being lighter than the electrolyte, the microspheres stay afloat at the electrolyte–air interface for the duration of the experiment. The microspheres are illuminated with white light-emitting diodes placed near the endwalls, outside the container holding the fluids. The flow is imaged at 15 Hz with a camera (DMK 31BU03, manufactured by The Imaging Source) placed directly above the set-up. This camera has a charge-coupled device (CCD) sensor with a resolution of
$1024\times 768$
pixels, which results in an adequate resolution of approximately 53 pixels per magnet width. The flow velocities are calculated using the PRANA particle image velocimetry (PIV) package (Eckstein & Vlachos Reference Eckstein and Vlachos2009; Drew, Charonko & Vlachos Reference Drew, Charonko and Vlachos2013). This software employs a multigrid PIV algorithm that deforms images to better resolve flows with high shear. The velocity field is resolved on a
$169\times 126$
grid, with approximately 9 points per magnet width.
For the experimental measurements listed above, we obtain the following depth-averaged values for the parameters in (1.4):
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
,
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
and
$\bar{\unicode[STIX]{x1D70C}}=976~\text{kg}~\text{m}^{-3}$
. These parameters were computed using the vertical profile
$P(z)$
that corresponds to the strictly sinusoidal flow (Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). The complexity of the flow in both the experiment and simulation is characterized by the Reynolds number, which we define as
$$\begin{eqnarray}Re=\frac{Uw}{\unicode[STIX]{x1D708}},\end{eqnarray}$$
where
$w=1.27~\text{cm}$
is the width of one magnet and
$U=\sqrt{\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{ u}\rangle }$
is the measured root-mean-square (r.m.s.) velocity, where
$\langle \cdot \rangle$
denotes spatial averaging over a subregion with dimensions
$10.16~\text{cm}\times 10.16~\text{cm}$
at the centre of the magnet array. Note that, from this point on, unless noted otherwise, the characteristic length scale
$w$
, velocity scale
$U$
and time scale
$w/U$
will be used for non-dimensionalization.
3 Numerical modelling
In this section we present a model of the magnetic field generated by the finite array of permanent magnets in the experiment. We then introduce three types of boundary conditions used in our numerical simulations of the flow.
3.1 Modelling the magnetic field
In the discussion so far, we have not addressed an important question of how the 2D forcing function
$\boldsymbol{f}$
in (1.4) relates to the 3D forcing
$\boldsymbol{F}$
in the experiment. For a 2D Kolmogorov flow, the forcing
$\boldsymbol{f}$
is sinusoidal, by definition. However, for Kolmogorov-like flows realized in electromagnetically driven shallow layers of electrolyte,
$\boldsymbol{f}$
needs to be computed from the 3D Lorentz force
$\boldsymbol{F}$
arising from the interaction of the magnetic field
$\boldsymbol{B}$
produced by a finite magnet array with a current density
$\boldsymbol{J}$
(Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). The current density is easily calculated from geometrical considerations, but the magnetic field generated by the array of permanent magnets is quite complicated. For
$\boldsymbol{J}=J\hat{\boldsymbol{y}}$
, the Lorentz force density at any location
$(x,y,z)$
within the electrolyte layer is given by
$\boldsymbol{F}=\boldsymbol{J}\times \boldsymbol{B}=JB_{z}\hat{\boldsymbol{x}}-JB_{x}\hat{\boldsymbol{z}}$
. Here,
$B_{x}$
and
$B_{z}$
are the
$x$
- and
$z$
-components of the magnetic field, respectively, which vary along all three coordinates
$x$
,
$y$
and
$z$
. Experimental measurements show that the typical value of
$B_{x}$
is less than 3 % of the value of
$B_{z}$
at any given location within the electrolyte. Furthermore, the vertical component of the Lorentz force along with the gravitational force will be balanced by the vertical gradient of the pressure. Hence, the Lorentz force density for all practical purposes can be approximated as
$\boldsymbol{F}\approx JB_{z}\hat{\boldsymbol{x}}$
. One can then compute
$\boldsymbol{f}$
using the expression
$$\begin{eqnarray}\boldsymbol{f}=\frac{1}{\bar{\unicode[STIX]{x1D70C}}}\int _{h_{d}}^{h_{d}+h_{e}}\frac{JB_{z}(x,y,z)\,\text{d}z}{h_{d}+h_{e}}\hat{\boldsymbol{x}},\end{eqnarray}$$
where
$h_{e}$
and
$h_{d}$
are the thicknesses of the electrolyte and dielectric layers, respectively.
Figure 2. The
$z$
-component of the magnetic field
$B_{z}$
(a) at the longitudinal centre of the domain (
$x=0$
) and (b) along the magnet centrelines at
$y=\pm \{0.5,1.5,2.5,3.5,4.5,5.5,6.5\}$
. In (a), the experimental measurements at a height
$z=0.265$
(just above the dielectric–electrolyte interface) and at
$z=0.438$
(just below the electrolyte free surface) are shown, respectively, as open squares and filled circles. In (b), the black symbols indicate the experimental measurements at a height
$z=0.438$
along the magnet centrelines. A least-squares fit has been performed using the data in (a) to determine the scaling factor for the dipole summation; the scaled dipole summation magnetic field is shown as the red full lines. The experimental uncertainties are the size of the symbols or smaller.
The black symbols in figure 2(a) show the experimental measurements of
$B_{z}$
along the line
$x=0$
, passing above the centre of the magnet array at two different heights. Clearly, the magnetic field profile deviates significantly from that of a pure sinusoid. Furthermore, one cannot ignore the fringe fields near the edges of the array. To obtain a magnetic field profile that closely resembles the one in the experiment, one could measure the
$z$
-component of the magnetic field (
$B_{z}$
) across the entire flow domain at various heights above the magnet array. Using the measured field, one could then compute the depth-averaged forcing profile using (3.1) (Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). However, since measuring
$B_{z}$
on a 3D grid is an extremely tedious process, we circumvent the labour by numerically modelling the magnet array as described below.
The magnets in the array are arranged such that adjacent ones have magnetization pointing in opposite directions, along
$\pm \hat{\boldsymbol{z}}$
. To obtain a magnetic field that closely resembles the one due to this array, we model each magnet as a uniformly magnetized medium, i.e. as a 3D cubic lattice of identical dipoles, each with a moment
$m\hat{\boldsymbol{z}}$
. Changing the sign of
$m$
across adjacent magnets accounts for the alternating direction of magnetization. The magnetic field at any location
$(x,y,z)$
above the array is then approximated using the linear superposition of the field contribution from all of the dipoles modelling the array. Hence, we refer to this model as the ‘dipole summation’. Since the strength of the dipole
$m$
cannot be measured experimentally, a single scaling parameter is calculated from a least-squares fit with the experimental measurements, taken at two heights. The rescaled dipole summation magnetic field is shown in figure 2(a) (red lines), along with the experimental measurements of
$B_{z}$
(black symbols), corresponding to the line
$x=0$
at heights
$z=0.265$
and
$z=0.438$
. Figure 2(b) shows the magnetic field comparison at
$z=0.438$
along the magnet centrelines. Note that the electrolyte layer in the experiment is bounded by the planes
$z=0.236$
and
$z=0.472$
. Hence, we compute the magnetic field
$B_{z}(x,y,z)$
using the dipole summation at various heights, in steps of
$0.0197$
, in the region
$0.236<z<0.472$
and depth-average it using a discrete version of the expression (3.1).
3.2 Boundary conditions for direct numerical simulations
In the experimental Kolmogorov-like flow, vertical solid walls serve as the lateral boundaries, resulting in a no-slip boundary condition for the velocity. However, for reasons of analytical and computational feasibility, Kolmogorov flow has been studied almost exclusively using unbounded or periodic domains. Neither an infinite lateral extent nor periodicity offer a realistic representation of the effect of boundary conditions in the experiment, as far as the flow’s structure and its stability are concerned. To explore the role of boundaries, we compare the experiment to numerical simulations using computational domains with increasing degrees of confinement. The three different computational domains we study are described below.
-
(i) Doubly periodic domain: This computational domain is chosen to coincide with the central
$8w\times 8w$
region of the experimental domain (
$|x|\leqslant 4$
and
$|y|\leqslant 4$
in non-dimensional units). The simulated flow is constrained to be periodic in both the longitudinal and transverse directions, i.e.
$\boldsymbol{u}(x=-4,y)=\boldsymbol{u}(x=4,y)$
and
$\boldsymbol{u}(x,y=-4)=\boldsymbol{u}(x,y=4)$
. Along the transverse direction it spans a width equalling that of eight magnets. The 2D forcing profile
$\boldsymbol{f}=f_{x}(y)\hat{\boldsymbol{x}}$
over this doubly periodic domain is constructed from the depth-averaged magnetic field presented in § 3.1 by retaining only the two dominant Fourier modes,
$\sin (\unicode[STIX]{x1D705}y)$
and
$\sin (3\unicode[STIX]{x1D705}y)$
, along the
$y$
-direction, where
$\unicode[STIX]{x1D705}=\unicode[STIX]{x03C0}$
in dimensionless units. Along the
$x$
-direction the profile is uniform:
$f_{x}(y)=1.05\sin (\unicode[STIX]{x1D705}y)+0.05\sin (3\unicode[STIX]{x1D705}y)$
. -
(ii) Singly periodic domain: This computational domain coincides with the region
$|x|\leqslant 7$
and
$|y|\leqslant 4$
. The longitudinal dimension is the same as that of the experiment, while the transverse one spans a width equalling that of eight magnets, like in the doubly periodic domain. No-slip boundary conditions are imposed at the endwalls, i.e.
$\boldsymbol{u}(x=\pm 7,y)=0$
, while periodic boundary conditions are imposed along the transverse direction, i.e.
$\boldsymbol{u}(x,y=4)=\boldsymbol{u}(x,y=-4)$
. The 2D forcing profile
$\boldsymbol{f}=f_{x}(x,y)\hat{\boldsymbol{x}}$
over this singly periodic domain is constructed as a product of two one-dimensional (1D) profiles
$f_{x}(x,y)=\unicode[STIX]{x1D712}(x)\unicode[STIX]{x1D713}(y)$
. Along the
$y$
-direction the profile is once again constructed by retaining only two dominant Fourier modes of the depth-averaged magnetic field,
$\unicode[STIX]{x1D713}(y)=1.05\sin (\unicode[STIX]{x1D705}y)+0.05\sin (3\unicode[STIX]{x1D705}y)$
. Along the
$x$
-direction the profile
$\unicode[STIX]{x1D712}(x)$
is chosen to be the depth-averaged magnetic field profile from the dipole summation along the magnet centreline
$y=0.5$
. We note that the effect of transverse confinement has been studied by Thess (Reference Thess1992), and therefore is not investigated here separately. -
(iii) Non-periodic domain: This computational domain is identical to the experimental one in both lateral dimensions, i.e.
$|x|\leqslant 7$
and
$|y|\leqslant 9$
, with no-slip boundary conditions imposed at both the endwalls and sidewalls, i.e.
$\boldsymbol{u}(x=\pm 7,y)=0$
and
$\boldsymbol{u}(x,y=\pm 9)=0$
. As mentioned in § 3.1, the forcing over this domain is computed by depth-averaging the dipole summation.
To compare the experimental observations with those predicted by (1.4) with the three types of boundary conditions described above, we have performed DNS. The flow over the doubly periodic domain is simulated using a pseudo-spectral method in the vorticity–streamfunction formulation, as described in Mitchell (Reference Mitchell2013). This simulation is henceforth referred to as the ‘doubly periodic simulation’ (DPS). For the singly periodic and the non-periodic domains, numerical simulations have been performed using a finite-difference scheme, described in Armfield & Street (Reference Armfield and Street1999). These simulations are hereafter referred to as the ‘singly periodic simulation’ (SPS) and the ‘non-periodic simulation’ (NPS), respectively. Details of the spatiotemporal discretization and the integration schemes employed in all the numerical simulations can be found in appendix B.
4 Comparison of experiment and simulations
In this section we present the results of our comparison between the experiment and the numerical simulations on the three domains described above. First, we discuss the straight uniform flow found at lower Reynolds numbers, with a special emphasis on the effect of boundaries. We then perform linear stability analysis to demonstrate that (1.4) describes the primary instability more accurately than (1.2) on an unbounded domain, but still substantially underpredicts the experimentally observed critical Reynolds number. Next, we describe and compare the steady flow states found in the experiments and simulations above the primary instability. Finally, we discuss the secondary instability that gives rise to a time-periodic flow.
4.1 Straight flow
For weak driving, the flow mimics the forcing closely, with spatially alternating bands of fluid flow along the
$\pm x$
-directions, as can be seen in figure 3 for
$Re=8.1$
. In this figure, black vectors represent the velocity field
$\boldsymbol{u}$
and the colour indicates the vorticity
$\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D735}\times \boldsymbol{u})\boldsymbol{\cdot }\hat{\boldsymbol{z}}$
. For the experiment (figure 3
d), the
$y$
-component of the velocity measured in the central region of the domain is close to zero. However, there are regions of strong recirculation near the endwalls, characterized by a non-zero
$y$
-component of velocity. A closer inspection of the flow shows a slight tilt in the alignment of the flow bands. This tilt is due to the global circulation, resulting from confinement and the fluid flowing in opposite directions over the end magnets at
$y=\pm 6.5$
. Figure 3(a,b) shows the straight flows found in the DPS and SPS. It can be seen that flow fields in the DPS and SPS reproduce the experimental flow qualitatively away from the lateral walls. Furthermore, the SPS captures the turnaround flow near the endwalls. However, neither the SPS nor the DPS displays the tilt of the flow bands observed in the experiment, since the periodic flows are devoid of global circulation. In contrast, the NPS generates a flow field that looks indistinguishable from the experimental one (cf. figure 3
c).
Figure 3. Straight flow fields at
$Re=8.1$
with
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
for the (a) DPS, (b) SPS, (c) NPS and (d) experiment. The dashed lines in (d) indicate the locations of velocity profiles in the experiment that are compared to the simulations. The vorticity colour scale plotted for (a) also applies to (b–d). The velocity vectors are downsampled in each direction by a factor of 8 for the simulations and 4 for the experiment.
For a quantitative description of the straight flow profile, we have plotted in figure 4(a) the longitudinal component
$u_{x}^{exp}$
of the velocity along the line
$x=0$
in the experiment. The location of this cross-section is indicated by the vertical dashed line in figure 3(d). The difference in
$u_{x}$
between the experiment and the numerical simulations along this line is shown in figure 4(b). As can be seen, the DPS and SPS, which are only defined for
$|y|\leqslant 4$
, show systematic deviation from the experiment as high as 18 %, since they do not capture global circulation. In comparison, the NPS agrees to within approximately 5 % over the same region, with no clear systematic deviation. The disagreement between the experiment and NPS in this region, we believe, is a result of the dipole summation not accounting for the variation in the strength of each individual magnet. Closer to the boundaries, at
$y\approx 7$
and
$y\approx -6$
, the largest difference between the NPS and the experiment is around 12 %. The origin of this error is quite subtle and we shall defer its analysis to appendix C.
Figure 4. Profiles of the longitudinal velocity and longitudinal velocity differences at
$Re=8.1$
with
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
. (a) Plot of
$u_{x}^{exp}$
as a function of
$y$
at the longitudinal centre (
$x=0$
). (b) The difference between the longitudinal velocity in the simulations and the experiment,
$u_{x}^{sim}-u_{x}^{exp}$
, as a function of
$y$
at the longitudinal centre (
$x=0$
). Note that the curves corresponding to the DPS and SPS are virtually indistinguishable. (c) Plot of
$u_{x}^{exp}$
as a function of
$x$
at the centreline of a middle magnet (
$y=-0.5$
). (d) The difference between the longitudinal velocity of the simulations and the experiment,
$u_{x}^{sim}-u_{x}^{exp}$
, as a function of
$x$
at the centreline of a middle magnet (
$y=-0.5$
). Note that the curves corresponding to the DPS and SPS are virtually indistinguishable in the region
$-4<x<4$
, where the DPS is defined. Experimental uncertainties are the size of the symbols or smaller.
The experimental longitudinal velocity component
$u_{x}^{exp}$
at
$y=-0.5$
(along a central magnet centreline) is shown in figure 4(c). The very slight asymmetry in the longitudinal velocity is a result of the global circulation. In contrast, the flow in the DPS is perfectly uniform and thus does not capture this asymmetry, as can be seen from the plot of its difference with the experimental profile in figure 4(d). The SPS, which is defined all the way to the endwalls, also does not capture this asymmetry due to the lack of global circulation. The NPS produces the closest agreement: the corresponding flow displays the asymmetry observed in the experiment, with no significant systematic deviation. In summary, the NPS succeeds in capturing the effects of confinement in the experiment with good accuracy, while the DPS and SPS show significant systematic deviations.
4.2 Linear stability analysis of the straight flow
As the strength of the forcing increases, the flow in the experiment undergoes a qualitative change at
$Re_{c}=11.07\pm 0.05$
, with uniform flow bands (cf. figure 4
c) developing modulation that eventually gives rise to distinct stationary vortices. Hence, we shall refer to this flow as the ‘modulated flow’. Several previous experimental studies have reported this transition and have characterized it using the critical Reynolds number (
$Re_{c}^{exp}$
) and wavenumber (
$k_{c}^{exp}$
) of the modulation (Bondarenko et al.
Reference Bondarenko, Gak and Dolzhanskiy1979; Batchaev & Dowzhenko Reference Batchaev and Dowzhenko1983; Obukhov Reference Obukhov1983). In our experiments, the wavenumber just above this transition was measured to be
$k_{c}^{exp}=0.50\unicode[STIX]{x1D705}$
, where
$\unicode[STIX]{x1D705}$
is the wavenumber associated with the forcing. In virtually all previous studies, theoretical estimates for these critical parameters have been obtained by using (1.2) and modelling the straight flow in experiment as a strict sinusoid
$\boldsymbol{u}_{s}\propto \sin (\unicode[STIX]{x1D705}y)$
; the flow stability is then analysed with respect to perturbations
$\unicode[STIX]{x1D6FF}\boldsymbol{u}(y)\text{e}^{\text{i}kx}$
in the transverse component of the velocity. In this section, we revisit this analytical approach for (1.4) to provide estimates for the critical parameters. Many previous studies used a different non-dimensionalization, which corresponds to setting the non-dimensional forcing wavenumber
$\unicode[STIX]{x1D705}$
to unity. To make comparison easier, we will introduce a scaled wavenumber
$q=k/\unicode[STIX]{x1D705}$
, which corresponds to the convention used in those studies.
The strictly sinusoidal straight flow governed by (1.4) on an unbounded domain becomes unstable with respect to perturbations with wavenumber
$q$
above the Reynolds number
$Re=Re_{n}(q)$
, which to a very good accuracy is given by
$$\begin{eqnarray}Re_{n}(q)=\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FD}}\frac{1}{q}\sqrt{\frac{(1+q^{2})}{(1-q^{2})}\left(q^{2}+\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}}\right)\left(1+q^{2}+\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}^{2}}\right)}.\end{eqnarray}$$
This expression was computed by linearizing (1.4) around
$\boldsymbol{u}_{s}$
and calculating its stability with respect to perturbations including three dominant modes,
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}\boldsymbol{u}(y)\text{e}^{\text{i}kx}=\mathop{\sum }_{n=-1,0,1}\unicode[STIX]{x1D716}_{n}\text{e}^{\text{i}\unicode[STIX]{x1D705}(ny+qx)}.\end{eqnarray}$$
A detailed discussion of the stability analysis and the analytical expression for the neutral stability curve, similar in form to that in (4.1), can be found in appendix B of Dolzhansky (Reference Dolzhansky2013). The critical Reynolds number
$Re_{c}=\min _{q}Re_{n}(q)$
and the corresponding critical wavenumber
$k_{c}=\unicode[STIX]{x1D705}q_{c}$
computed using (4.1) can be compared with experimental observations. It is worth noting that this linear stability analysis is based on a purely sinusoidal forcing profile that does not contain any harmonics. Hence, this is one source of discrepancy between these analytical results and the experiment and simulations discussed below, whose forcing profiles are more complicated.
Figure 5. Neutral stability curves (4.1) describing the primary instability. The blue dashed line corresponds to
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
and
$\unicode[STIX]{x1D6FD}=0.83$
, while the green dot-dashed line corresponds to
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
and
$\unicode[STIX]{x1D6FD}=1.0$
. The measurement from the experiment (NPS) is plotted as a black dot (red square); note that the uncertainties in
$Re_{c}$
are smaller than the size of the symbols. In all the cases,
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
is held constant.
The neutral stability curve (blue dashed line) which corresponds to the experimental values of parameters
$\unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D6FD}$
and
$\unicode[STIX]{x1D708}$
is shown in figure 5. The minimum of this neutral stability curve yields a critical Reynolds number
$Re_{c}=9.16$
and an associated critical wavenumber
$q_{c}=0.465$
. The black dot on the plot indicates the critical values
$Re_{c}^{exp}=11.07$
and
$q_{c}^{exp}=0.50$
, corresponding to the instability we observe in the experiment. The relative difference
$(Re_{c}^{exp}-Re_{c})/Re_{c}^{exp}$
between the theoretical estimate for the critical Reynolds number and that measured in experiment is approximately 17 %. The critical wavenumber, however, is in better agreement with the experimentally measured one, with a 7 % relative error.
While the critical Reynolds number obtained from the linear stability analysis clearly disagrees with the experimentally observed one, it is still a significant improvement over analytical estimates for a flow modelled using (1.2), which corresponds to setting
$\unicode[STIX]{x1D6FD}=1$
in (1.4). The corresponding neutral stability curve is indicated by the green dot-dashed line in figure 5. From (4.1) it can be seen that the entire neutral stability curve scales as
$1/\unicode[STIX]{x1D6FD}$
. This implies that the critical wavenumber (
$q_{c}=0.465$
) is independent of
$\unicode[STIX]{x1D6FD}$
, while the predicted critical Reynolds number for
$\unicode[STIX]{x1D6FD}=1$
is
$Re_{c}=7.60$
. This is a 31 % discrepancy with the experimental value, which is comparable to the 30 % discrepancy reported by Bondarenko et al. (Reference Bondarenko, Gak and Dolzhanskiy1979) in a study based on (1.2).
As we discussed previously, the parameter
$\unicode[STIX]{x1D6FD}$
describes the effect of the vertical variation in the magnitude of the horizontal velocity on the effective inertia and nonlinearity of the flow. Equation (1.2) does not account for this effect, so it is natural that its predictions are substantially less accurate.
4.3 Modulated flow
Figure 6(a–d) shows the modulated flow fields corresponding to the DPS, SPS, NPS and experiment, respectively, at
$Re=14$
. At this Reynolds number, the modulated flow is well developed and is visually quite distinct from the straight flow. The anticlockwise global circulation in the experiment strongly affects the alignment of the vortices (see figure 6
d), as can be seen by comparing the modulated flows in the DPS and SPS with the relevant regions of the experimental flow. Unlike the DPS and SPS, the flow field in the NPS captures the features observed in the experiment remarkably well. This unambiguously demonstrates the importance of properly modelling the confinement effects in both the longitudinal and the transverse direction to reproduce the features of the flow in the experiment.
Figure 6. Modulated flow fields at
$Re=14$
with
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
for the (a) DPS, (b) SPS, (c) NPS and (d) experiment. The vorticity colour scale plotted for (a) also applies to (b–d). The velocity vectors are downsampled in each direction by a factor of 8 for the simulations and a factor of 4 for the experiment.
The onset of the modulated flow is characterized by the appearance of the transverse component
$u_{y}$
of the velocity throughout the flow domain. As the driving is increased, the magnitude of
$u_{y}$
also increases. A bifurcation diagram characterizing the transition from the straight to the modulated flow is shown in figure 7(a). We use the spatial mean-square transverse velocity,
$\langle u_{y}^{2}\rangle$
, as the order parameter and plot it as a function of
$Re$
. The spatial average is computed over the central region
$|x|\leqslant 4$
and
$|y|\leqslant 4$
for all simulations and experiment. In comparison to the experimental value of
$Re_{c}^{exp}=11.07$
, the primary instability in the DPS and SPS occurs at much lower Reynolds numbers,
$Re_{c}=9.39$
and
$Re_{c}=9.53$
, respectively. In contrast, by imposing the correct (no-slip) boundary conditions in both the longitudinal and transverse directions, in addition to using a realistic model of the magnetic field, the transition can be predicted quite accurately. The straight to modulated transition in the NPS occurs at
$Re_{c}=10.49$
(red square in figure 5), which is within 5.2 % of
$Re_{c}^{exp}$
. Finally, we note that setting
$\unicode[STIX]{x1D6FD}=1$
results in a poor prediction
$Re_{c}=8.71$
even in the NPS, which corresponds to a 21 % error.
Figure 7. Primary instability for
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
. (a) A bifurcation diagram and (b) the average wavelength of the pattern,
$\bar{\unicode[STIX]{x1D706}}_{x}$
, as a function of
$Re$
for the modulated flow regime. At each
$Re$
, wavelength measurements are made in the central region
$|y|\leqslant 4$
, then averaged; the uncertainty bars indicate one standard deviation in the spatial measurements.
Given that the pattern of vortices observed in the experiment lacks perfect periodicity, we compute the average longitudinal wavelength
$\bar{\unicode[STIX]{x1D706}}_{x}$
using the spatial average of the separation between adjacent vortex centres in the central region
$|y|\leqslant 4$
; the vortex centres are identified by locating local minima in the velocity magnitude. Just above onset, the vortices in the experiment form a lattice with a fairly uniform separation,
$\bar{\unicode[STIX]{x1D706}}_{x}^{exp}\approx 4.0$
(
$q_{c}^{exp}=0.50$
). As the forcing is increased, the mean separation between the vortices increases, as can be seen from the plot of
$\bar{\unicode[STIX]{x1D706}}_{x}$
versus
$Re-Re_{c}$
shown in figure 7(b). Additionally, the vortex lattice becomes spatially irregular, as can be seen in figure 6(d). This spatial variation is quantified in the plot in figure 7(b) wherein the uncertainty bars indicate one standard deviation in the spatial variation of the separation between adjacent vortices. Note that, immediately above onset, accurate identification of the vortex centres in the experiment is not possible because of the very weak modulation; hence, experimental measurements are only plotted for
$Re-Re_{c}\geqslant 1$
.
For comparison, figure 7(b) also shows the average wavelength of the flow pattern in the DPS, SPS and NPS. Finer spatial resolution in the simulation, compared to that in experiment, facilitates measuring
$\bar{\unicode[STIX]{x1D706}}_{x}$
closer to onset with greater accuracy. In the DPS, the size of the domain along
$x$
was chosen a posteriori to be commensurate with the critical wavelength at onset in the experiment. Despite this, neither the spatial variation of the wavelength nor its variation with
$Re-Re_{c}$
observed in the experiment are captured. The SPS, however, shows a qualitatively similar trend for the dependence of
$\bar{\unicode[STIX]{x1D706}}_{x}$
on
$Re-Re_{c}$
. The periodicity in the transverse direction results in a uniform vortex pattern with smaller spatial variation in the separation between vortices compared to the experiment. In contrast, the NPS captures both the spatial variation of the wavelength and the distortion of the lattice with increasing forcing quite satisfactorily.
At
$Re-Re_{c}\approx 1$
, the discrepancy is much smaller than the uncertainty bars, but for
$Re-Re_{c}\gtrsim 1.5$
, the NPS overestimates the wavelength compared to what is observed in the experiment. The largest discrepancy, which is 0.46 (a 10 % relative error), occurs around
$Re-Re_{c}=3.7$
. The difference in the flow patterns in the NPS and the experiment is due to the deviation of the latter from being perfectly Q2D. The analysis in appendix C shows that the wavelength of the pattern sensitively depends on relatively minor changes in the forcing profile, which is responsible for the observed discrepancy between the numerics and experiment.
While the NPS provides a reasonably accurate description of the transition from the straight to the modulated flow in the experiment, as we mentioned previously, it somewhat underestimates the critical Reynolds number. To resolve this discrepancy, we tested the sensitivity of the transition in the NPS to changes in the forcing profile as well as variations in parameters
$\unicode[STIX]{x1D6FD}$
,
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D6FC}$
. In particular, we found that
$Re_{c}$
is fairly insensitive to spatial variation in the strength of the magnets in the array. Consequently, we turned our attention to studying the sensitivity of
$Re_{c}$
to the values of the parameters
$\unicode[STIX]{x1D6FD}$
,
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D6FC}$
. In order to match
$Re_{c}$
in the NPS and experiment, we had either to decrease
$\unicode[STIX]{x1D6FD}$
by 6 %, increase
$\unicode[STIX]{x1D708}$
by 7 %, or increase
$\unicode[STIX]{x1D6FC}$
by 22 %. Figure 8(a,b) shows that the variation of parameters has a fairly weak effect on both the amplitude of the modulation and the wavelength of the pattern, which suggests that the disagreement between the simulation and experiment is primarily due to the deviation of the flow and/or forcing from quasi-two-dimensionality, which is discussed in appendix C.
Note that the Reynolds number
$Re=Uw/\unicode[STIX]{x1D708}$
is defined using the measured r.m.s. velocity
$U$
and parameter
$\unicode[STIX]{x1D708}$
, which cannot be measured, but has to be computed. While in the simulation the value of
$\unicode[STIX]{x1D708}$
is well defined (it is one of the parameters of the model), in the experiment it is not, so the corresponding
$Re$
depends on the choice of
$\unicode[STIX]{x1D708}$
. Hence, to enable a proper comparison of experiment with numerics, we defined
$Re$
in both cases using the analytically computed depth-averaged value
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
(Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014), regardless of the actual value of
$\unicode[STIX]{x1D708}$
used in the simulation. Matching
$Re$
using this convention is effectively equivalent to matching the r.m.s. velocity
$U$
.
Figure 8. The effect of the variation in model parameters. (a) A bifurcation diagram and (b) the average wavelength of the pattern in the modulated flow regime. The numerical results correspond to either a 7 % increase in
$\unicode[STIX]{x1D708}$
, a 22 % increase in
$\unicode[STIX]{x1D6FC}$
, or a 6 % decrease in
$\unicode[STIX]{x1D6FD}$
compared with the depth-averaged values for the straight flow (
$\unicode[STIX]{x1D6FC}=0.064~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
). The uncertainty bars in (b) are only shown for every other data point for clarity. In all cases,
$Re$
is defined using
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
.
4.4 Secondary instability
As we increase the forcing further, the modulated flow in the experiment becomes unstable, giving way to a time-periodic flow with a period
$T_{p}=42.8\pm 0.4$
(
$120\pm 1$
s in dimensional units) at onset, which corresponds to
$Re_{p}=17.6\pm 0.1$
. The modulated state in the NPS, with the depth-averaged (
$\unicode[STIX]{x1D6FC}=0.064$
s
$^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
) as well as the adjusted parameters, undergoes a Hopf bifurcation as the forcing is increased. Table 1 compares
$Re_{p}$
and
$T_{p}$
in the experiment with those from the NPS for the different parameter sets. A video comparing the time-periodic flow in the NPS (with depth-averaged parameters) and experiment is included as online supplementary material available at https://doi.org/10.1017/jfm.2017.553.
Table 1. Critical transition parameters characterizing the stable periodic regime for the experiment and the NPS with different sets of parameters.
In the simulation
$Re_{p}$
and
$T_{p}$
are within
$15$
% of the experimental measurements for the depth-averaged values of parameters computed using the vertical profile
$P(z)$
that corresponds to the straight flow. However, the values of
$\unicode[STIX]{x1D708}$
,
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
should vary slowly with
$Re$
, since
$P(z)$
is weakly dependent on the horizontal flow profile. Hence, a different set of parameters is required to describe the two instabilities and, more generally, there is no universal set of parameters
$\unicode[STIX]{x1D6FD}$
,
$\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D6FC}$
that correctly describes the experimental flow at all
$Re$
. From table 1 we see that
$Re_{p}$
and
$T_{p}$
show very different sensitivities to changes in each of the parameters. Hence, while separately modifying
$\unicode[STIX]{x1D708}$
,
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
shows some improvement in matching either
$Re_{p}$
or
$T_{p}$
, it should be possible to obtain even better agreement by modifying all the model parameters simultaneously, each by only a few per cent.
The necessity for modifying parameters across different dynamical regimes also raises the question of how robust
$\unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D6FD}$
and
$\unicode[STIX]{x1D708}$
are to changes in the (local) wavenumber of the flow. To test this, we have recomputed the parameters using the wavenumber
$k\approx \sqrt{5/4}\,\unicode[STIX]{x1D705}$
associated with the modulated flow. We found that
$\unicode[STIX]{x1D6FD}$
and
$\unicode[STIX]{x1D708}$
change by less than 1 %, and
$\unicode[STIX]{x1D6FC}$
by approximately 3.5 %, compared to those computed using
$k=\unicode[STIX]{x1D705}$
. This robustness suggests that, once adjusted to match the experiment, the 2D model should provide a reasonably accurate description of the dynamics even in the weakly turbulent regime where the wavenumber may vary in space and time (Suri et al.
Reference Suri, Tithof, Grigoriev and Schatz2017).
5 Nature of the primary instability
An important consequence of confining the flow in the longitudinal or transverse directions is that we restrict the set of coordinate transformations (symmetries) that leave the governing equation (1.4) equivariant. The symmetries of the governing equation, in turn, determine the number of, and the relation between, distinct modulated flow solutions created as a result of the primary bifurcation. Below, we discuss each of the different flow domains, in the order of decreasing symmetry.
5.1 Doubly periodic simulation
On an unbounded or a doubly periodic domain, (1.4) is equivariant under the following symmetry operations (Chandler & Kerswell Reference Chandler and Kerswell2013):
-
(i) continuous shift by
$\unicode[STIX]{x1D6FF}x$
in
$x$
,
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}(x,y)\rightarrow (x+\unicode[STIX]{x1D6FF}x,y)$
; -
(ii) reflection in
$x$
combined with a discrete shift of half a period in
$y$
,
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}(x,y)\rightarrow (-x,y+w)$
; -
(iii) reflections in both
$x$
and
$y$
,
${\mathcal{R}}_{x}{\mathcal{R}}_{y}(x,y)\rightarrow (-x,-y)$
.
Note that the double reflection
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
is equivalent to a rotation by angle
$\unicode[STIX]{x03C0}$
about the
$z$
-axis, while the square of the symmetry operation
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
corresponds to a discrete shift
${\mathcal{T}}_{y}^{2w}$
in the
$y$
-direction, i.e.
$({\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w})^{2}={\mathcal{T}}_{y}^{2w}$
. For the DPS, the corresponding symmetry group is the semidirect product
$D_{n}\rtimes SO(2)$
(Armbruster et al.
Reference Armbruster, Nicolaenko, Smaoui and Chossat1996), where
$n$
is the number of magnets (here,
$n=8$
).
The above transformations that leave the governing equation equivariant, however, need not leave the flow fields invariant. When a flow field does not share a certain symmetry of the governing equation, one can generate – by applying the corresponding coordinate transformation – a dynamically equivalent symmetry-related copy of the flow. The straight flow in figure 3(a) remains unchanged when an arbitrary translation
$\unicode[STIX]{x1D6FF}x\in [0,L_{x}]$
is applied along the
$x$
-direction, so there is a unique solution
$\boldsymbol{u}_{s}$
. However, since the primary instability breaks the translational symmetry, there is a continuum of distinct modulated flow solutions
$\boldsymbol{u}_{m}$
related by translations in the
$x$
direction. This instability therefore corresponds to a circle pitchfork bifurcation (cf. figure 9
a).
Figure 9. A schematic showing the bifurcations corresponding to the primary instability: (a) circle pitchfork in the DPS, (b) sequence of pitchfork bifurcations in the SPS, (c) imperfect pitchfork bifurcation in the NPS. Solid (dashed) lines indicate stable (unstable) solution branches. The vertical and out-of-plane axes correspond to deviations of the flow from straight that are invariant under
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
and
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
, respectively.
The equivariance of the governing equation under
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}$
with arbitrary
$\unicode[STIX]{x1D6FF}x$
makes the choice of the coordinate origin
$x=0$
for a modulated flow arbitrary. We fix it by requiring that
$\boldsymbol{u}_{m}^{1}={\mathcal{R}}_{x}{\mathcal{R}}_{y}\boldsymbol{u}_{m}^{1}$
for a particular modulated flow solution
$\boldsymbol{u}_{m}^{1}$
. Since both the discrete symmetries
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
and
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
include reflection of the flow about the line
$x=0$
, the choice of the origin determines whether a particular solution remains invariant under either of these discrete symmetries. Figure 10 shows the four distinct solutions related by discrete translations
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}$
with
$\unicode[STIX]{x1D6FF}x=L_{x}/8$
:
$$\begin{eqnarray}\boldsymbol{u}_{m}^{3}={\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}\boldsymbol{u}_{m}^{1},\quad \boldsymbol{u}_{m}^{2}={\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}\boldsymbol{u}_{m}^{3},\quad \boldsymbol{u}_{m}^{4}={\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}\boldsymbol{u}_{m}^{2},\quad \boldsymbol{u}_{m}^{1}={\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}\boldsymbol{u}_{m}^{4}.\end{eqnarray}$$
Each of these four solutions is invariant under
${\mathcal{T}}_{y}^{2w}$
and either
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
or
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
. In particular,
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
are invariant under
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
, while
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
are invariant under
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
. Furthermore, the states
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
are related to each other via
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
, i.e.
$\boldsymbol{u}_{m}^{1}={\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}\boldsymbol{u}_{m}^{2}$
. Similarly,
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
are related via
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
, i.e.
$\boldsymbol{u}_{m}^{3}={\mathcal{R}}_{x}{\mathcal{R}}_{y}\boldsymbol{u}_{m}^{4}$
. Note that the operator
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
contains a single reflection, which causes the sign of the vorticity to change. In summary, by virtue of the continuous translational symmetry of the governing equation, the laminar flow in the DPS undergoes a circle pitchfork bifurcation with an infinite number of translation-related copies of a modulated flow. Only four of these copies, however, remain invariant under the discrete symmetries involving the reflection
${\mathcal{R}}_{x}$
.
Figure 10. Modulated flow fields (a)
$\boldsymbol{u}_{m}^{1}$
, (b)
$\boldsymbol{u}_{m}^{3}$
, (c)
$\boldsymbol{u}_{m}^{2}$
and (d)
$\boldsymbol{u}_{m}^{4}$
at
$Re=14$
in the DPS. The vorticity colour scale is the same as that in figure 6.
5.2 Singly periodic simulation
The no-slip boundary condition at
$x=\pm L_{x}/2$
in the SPS destroys the equivariance under translation
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}$
, reducing the symmetry group to
$D_{n}$
. The governing equation, however, still remains equivariant under each of the discrete transformations
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
and
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
. The loss of equivariance under
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}$
, which connected the states
$\boldsymbol{u}_{m}^{1}$
,
$\boldsymbol{u}_{m}^{2}$
with
$\boldsymbol{u}_{m}^{3}$
,
$\boldsymbol{u}_{m}^{4}$
in the DPS, implies either of
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
or
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
is broken in the straight to modulated transition in the SPS. Breaking either of the discrete symmetries,
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
or
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
, should generate (only) two branches, i.e. should result in a pitchfork bifurcation, since the modulated flow states in the SPS should remain symmetric with respect to
${\mathcal{T}}_{y}^{2w}$
, which is not affected by confinement in
$x$
. Consequently, of the infinite number of modulated states in DPS, only four, the counterparts of those shown in figure 10, will survive in the SPS, and should be formed via two distinct pitchforks.
This is indeed what we observe in the simulations, wherein two pairs of distinct solutions, shown in figure 11, are formed via two distinct pitchfork bifurcations of the straight flow. Just as in the DPS,
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
are invariant under
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
, while
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
are invariant under
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
. In figure 9(b) the
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
branches are plotted to lie in a plane perpendicular to that containing
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
. Since the bifurcations break either
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
or
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
, each pair of branches is related via the broken symmetry, i.e.
$\boldsymbol{u}_{m}^{1}={\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}\boldsymbol{u}_{m}^{2}$
and
$\boldsymbol{u}_{m}^{3}={\mathcal{R}}_{x}{\mathcal{R}}_{y}\boldsymbol{u}_{m}^{4}$
.
Unlike the DPS, where
${\mathcal{T}}_{x}^{\unicode[STIX]{x1D6FF}x}$
relates all the distinct solutions corresponding to the modulated flow (cf. (5.1)), there is no coordinate transformation that maps
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
to
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
. On an infinite domain, all four branches of the modulated flow are created at exactly the same
$Re$
(as in the DPS); however, on a finite domain, the pitchfork bifurcations that produce the two pairs of solutions would generally happen at different
$Re$
(cf. figure 9
b) that depend on the confinement in the
$x$
-direction, i.e. on
$L_{x}$
. For
$L_{x}=14$
, chosen from experimental considerations, the bifurcation that gives rise to
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
occurs at a higher
$Re$
than the bifurcation that gives rise to
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
. For other choices of
$L_{x}$
, the sequence may reverse. Note that the two modulated flow branches (
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
) that are formed from the second pitchfork are initially unstable, because they bifurcate off the unstable straight flow solution.
Figure 11. Modulated flow fields (a)
$\boldsymbol{u}_{m}^{1}$
, (b)
$\boldsymbol{u}_{m}^{3}$
, (c)
$\boldsymbol{u}_{m}^{2}$
and (d)
$\boldsymbol{u}_{m}^{4}$
at
$Re=14$
in the SPS. Vertical black lines indicate the central region, which is analogous to the flow fields shown in figure 10. The vorticity colour scale is the same as that of figure 6.
5.3 Non-periodic simulation
In the NPS, the additional no-slip boundary condition at
$y=\pm L_{y}/2$
breaks the equivariance of the problem under
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
, leaving the governing equation equivariant only under
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
. The pitchfork bifurcation that gives rise to the rotationally invariant solutions
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
in the SPS is associated with breaking of the
${\mathcal{R}}_{x}T_{y}^{w}$
symmetry. This symmetry is only approximate in the NPS for all
$Re$
, so one finds an imperfect pitchfork bifurcation instead, as shown in figure 9(c). The straight flow
$\boldsymbol{u}_{s}^{1}$
at lower
$Re$
smoothly transitions to the modulated flow
$\boldsymbol{u}_{m}^{1}$
at higher
$Re$
without an instability taking place, i.e. the real part of the leading eigenvalue of the straight flow does not change sign as we increase
$Re$
in the NPS. The shapes of the bifurcation curves close to
$Re_{c}$
(figure 7
a) showcase the difference in the nature of the primary instability between the NPS and the two periodic simulations.
In the NPS, the
$\boldsymbol{u}_{m}^{2}$
branch and the higher-
$Re$
branch of the straight flow
$\boldsymbol{u}_{s}^{2}$
are created in a saddle-node bifurcation at
$Re=10.72$
. The states
$\boldsymbol{u}_{m}^{1}$
and
$\boldsymbol{u}_{m}^{2}$
, both of which are symmetric with respect to
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
, are shown in figure 12. While
${\mathcal{R}}_{x}T_{y}^{w}$
is not an exact symmetry in the NPS, given the large transverse extent of the domain compared with the period of the forcing (
$L_{y}/2w$
= 9), near the centre of the domain this approximate symmetry holds and consequently
$\boldsymbol{u}_{m}^{2}\approx {\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}\boldsymbol{u}_{m}^{1}$
. However, unlike
$\boldsymbol{u}_{m}^{1}$
, which remains stable up to
$Re=15.4$
in the NPS,
$\boldsymbol{u}_{m}^{2}$
is unstable over the entire range of
$Re$
where it exists (
$Re\geqslant 10.72$
). This explains why our numerical simulations starting from randomized initial conditions have always converged to the modulated flow
$\boldsymbol{u}_{m}^{1}$
and why
$\boldsymbol{u}_{m}^{2}$
was never found in the numerical simulations or observed in the experiment.
Figure 12. Modulated flow fields at
$Re=11.6$
in the NPS beyond the imperfect pitchfork bifurcation shown in figure 9(c). The flow fields shown here are (a)
$\boldsymbol{u}_{m}^{1}$
, which emerges smoothly from the straight flow, and (b)
$\boldsymbol{u}_{m}^{2}$
, which is formed through a saddle-node bifurcation.
The pitchfork bifurcation, which gave rise to the branches
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
in the SPS, also does not carry over into the NPS. Instead,
$\boldsymbol{u}_{s}^{2}$
undergoes a Hopf bifurcation at
$Re=12.6$
. This change in the nature of the bifurcation is probably caused by transverse confinement, which has a more prominent effect on
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
: these states are invariant under the
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
symmetry in the SPS, but this symmetry is broken in the NPS. While we fail to observe
$\boldsymbol{u}_{m}^{3}$
and
$\boldsymbol{u}_{m}^{4}$
in the NPS, the analogues of these solutions may appear for a different set of model parameters, forcing profile and/or degree of confinement. Details regarding the computation of the unstable branches associated with the various bifurcations are included in § B.3.
6 Conclusions
In this article, we have presented a combined experimental and numerical study of bifurcations in a Q2D Kolmogorov-like flow. This flow is realized in the laboratory by electromagnetically driving a stratified layer of electrolyte above an immiscible layer of dielectric. This Q2D flow is described using a 2D model (1.4) derived from first principles by depth-averaging the 3D Navier–Stokes equation. In contrast, virtually all previous studies have modelled Q2D flows using (1.2), a semi-empirical variation of the 2D Navier–Stokes equation with the addition of a linear friction. Also unlike previous studies of Kolmogorov-like flows, which have assumed a perfectly sinusoidal forcing profile, we have introduced a realistic model of the forcing, which has been validated against 3D experimental measurements.
To test the importance of lateral confinement, we have compared experimental measurements with numerical simulations using different boundary conditions. We have found that, by incorporating realistic, no-slip boundary conditions at all lateral boundaries and a realistic forcing profile, quantitative agreement between the experiment and simulation can be achieved with no adjustable parameters. In particular, the Reynolds number
$Re_{c}$
for the primary instability can be predicted to within approximately 5 % and the critical wavenumber
$k_{c}$
can be predicted to an accuracy higher than the measurement accuracy. These are significant improvements compared with previous studies, none of which were able to predict both
$Re_{c}$
and
$k_{c}$
with this level of accuracy, despite using adjustable parameters.
We have also performed a systematic study of how lateral confinement affects the nature of the bifurcation describing the transition from the straight flow to the modulated flow. Previous studies have characterized this transition in the experiment as a pitchfork bifurcation, using analytical computations on a periodic domain. We have shown that, because of confinement, an imperfect pitchfork bifurcation is found instead. We have also numerically computed the two unstable branches of the imperfect pitchfork bifurcation describing flows that are not observed in either experiment or simulations under normal conditions.
Furthermore, we have demonstrated that the model reasonably accurately predicts the modulated flow pattern (both the wavenumber and the amplitude) beyond the onset of the primary instability. Moreover, this is the first study, experimental or theoretical, to provide a quantitative analysis of the secondary instability of a Kolmogorov-like flow that generates a time-dependent pattern of vortices. Even for the secondary instability, the numerical predictions of the critical Reynolds number
$Re_{p}$
and the critical period
$T_{p}$
are in general agreement with the experiment, although the accuracy of the numerical predictions decreases with increasing
$Re$
.
The discrepancy between the numerical predictions and experiments has been traced back to the variation of the forcing profile with height. This points to the limitations of a 2D model of what in reality is a 3D flow, albeit with a strongly suppressed vertical component of the velocity. Nonetheless, for a select range of
$Re$
, the experimental flow can be reproduced with extremely good quantitative accuracy by making fairly small adjustments to the model parameters, compared with their depth-averaged values computed for the simple straight flow. The ability of the model to closely reproduce the experimental flow is crucial for the utility of Q2D flows for testing the geometrical description of weakly turbulent flows and studying the dynamical role of exact coherent structures (Suri et al.
Reference Suri, Tithof, Grigoriev and Schatz2017). Such tests will be the main focus of follow-up studies.
Acknowledgements
We thank R. Kerswell for useful discussions and D. Borrero for his useful suggestions and insight. J.T. is grateful to S. Raben for his help with the PRANA PIV software package. This work was supported in part by the National Science Foundation under grants no. CMMI-1234436 and DMS-1125302.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2017.553.
Appendix A. Scaling and non-dimensionalization
The governing equation (1.4) was presented in dimensional form to highlight the dependence of parameters
$\unicode[STIX]{x1D6FC}$
,
$\unicode[STIX]{x1D6FD}$
and
$\unicode[STIX]{x1D708}$
on the properties of the fluid layers. The dimensional form also makes it easier to explore the sensitivity of the dynamics to changes in these parameters. To simplify comparison of our results with other studies, it is helpful to non-dimensionalize this equation. Choosing the width of a magnet
$w$
as the length scale, the r.m.s. velocity computed over the central region
$|x|\leqslant 4w$
,
$|y|\leqslant 4w$
as the velocity scale
$U$
, and the ratio of these two scales as the time scale, one obtains the following non-dimensional equation:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D6FD}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p_{0}+\frac{1}{Re}(\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D6FE}\boldsymbol{u})+\boldsymbol{f}_{0},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FC}w^{2}/\unicode[STIX]{x1D708}$
describes the relative strength of the Rayleigh friction and viscous terms in (1.4). Finally,
$p_{0}=p/(U^{2}\bar{\unicode[STIX]{x1D70C}})$
is the non-dimensional pressure and
$\boldsymbol{f}_{0}=w\boldsymbol{f}/U^{2}$
is the non-dimensional forcing profile.
It is possible to eliminate the parameter
$\unicode[STIX]{x1D6FD}$
from (A 1) by making the length, time and velocity scales independent. If we again choose the width of a magnet
$w$
as the length scale, the r.m.s. velocity as the velocity scale
$U$
, and
$w/(\unicode[STIX]{x1D6FD}U)$
as the time scale, we instead obtain the following non-dimensional equation:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p_{1}+\frac{1}{Re^{\prime }}(\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\unicode[STIX]{x1D6FE}\boldsymbol{u})+\boldsymbol{f}_{1},\end{eqnarray}$$
where
$Re^{\prime }=Re/\unicode[STIX]{x1D6FD}$
,
$p_{1}=\unicode[STIX]{x1D6FD}p/(U^{2}\bar{\unicode[STIX]{x1D70C}})$
is the non-dimensional pressure, and
$\boldsymbol{f}_{1}=w\unicode[STIX]{x1D6FD}\boldsymbol{f}/U^{2}$
is the non-dimensional forcing profile. Although (A 2) does not contain
$\unicode[STIX]{x1D6FD}$
explicitly, the Reynolds number is rescaled by
$\unicode[STIX]{x1D6FD}$
. Hence, the non-dimensional equations (A 1) and (A 2) as well as the dimensional equations (1.2) and (1.4) predict an identical sequence of bifurcations. However, the critical values of
$Re$
scale as
$1/\unicode[STIX]{x1D6FD}$
, as we have found explicitly in (4.1).
Appendix B. Numerical methods
In this appendix, we present the details of discretization methods and numerical integration schemes employed in the NPS, SPS and DPS. Additionally, we also detail the computation of the unstable branches associated with the pitchfork bifurcation that cannot be obtained from simple numerical integration.
B.1 Non-periodic and singly periodic simulations
Since the NPS, as well as the SPS, requires prescribing no-slip (e.g. Dirichlet) boundary conditions on the velocity field
$\boldsymbol{u}$
, numerical simulations are performed using the primitive variable (
$u_{x}$
,
$u_{y}$
and
$p$
) formulation by employing a semi-implicit fractional-step method detailed in Armfield & Street (Reference Armfield and Street1999). Temporal discretization of (1.4) is performed using the following difference scheme:
$$\begin{eqnarray}\frac{\boldsymbol{u}_{n+1}-\boldsymbol{u}_{n}}{\unicode[STIX]{x0394}t}+\frac{3}{2}{\mathcal{N}}\boldsymbol{u}_{n}-\frac{1}{2}{\mathcal{N}}\boldsymbol{u}_{n-1}=-\frac{1}{\bar{\unicode[STIX]{x1D70C}}}\unicode[STIX]{x1D735}p_{n+1}+\frac{1}{2}{\mathcal{L}}(\boldsymbol{u}_{n+1}+\boldsymbol{u}_{n})+\boldsymbol{f}.\end{eqnarray}$$
In the above equation,
$\boldsymbol{u}_{n}$
and
$p_{n+1}$
are the velocity and pressure fields, with the subscript
$n$
indicating a discrete time instant
$t_{n}=n\unicode[STIX]{x0394}t$
, where
$\unicode[STIX]{x0394}t$
is the time step for the update. For purposes of brevity, we have used the notation
${\mathcal{N}}\boldsymbol{u}_{n}=\unicode[STIX]{x1D6FD}\boldsymbol{u}_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{n}$
and
${\mathcal{L}}\boldsymbol{u}_{n}=\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{n}-\unicode[STIX]{x1D6FC}\boldsymbol{u}_{n}$
to represent the nonlinear and linear terms, respectively. The above discretization is a semi-implicit approximation of (1.4), where the linear terms in the update are treated implicitly using the Crank–Nicolson scheme, while the nonlinear term is handled explicitly using the Adams–Bashforth scheme. The velocity field
$\boldsymbol{u}_{n+1}$
at every instant satisfies the incompressibility condition:
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{n+1}=0,\end{eqnarray}$$
which is enforced on each update through the three-fractional-step P2 (pressure correction) projection method discussed in Armfield & Street (Reference Armfield and Street1999).
Spatial discretization of the velocity and pressure fields is carried out using the standard marker and cell (MAC) staggered grid (Harlow & Welch Reference Harlow and Welch1965). The spatial derivatives in (B 1) are approximated using finite central differences; the 2D Laplacian operator (
$\unicode[STIX]{x1D6FB}^{2}$
) uses a five-point stencil formula; and the nonlinear term uses the three-point central difference formula.
For both the NPS and the SPS, we have chosen 20 cells per magnet width
$w$
to discretize the velocity and pressure fields. Since the dimensions of the NPS are identical to the lateral dimensions of the experiment, i.e.
$14w\times 18w$
, a total of
$280\times 360$
cells were used to sample the flow domain. The SPS, however, corresponds to a domain of dimensions
$14w\times 8w$
, which maps to a region including the central eight magnets in the experiment. Hence, a total of
$280\times 160$
cells were used to discretize the SPS domain. For both the SPS and NPS, a time step of
$\unicode[STIX]{x0394}t=1/40$
s was used for all the numerical simulations.
To test the adequacy of the spatial resolution, the velocity field corresponding to the modulated flow at
$Re\approx 15.5$
was recomputed by doubling the resolution, i.e. using 40 cells per magnet width. To compare this velocity field (
$\boldsymbol{u}_{40}$
) with the one computed on the 20-cell grid (
$\boldsymbol{u}_{20}$
), we interpolated
$\boldsymbol{u}_{40}$
onto the 20-cell grid to obtain
$\boldsymbol{u}_{interp}$
(interpolation was required due to the staggered nature of the grid). The difference between
$\boldsymbol{u}_{interp}$
and
$\boldsymbol{u}_{20}$
, computed as
$\Vert (\boldsymbol{u}_{interp}-\boldsymbol{u}_{20})\Vert /\Vert \boldsymbol{u}_{20}\Vert$
, was 1.2 %. Since the interpolation introduces error, it can be concluded that the actual error should be less than 1.2 %. We find that global measures, such as
$Re$
or
$\langle u_{y}^{2}\rangle$
(used to characterize the primary instability), computed directly using
$\boldsymbol{u}_{40}$
and
$\boldsymbol{u}_{20}$
differed by less than 0.2 %. These tests confirm that a resolution of 20 cells per magnet width is sufficient to simulate the flow and characterize the bifurcations accurately.
B.2 Doubly periodic simulation
Simulations on the doubly periodic domain can be sped up significantly using a spectral method (Canuto et al.
Reference Canuto, Hussaini, Quarteroni and Zang1988). Since solving linear equations involving the Laplacian is very cheap in the spectral method, it is convenient to use the vorticity–streamfunction formulation instead of the velocity–pressure formulation. Taking the curl of (1.4), we obtain the following equation for the
$z$
-component of vorticity
$\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D735}\times \boldsymbol{u})\boldsymbol{\cdot }\hat{\boldsymbol{z}}$
:
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FD}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D714}+W,\end{eqnarray}$$
where
$W=(\unicode[STIX]{x1D735}\times \boldsymbol{f})\boldsymbol{\cdot }\hat{\boldsymbol{z}}$
. The horizontal components of the velocity field
$u_{x}=\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}y$
and
$u_{y}=-\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}x$
can be computed using the streamfunction
$\unicode[STIX]{x1D713}$
, which satisfies the Poisson equation
$\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D714}$
.
The vorticity field
$\unicode[STIX]{x1D714}$
is discretized in the Fourier space using
$128$
modes along each of the
$x$
- and
$y$
-directions. Since the lateral dimensions of the periodic domain are
$8w\times 8w$
units, the spatial resolution associated with the Fourier grid corresponds to 16 grid points per magnet width
$w$
. Taking the Fourier transform of (B 3), we obtain
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6FA}=-\unicode[STIX]{x1D6FD}{\mathcal{F}}[\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}]+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6FA}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FA}+{\mathcal{F}}[W],\end{eqnarray}$$
where
${\mathcal{F}}[\cdot ]$
represents the Fourier transform and
$\unicode[STIX]{x1D6FA}={\mathcal{F}}[\unicode[STIX]{x1D714}]$
.
Equation (B 4) is stepped forward in time (
$t\rightarrow t+\unicode[STIX]{x0394}t$
) using a three-substep semi-implicit Strang–Marchuk splitting algorithm (Ascher, Ruuth & Wetton Reference Ascher, Ruuth and Wetton1995; Mitchell Reference Mitchell2013), where the first and last substeps advance the vorticity field using the nonlinear term by means of a second-order explicit Runge–Kutta scheme (using a time step
$\unicode[STIX]{x0394}t/2$
), while the intermediate substep advances the vorticity field using the Crank–Nicolson scheme (using a time step
$\unicode[STIX]{x0394}t$
). We have used the time step
$\unicode[STIX]{x0394}t=1/32~\text{s}$
.
B.3 Computing unstable branches in the pitchfork bifurcation
The schematic depicting the pitchfork bifurcation in figure 9 was constructed following the computation of all the stable and unstable states using the matrix-free Newton–Krylov solver (Kelley Reference Kelley2003). Guesses for the stable states, to initialize the Newton solver, can be easily obtained using numerical integration. However, those for the unstable states should be constructed using continuation or using the eigenmode that goes unstable at the bifurcation.
To begin with, initial guesses for the unstable straight flow branches in the SPS and NPS simulations were constructed by extrapolating the stable straight solutions in Reynolds number
$Re$
, i.e.
$$\begin{eqnarray}\boldsymbol{u}_{s}^{2}(Re_{c}+\unicode[STIX]{x1D716})\approx \boldsymbol{u}_{s}^{1}(Re_{c})+\unicode[STIX]{x1D716}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{s}^{1}}{\unicode[STIX]{x2202}Re}\right)_{Re_{c}},\end{eqnarray}$$
where the derivative
$(\unicode[STIX]{x2202}\boldsymbol{u}_{s}^{1}/\unicode[STIX]{x2202}Re)$
was approximated using finite differences,
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{s}^{1}}{\unicode[STIX]{x2202}Re}\right)_{Re}\approx \frac{\boldsymbol{u}_{s}^{1}(Re)-\boldsymbol{u}_{s}^{1}(Re-\unicode[STIX]{x0394}Re)}{\unicode[STIX]{x0394}Re}.\end{eqnarray}$$
This method proved particularly useful in obtaining a good initial guess for the unstable straight flow
$\boldsymbol{u}_{s}^{2}$
in the NPS, since it is disconnected from
$\boldsymbol{u}_{s}^{1}$
, as shown in figure 9(c). In the NPS,
$Re_{c}\approx 10.5$
is not estimated by identifying the instability of
$\boldsymbol{u}_{s}^{1}$
, since there exists none. Instead, it is computed using the intercept of a linear fit of the amplitude
$\langle u_{y}^{2}\rangle$
versus
$Re$
, shown in figure 7, close to the onset of modulation. The initial guess at
$Re\approx 10.75$
was constructed by extrapolating
$\boldsymbol{u}_{s}^{1}$
from
$Re\approx 10.25$
by choosing
$\unicode[STIX]{x1D716}=0.5$
in (B 5).
For the unstable modulated branches emerging from the second pitchfork in the SPS, a good initial guess for
$\boldsymbol{u}_{m}^{3}$
(
$\boldsymbol{u}_{m}^{4}$
) is constructed using
$\boldsymbol{u}_{m}^{3}\approx \boldsymbol{u}_{s}^{2}\pm p\,\hat{\boldsymbol{e}}_{2}$
. Here
$\hat{\boldsymbol{e}}_{2}$
is the second unstable eigenvector of the straight flow which has the symmetry
${\mathcal{R}}_{x}{\mathcal{T}}_{y}^{w}$
. Since amplitude
$p$
is not known a priori, convergence to
$\boldsymbol{u}_{m}^{3}$
is tested by incrementing
$p$
. In the NPS, the initial guess for the unstable branch
$\boldsymbol{u}_{m}^{2}$
was similarly constructed,
$\boldsymbol{u}_{m}^{2}\approx \boldsymbol{u}_{s}^{2}-p\,\hat{\boldsymbol{e}}_{1}$
, using the eigenvector
$\hat{\boldsymbol{e}}_{1}$
with the
${\mathcal{R}}_{x}{\mathcal{R}}_{y}$
symmetry. However,
$p$
in the NPS case can be estimated using
$p=\langle \hat{\boldsymbol{e}}_{1}\mid \boldsymbol{u}_{m}^{1}-\boldsymbol{u}_{s}^{2}\rangle$
, since the stable modulated flow
$\boldsymbol{u}_{m}^{1}$
is known from numerical integration and
$\boldsymbol{u}_{s}^{2}$
is computed from extrapolation.
Appendix C. Inherent three-dimensionality of the forcing in the experiment
In § 4.3 we compared measurements of the average longitudinal wavelength
$\bar{\unicode[STIX]{x1D706}}_{x}$
of the modulated flow in the experiment and numerical simulations. The comparison between the experiment and the NPS (cf. figures 7 and 8) showed systematic differences in
$\bar{\unicode[STIX]{x1D706}}_{x}$
for both the depth-averaged and modified parameters, with the maximum difference being approximately 10 % at
$Re-Re_{c}=3.7$
(cf. figure 7). Here we show that this deviation is probably due to the inherent three-dimensionality of the experiment, not captured by a strictly 2D model, rather than the choice of the model parameters.
As we discussed in § 3.1, the Lorentz force density due to the specific arrangement of magnets employed in the experiment is to a very good approximation given by
$\boldsymbol{F}=JB_{z}\hat{\boldsymbol{x}}$
, where
$J$
is the magnitude of current density and
$B_{z}(x,y,z)$
is the
$z$
-component of the magnetic field at any given location within the electrolyte. In deriving (1.4) it was assumed (Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014) that
$B_{z}$
can be decomposed as the product of a 2D horizontal profile
$B_{2D}(x,y)$
, which depends exclusively on the extended coordinates
$(x,y)$
, and a 1D vertical profile
$D(z)$
, which captures the variation of the magnetic field above the magnet array, i.e.
$B_{z}(x,y,z)=D(z)B_{2D}(x,y)$
. This implies that, when normalized, the planar magnetic field profiles at various heights
$z$
within the electrolyte are identical. Such a magnetic field, which we call ‘Q2D’, facilitates the decomposition of the plane-parallel Q2D velocity field (1.3) which underpins the strictly 2D model (1.4).
A magnetic field that is truly Q2D, however, cannot be created using a magnet array with finite dimensions, i.e. the shape of the magnetic field profile generated by permanent magnets in the laboratory always changes with the vertical height
$z$
to some extent. Experimental measurements of the magnetic field from previous studies have shown such changes in the shape of the field profile as a function of
$z$
(Dovzhenko et al.
Reference Dovzhenko, Krymov and Ponomarev1984; Suri et al.
Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014). This is very much the case in our experiment as well, as can be seen from the magnetic field profiles shown in figure 2(a). For instance, if one rescales the transverse magnetic field profiles at heights
$z=0.438$
and
$z=0.265$
such that they match near the centre of the array, we see that these profiles would not match near the end magnets. This is most apparent by comparing the relative heights of the peaks at
$y=-6.5$
and
$y=-4.5$
for the two profiles in figure 2(a). A quantitative estimate of the deviation from quasi-two-dimensionality can be obtained by comparing the magnetic field profiles computed using the dipole summation at the bottom
$B_{b}=B_{z}(x,y,z=0.236)$
and the top of the electrolyte layer
$B_{t}=B_{z}(x,y,z=0.472)$
. Normalizing
$B_{b}$
and
$B_{t}$
separately, using their respective spatial r.m.s. values computed over the entire lateral extent of the domain, we estimate the largest difference between the profiles to be approximately 12 %. This difference is fairly localized towards the ends of the magnet array and is likely to be the reason behind the larger discrepancy in longitudinal velocity measurements over the end magnets (cf. figure 4
b).
To demonstrate the impact of the
$z$
dependence of the forcing profile on the flow, we have recomputed the straight and modulated flow fields in the NPS using
$B_{2D}=B_{b}$
and
$B_{2D}=B_{t}$
, in addition to the depth-averaged magnetic field
$B_{da}$
. The value of
$\unicode[STIX]{x1D6FC}$
was increased by 22 % relative to the depth-averaged value to reduce the influence of the uncertainty in the model parameters on the flow pattern. This choice also yields the best agreement between the average wavelengths of the flow pattern in the simulation and experiment for
$Re_{c}<Re<Re_{p}$
(see figure 8
b). As figure 13(a) shows, the forcing profile strongly affects both
$Re_{c}$
and the amplitude of the modulation of the flow for
$Re>Re_{c}$
. It also shows that
$B_{da}$
produces substantially better agreement with experiment than either
$B_{b}$
or
$B_{t}$
. Similarly, we find that the forcing profile strongly influences the modulation wavelength. As figure 13(b) shows, for both
$B_{b}$
or
$B_{da}$
the average wavelength agrees reasonable well with experiment, while
$B_{t}$
produces a very poor agreement. Similar results (not shown) are obtained if, instead of
$\unicode[STIX]{x1D6FC}$
, either
$\unicode[STIX]{x1D6FD}$
or
$\unicode[STIX]{x1D708}$
is modified to match
$Re_{c}$
.
The above analysis shows that, although the NPS with the depth-averaged magnetic field profile captures the salient features of the dynamics fairly well, the flow pattern depends fairly sensitively on the details of the forcing. Hence, one should expect systematic deviations between the 2D model derived for a Q2D flow and the experiment where quasi-two-dimensionality is broken by the forcing. It should be mentioned that the wavelength of the modulated flow measured by either seeding the dielectric–electrolyte interface or the top surface of the electrolyte are virtually identical. This implies that viscous coupling across the fluid layers produces a Q2D flow despite the fact that the forcing profile driving the flow is not perfectly Q2D.
Figure 13. Sensitivity to the magnetic field profile for
$\unicode[STIX]{x1D6FC}=0.078~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FD}=0.83$
and
$\unicode[STIX]{x1D708}=3.26\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$
. (a) A bifurcation diagram for the primary instability and (b) the average wavelength of the pattern in the modulated regime. The simulations were performed with
$\unicode[STIX]{x1D6FC}$
increased by 22 % relative to the depth-averaged value for the straight flow and used either
$B_{da}$
,
$B_{b}$
, or
$B_{t}$
.
