## Abstract

We investigate electrostatically induced interfacial instabilities and subsequent generation of nonlinear coherent structures in immiscible, viscous, dielectric multi-layer stratified flows confined in small-scale channels. Vertical electric fields are imposed across the channel to produce interfacial instabilities that would normally be absent in such flows. In situations when the imposed vertical fields are constant, interfacial instabilities emerge due to the presence of electrostatic forces, and we follow the nonlinear dynamics via direct numerical simulations. We also propose and illustrate a novel pumping mechanism in microfluidic devices that does not use moving parts. This is achieved by first inducing interfacial instabilities using constant background electric fields to obtain fully nonlinear deformations. The second step involves the manipulation of the imposed voltage on the lower electrode (channel wall) to produce a spatio-temporally varying voltage there, in the form of a travelling wave with pre-determined properties. Such travelling wave dielectrophoresis methods are shown to generate intricate fluid–surface–structure interactions that can be of practical value since they produce net mass flux along the channel and thus are candidates for microfluidic pumps without moving parts. We show via extensive direct numerical simulations that this pumping phenomenon is a result of an externally induced nonlinear travelling wave that forms at the fluid–fluid interface and study the characteristics of the generated velocity field inside the channel.

## 1. Introduction

The control of fluid interfaces in multi-layer flows has been the subject of numerous scientific and industrial studies in recent years. Applications in microfluidics, polymer self-assembly and lab-on-a-chip devices are only a few of the topics that have received much attention both theoretically and experimentally (see [1–3] and references therein).

Fluid dynamics on the microscale poses numerous engineering challenges. On the one hand, the absence or near absence of inertia can render such flows stable thus preventing the efficient utilization of physical processes, such as mass and heat transfer. On the other hand, the presence of multiple interfaces in immiscible viscous multi-layer flows can produce instabilities even at zero Reynolds numbers [4], and this may be undesirable in applications such as cell-sorting [5].

It is well known from the seminal work of Melcher & Taylor [6–8] that a vertical electric field acting on an otherwise quiescent viscous liquid layer induces long-wave instabilities with a band of unstable waves (and corresponding maximum growth rate) that increases as the electric field strength is increased. This phenomenon has led to several technological soft-lithography applications that are based on the appearance and control of morphological instabilities (see [9] and references therein). When the field acts horizontally, however, it can stabilize phenomena such as the Rayleigh–Taylor instability when a heavy fluid lies above lighter fluid. Linear theory and direct numerical computations based on the Navier–Stokes equations [10] have illustrated the stabilization and have also provided a method of generating and sustaining nonlinear time-periodic oscillations of the unstable interface by switching the field on and off sequentially. Additional work remains to be done especially in identifying electric field modulations that can produce optimal mixing, for example. For completeness, we mention analogous linear and nonlinear stability results in electrified shear flows, where two-layer Couette or falling film flows can be destabilized by vertical electric fields [11–13].

An important problem on the microscale is that of mixing, and typically this is achieved by applying time-pulsing in T-junctions as in the T-mixer in [14]. In addition, it has been shown experimentally that the instabilities due to vertical electric fields in two-layer shear flows in channels can be used to generate a monodisperse distribution of droplets of one fluid suspended in a second fluid [15]. Simulations that follow the instabilities into the nonlinear regime and interfacial touchdown (the crucial mechanism that generated the encapsulated droplets) have been carried out and the reduction in droplet volume as the field is increased is supported by them [16]. Allied direct numerical simulations on the effect of electric fields on both buoyantly driven or shear drops can be found in [17,18].

In this work, we propose a novel approach for the generation of nonlinear travelling waves inside a channel, which does not require the manipulation of an imposed velocity field or any moving parts. The electric field dynamics is the sole driving force in the system and we show how the motion of the interface can be accurately controlled in order to achieve directional velocity fluxes from a state which is originally at rest. Initially, a uniform vertical electric field is used (optimized using linear stability results) to steer the fluid–fluid interface towards the vicinity of the electrodes. This is followed by superimposing a travelling wave-like component for the electric field at one of the boundaries, which subsequently affects the nonlinear interfacial dynamics in such a way as to generate a controlled horizontal motion. A nonlinear travelling wave forms, which can be sustained as long as the travelling wave electric field boundary condition is kept switched on. This in turn affects the flow-field globally and indeed generates a net mass flux in one direction. A travelling wave electric field boundary condition has been used successfully in electrophoretically transporting particles inside microfluidic channels [19–21].

The underlying phenomena considered here are inherently nonlinear, and we undertake direct numerical simulations to describe them (linear theory is also included in order to guide and evaluate the accuracy of the computations). We use the GERRIS [22] volume-of-fluid (VOF) flow solver that has been tested extensively [23–25], with excellent results. The electrodynamics package is based on the work of López-Herrera *et al*. [25].

The remainder of this paper is organized as follows. Section 2 describes the mathematical model, highlights the geometry of the problem and presents the governing equations; it also presents linear stability results and the effect of the electric field on stability. Section 3 is dedicated to direct numerical simulations, starting with instabilities under uniform electric fields (along with comparisons with linear theory); it then describes the travelling wave electric field boundary condition used in the simulations; finally, it concludes with the presentation of extensive numerical results, including the generation of nonlinear interfacial travelling waves, a description of the flow field and the emergence of net horizontal flow, and finally a quantification of the net horizontal flux as electric field parameters vary. Section 4 is dedicated to our conclusion.

## 2. Mathematical model

We consider two incompressible, immiscible, viscous fluids in a two-dimensional setting as shown in the schematic of the domain presented in figure 1. A Cartesian coordinate system is used with channel walls parallel to the horizontal *x*-axis. The channel has height *L* and is taken to be sufficiently long—moderate to large aspect ratios are quite common in the experimental literature. In our analytical and computational work, periodic boundary conditions are imposed in the horizontal direction. By modelling the channel walls as electrodes and imposing constant voltages there, a potential difference is maintained laterally across the channel and a constant electric field of magnitude is generated vertically.

The interface is denoted by *y*=*S*(*x*,*t*) and the lower fluid 1 occupies the region −*L*/2<*y*<*S*(*x*,*t*), while the upper fluid 2 lies inside the domain *S*(*x*,*t*)<*y*<*L*/2. We use subscripts 1,2 to refer to quantities in each corresponding region. With this notation, the velocity fields and pressures are denoted by **u**_{1,2}=(*u*_{1,2},*v*_{1,2}) and *p*_{1,2}, while the relevant fluid properties are the constant densities *ρ*_{1,2}, constant viscosities *μ*_{1,2} and constant permittivities *ϵ*_{1,2} (the fluids are assumed to be perfect dielectrics). Surface tension is also present and has constant coefficient *σ*. To retain generality, we keep the effects of gravity via the acceleration defined by *g*.

Electric fields are governed by the electrostatic approximation of Maxwell's equations (this assumption is valid since induced magnetic fields in the problems studied here are negligible). In this limit, Maxwell's equations become ∇×**E**_{1,2}=0, ∇⋅(*ϵ*_{1,2}**E**_{1,2})=0 and the former equation implies that the electric fields can be written as **E**_{1,2}=−∇*V* _{1,2}. It follows from the second equation (Gauss's law) that away from the interface the voltage potentials satisfy Laplace's equation in each region of interest.

### (a) Governing equations

The dimensional momentum and continuity equations are
2.1
2.2
2.3
Note that Lorentz forces are absent in the momentum equations since the fluids have constant electrical properties and are dielectrics with no charges present. The coupling between hydrodynamics and electrostatics enters through the nonlinear interfacial boundary conditions discussed in the following paragraphs. As previously mentioned, the electric potential equations in each fluid are
2.4
The boundary conditions at the walls *y*=±*L*/2 are those of no slip for **u** and for the potentials we impose
2.5
where is a constant. At the sharp interface *y*=*S*(*x*,*t*), we must impose kinematic conditions, continuity of normal and tangential stress balances, continuity of velocities, continuity of the voltage displacement field (Gauss's law) and continuity of voltage potentials. These read
2.6
2.7
2.8
2.9
2.10
2.11
where represents the jump across the interface, , are the unit normal and tangent to the interface, respectively, and is the interfacial curvature. The stress tensor contains hydrodynamic and electrostatic contributions and is given by
2.12

and it is understood that variables on the right-hand side of (2.12) are evaluated in each fluid region *i*=1,2. Scaling lengths by *L*, velocities by *U* (this is a reference value, e.g. for capillary scales *U*=*σ*/*μ*_{1}), time by *L*/*U* and pressures by *ρ*_{1}*U*^{2} gives rise to the following set of dimensionless parameters:
2.13
where *ϵ*_{0} is the free-space permittivity. The parameters in (2.13) correspond to a dimensionless viscosity (inverse Reynolds number), an inverse Weber number, an inverse square Froude number, a voltage scaling, and density, viscosity and permittivity ratios, respectively. The apparently unusual scaling for voltages is motivated by the numerical work and is selected so that the dimensionless parameter in front of the Maxwell stresses in (2.12) is unity. Note also that with *U*=*σ*/*μ*_{1}, the parameter , where *O*_{h} is the Ohnesorge number. Using these scalings provides the dimensionless systems
2.14
2.15
2.16
where **j**=(0,1) and tildes are used to denote dimensionless quantities. The Laplace equations (2.4) are unchanged, and the boundary conditions at the channel walls become
2.17
and
2.18
where . Finally, the dimensionless form of the interfacial boundary conditions (2.6) maintain their form with all variables and parameters gaining tilde decorations in becoming dimensionless; the dimensionless stress tensor reads
2.19
For details regarding the electrohydrodynamic coupling through the stress tensor (2.19), the reader is referred to [26,27].

### (b) Linear stability

In this section, we carry out a linear analysis to determine growth rates in the presence of destabilizing electric fields. Apart from fundamental understanding of the physics of the problem, the results also serve as a basis for comparison with fully nonlinear computations presented later.

The base state solution is given by a flat interface, zero velocities and a uniform vertical electric field in each of the two fluids. There is a constant undisturbed induced electric pressure difference at the interface, which can be calculated from the normal stress balance to find
2.20
Hence we write
2.21
and
2.22
and linearize the equations and boundary conditions for *δ*≪1. Normal modes are introduced in the form , , e^{ikx+ωt}+c.c., , where c.c. denotes complex conjugates, and relevant eigenfunctions are determined analytically, using appropriate boundary conditions. For example, the perturbation voltage eigenfunctions can be expressed in terms of the interface perturbation eigenfunction *S* as follows:
2.23
The boundary conditions are manipulated to generate a system of homogeneous linear equations for the unknown constants of integration. For prescribed densities, viscosities, permittivities, surface tension and electric potential difference, we obtain a transcendental equation for the dispersion relation *ω*(*k*), which is solved numerically to determine the most unstable waves for a given wavenumber *k* (a temporal stability analysis is appropriate here).

Sample stability results are given in figure 2 for increasing values of the imposed electric field and fixed values of *r*=6, *m*=0.1, *ϵ*=2 and , . The band of unstable wavenumbers as well as the maximum growth rates increase as increases. The figure also includes results from direct numerical simulations starting from small initial conditions where linear theory holds, and agreement between analysis and computation is excellent (a full discussion is postponed until §3*a*).

## 3. Direct numerical simulations

Our numerical methods are based on a VOF platform adapted from the GERRIS [22] suite of algorithms. We also make use of the electrohydrodynamics toolbox available in the code architecture in order to model the action of electrodes on our boundaries. Several computational features of the package, including adaptive meshing and accurate multi-phase fluids representation, make it suitable for the present problem (see [23,24,28] and related work validating and extending the implementation, as well as testing several applications).

Here, we only highlight the numerical set-up and we point the interested reader to the work of Popinet *et al.* [22,28]. The most important numerical element in multi-fluid systems is the algorithm used for the representation of interface. In VOF methods, relevant properties such as density, viscosity or permittivity are represented in terms of a volume fraction *c*(**x**,*t*), where *c* is a generic colour function which takes the value 0 in one fluid and 1 in the other. More specifically we write
3.1
3.2
3.3
Under this treatment, a density equation of the general form
3.4
becomes
3.5
which is solved for *c* and the results substituted into (3.1)–(3.3). The value of *c* is interpolated across the interface by introducing a small transition layer in its vicinity to smooth the variation of quantities from one region to the other. Considerable efforts have been concentrated on how to best achieve this interpolation, and of particular interest to us is the interpolation across the interface of permittivity. Following numerical experimentation and to achieve added accuracy, the GERRIS package incorporates the harmonic mean interpolation formula first tested by Tomar *et al.* [29]. Comparison against analytical solutions in benchmark tests is very promising and is a direct result of several implemented features, including dynamic spatial adaptivity. Direct numerical simulations such as those presented here converge to results for sharp interface models when denser meshes are generated around the interface along with practical mathematical optimization features. The execution time is also greatly accelerated by the implementation of a quadtree structure that allows for massively parallel computations and extended processor load control. By customizing the code to include all essential features of a well-designed multi-phase platform, we consider the selected algorithms to be ideal choices for our model.

### (a) Instabilities for uniform background electric fields

We begin by using direct numerical simulation to validate the linear stability theory of §2*b*. To fix matters, we consider a square geometry (non-dimensional size 1×1), in which two immiscible fluids are separated by a horizontal interface given by
3.6
where the wavenumber *q*≥1 is an integer so that periodic disturbances have wavelengths 1/*q*, and *A*_{i} is the initial amplitude of the perturbation; for comparisons with linear theory, we typically take . Since we only perturb the interfacial shape in the computations, a small transient time is required for the system to enter the linear growth regime. Once this is reached, we use a sliding least-squares fit in time to extract the growth rate from the position of the interface, tracked at various locations throughout the simulation.

We are interested in demonstrating the destabilizing effect of increasing electric fields and choose parameters for the fluids that are appropriate for small-scale desktop experiments. Densities are taken to be , for a strongly stably stratified system, and the dimensionless viscosity is , with a viscosity ratio of *m*=0.1 between the two fluids. The dimensionless surface tension coefficient is and the gravitational acceleration coefficient reduces to 1 in our non-dimensionalization. The permittivities are taken to be and . In order to evaluate the growth rates for a variety of imposed voltage potentials, we consider and 5.0 and plot the results of growth rate versus wavenumber *k*=2*πq* in figure 2. The figure presents growth rate curves given by linear theory for the range 0≤*k*≤25, along with growth rates estimated using direct numerical simulations as explained earlier; the latter are carried out for three initial perturbation wavenumbers *k*=2*π*,4*π*,6*π* so that there are 13 direct numerical simulation produced growth rates to compare with linear theory, spanning a reasonable range of unstable wavenumbers. Symbols are used for the direct numerical simulation growth rates and the agreement is excellent. The results also show a dramatic increase in the maximum growth rate and the band of unstable wavenumbers as the applied electric field is increased. For example, doubling from 2 to 4, produces a tenfold increase in the maximum growth rate and a fivefold increase in the band of unstable wavenumbers. The impressive level of accuracy achieved by the simulations enables us to focus attention on the nonlinear regime and in particular to conduct a series of numerical experiments using imposed travelling wave electric fields on the wall electrodes in order to induce directional fluid motions.

### (b) Voltage travelling wave boundary condition

We propose to drive the multi-layer flow from the boundaries by imposing a spatio-temporally varying voltage on the lower wall electrode in the form of a travelling wave. (Practically, this could be achieved by etching a series of parallel strip electrodes on the lower substrate, and producing a travelling voltage waveform there by switching the voltage from on to off on successive electrodes (e.g. [20]).) Mathematically, we model such a forcing by imposing , where
3.7
The form (3.7) represents a uniform voltage *C* with a superimposed hump of amplitude *A*_{r} and width *x*_{R}−*x*_{L}. The transition from the background value *C* to *C*+*A*_{r} takes place in the vicinity of *x*=*x*_{L}+*U*_{r}*t* and *x*=*x*_{R}+*U*_{r}*t*, and the size of the transition layer is of size *d*>0 (small values of *d* correspond to sharp transitions). The substrate voltage is a travelling wave of constant speed *U*_{r} (right/left moving depending on whether *U*_{r} is positive/negative). Note that we can add several humps with non-equal distances between them, and extend the function periodically for consistency with the periodic boundary conditions used in our simulations. In fact, multiple such humps are constructed and their widths selected so that the substrate forces the generation of the most unstable waves provided by linear theory for a given set of parameters.

In the computations that follow, we set *d*=10^{−3}, and pick *C* sufficiently large so that instability is present (guided by linear theory (e.g. figure 2)), yet the growth rates being relatively small so as to prevent the interface touching the wall. The parameters *A*_{r} and *U*_{r} must also be tuned carefully so that the wall-forced flow can compete and indeed overcome the vertical interface touchdown motion that a uniform background field would produce. Sample numerical results of the interior voltage distribution for a substrate voltage of the form (3.7), but having three hump regions, are given in figure 3. The figure shows the equipotentials at a fixed time for the parameters *A*_{r}=1.0, *x*_{R}−*x*_{L}=0.1; it can be seen that a periodic disturbance emerges with the field being stronger above the hump regions as expected. Such a distribution of the electric field is used in the numerical results described in §3*c*, with the hump regions initially aligned with interfacial minima of the nonlinear solutions as it evolves towards touchdown with the lower wall. We also note that the electrostatic approximation used here remains valid even though we have a time-dependent electric field imposed on the lower substrate. Firstly, given that the length scales of interest are of the order of a few centimetres, it is easily established that induced magnetic fields are of the order of 10^{−10} T and hence negligible. Secondly, the period of oscillation of the imposed field is of the order of the viscous timescale *L*/*U*=*Lμ*/*σ*, which for oils and centimetre-sized systems is approximately 0.1–1 s, giving a frequency of oscillation of 1–10 Hz (electrostatic approximations hold even for frequencies in the kHz range).

The forcing introduced in this section (along with its multi-hump variants) provides a direct way of inducing fluid–structure interactions in order to drive the system far from its equilibrium state or indeed other nonlinear strongly attracting states. Of particular interest is the exploration of such mechanisms in breaking the directional symmetry of the problem and producing a net flow in one direction, without the need of moving parts. Such micro-pump systems can be useful in a host of applications and in what follows we demonstrate their feasibility using direct numerical simulations.

### (c) Numerical results

#### (i) Generation of nonlinear travelling waves using the substrate voltage

The geometry and fluid parameters used are the same as those in §2*b* and figure 2, and the main modifications enter through as discussed next. We start with a relatively strong uniform electric field of imposed on the lower electrode (note that the direction of the electric field is immaterial—we chose to place the forcing on the lower wall after preliminary numerical experiments indicated that such configurations are much more efficient computationally in producing uni-directional interfacial travelling waves), so that an initial perturbation of wavenumber *k*=6*π* is unstable. The amplitude of the perturbation quickly increases and enters the nonlinear stage with the local extrema of the interface moving towards the electrodes (figures not shown for brevity).

At *t*=0.25, we alter the imposed electric field behaviour, reducing its strength so that parameter *C*=3.5 in equation (2.1), and this will now be the strength of the background field on top of which the travelling wave boundary condition is imposed. When , figure 2 shows that the *k*=6*π* perturbation represents the most unstable wavenumber. Even though secondary harmonics are present in the problem resulting in rich nonlinear dynamics, the excitation of the interface with the most unstable wavenumber provides a computationally more efficient way of quantifying the phenomena we are investigating. Next, we impose a number of humps equal to the number of local minima in the initial perturbation and align them with the interface minima—see figure 3 for an illustration of the imposed field at this time. For the remaining two parameters at our discretion, *A*_{r} and *U*_{r}, we choose a set of four values for each, i.e. we extract pairs of values from {0.75,1.0,1.25,1.50} for both *A*_{r} and *U*_{r}. The 16 different full direct numerical simulations were carried out and the resulting data carefully post-processed to study the evolution of relevant flow quantities. The rationale behind the chosen set of values is as follows. If the amplitude *A*_{r} of the imposed voltage hump is too low, its effect is negligible and the vertical motion caused by the background electric field is sufficiently strong and was found to cause interfacial wall touchdown. On the other hand, an unnecessarily high value of *A*_{r} was found to result in violent flow phenomena, such as interface rupture and droplet formation, which are undesirable in the context of our study. The interval of appropriate values for *A*_{r} was constructed empirically via extensive numerical experiments. We will see that, surprisingly, the effect of the wave velocity *U*_{r} is of secondary importance and that the main crucial flow control quantity is the amplitude *A*_{r}.

Sample results are provided in figure 4 to depict the steps mentioned earlier. The left column shows the initial destabilizing effect of the vertical electric field with , bringing the interface close to the channel walls. After the lower wall voltage travelling wave form is turned on and allowed to have an effect on the flow, the emerging nonlinear dynamics consist of a nonlinear interfacial wave travelling to the left (snapshots at *t*=2.8,2.9,3.0 are shown in the right column from top down). The travelling wave is not of permanent form but contains asymmetric modulations along with additional higher harmonics entering and becoming nonlinear. The sequence of images presented in figure 4 represents only a small time window. A much more complete picture of the flow that includes the genesis of the travelling wave and its subsequent nonlinear dynamics is also available in animation form as electronic supplementary material. Frames are extracted every 0.0025 time units and combined into a movie of the interfacial dynamics in the interval 0≤*t*≤10. We therefore illustrate the full simulation from the start, first noticing the effect of the uniform electric field, followed by the introduction of the travelling wave voltage boundary condition that creates the horizontal velocity flux we wish to study in more detail. In this set of computations, the imposed voltage travelling wave speed is *U*_{r}=1.0 and the hump amplitudes *A*_{r}=1.0.

#### (ii) Velocity fields and the emerging net horizontal flow

The simulations of §3*c*(i) show the spatio-temporal evolution of the interface under the action of an imposed voltage travelling wave and provide strong evidence of the feasibility of creating and sustaining a large amplitude travelling wave. Such depictions do not provide details of the flow velocities throughout the domain, and we turn our attention to this next. We produce flow fields on regular grids of 256^{2} points on the square [−0.49,0.49]×[−0.49,0.49], by interpolating the computed velocities and pressure (these are computed on adaptive grids which can be much finer in places) every 0.25 time units. Of particular interest is the emergence of a net flow in the horizontal direction (to the left in these particular simulations, but we have investigated other parameters which produce right-moving travelling waves), and hence we focus our attention on the horizontal velocity *u*. In figure 5, we show the variation of the horizontal velocity at the origin *x*=0 across the channel, i.e. we track *u*(0,*y*,*t*) for the parameter values *A*_{r}=*U*_{r}=1. The different panels in figure 5 correspond to different times, starting at *t*=2 with the leftmost panel and increasing time by 2 units until *t*=10 at the rightmost panel.

There are two notable effects revealed by these results. Firstly, from the very onset of the voltage travelling wave boundary condition, a strong negative horizontal velocity is present in the flow, in the immediate vicinity of the voltage potential perturbation (the extent of the voltage potential perturbation can be observed in the simulations of figure 3). The flow is affected by the nonlinear interfacial dynamics induced by the non-uniform voltage boundary condition, and eventually the system settles down into a lateral leftward motion. Secondly, as time increases the horizontal velocity increases in magnitude in the whole domain, with substantial velocities also found away from the lower wall. By *t*=10.0 an almost uniform plug flow is seen in about 60% of the channel height, with no-slip conditions near the walls modifying these structures. This phenomenon is quite important since it shows clearly that a net horizontal flux can be set up by using electrohydrodynamic fluid–structure interactions. The induced fluxes are considered in more detail next.

#### (iii) Net horizontal fluxes as the imposed voltage travelling wave parameters vary

Using results analogous to those of figure 5, we compute the horizontal flux at *x*=0 defined by for each of the 16 configurations introduced in §3*c*(i). As noted earlier, we consider data for *F*_{0}(*t*) in increments of 0.25 starting from *t*=0 and going above *t*=10 (the electric field on the lower wall is switched off at *t*=10). The results are collectively shown in figure 6 for each pair *A*_{r},*U*_{r} as labelled on the figure. The vertical dashed-dotted line on the left indicates the time *t*=0.25 when the non-uniform voltage is switched on at the lower wall, while the dashed-dotted line on the right indicates the time *t*=10 when the voltage on the lower wall is switched off altogether. The results have been grouped and colour coded according to the value of *A*_{r}; the four curves corresponding to the smallest value *A*_{r}=0.75 are colour coded black, the *A*_{r}=1.0 curves dark grey, the *A*_{r}=1.25 curves light grey and the largest value *A*_{r}=1.5 curves are colour coded white with open symbols. This grouping of the data immediately indicates that the effect of *U*_{r} is rather weak, at least for the ranges investigated here. On the other hand, the value of *A*_{r} is of central importance, particularly regarding the generation of interfacial travelling waves. For example, the first set of four direct numerical simulations for the smallest value *A*_{r}=0.75 are interrupted before *t*=10 and in fact the longest lived computation corresponds to *U*_{r}=0.75 and was stopped at *t*=8. The reason for halting these computations is that for such a relatively small value of *A*_{r} the electrodynamic effects of the boundary travelling wave voltage potential on the flow are too weak to prevent the interface from touching the wall (this study is not concerned with interface–wall touchdown and so we do not present or analyse results in such cases). Going to larger values of *A*_{r}, we first observe that the flow is controlled in the sense that interfacial touchdown is prevented. Furthermore, as *A*_{r} increases, the magnitude of the flux |*F*_{0}| also increases with the largest fluxes attained for *A*_{r}=1.5. In all cases that do not suffer interface–wall touchdown, the flux appears to decrease linearly with time until *t*=10 when the voltage is switched off; beyond this time, the flow eventually relaxes to its original stably stratified state with a flat interface. The computations presented provide strong evidence that voltage travelling wave boundary conditions can be used to achieve three phenomena: (i) prevent interface–wall touchdown, (ii) produce nonlinear travelling wave interfacial structures, and (iii) induce and sustain a net horizontal flux thus underpinning possible applications of microfluidic pumps with no moving parts.

## 4. Conclusion

We have analytically and computationally studied the dynamics of viscous, immiscible perfect dielectric multi-layer flows confined in channel geometries, when electric fields act across the channel. We began by identifying and computing the underlying instabilities into the fully nonlinear regime when the imposed voltage potential across the channel is constant (each wall electrode is held at constant but different voltage potential). In this paper, we also used linear stability theory to evaluate our direct numerical simulation numerical algorithms and found agreement between the two to be excellent (figure 2). If the constant applied potential difference is sufficiently large, the flow evolves to a nonlinear and quasi-static configuration, and if the potential difference is increased further the interface becomes attracted to the walls and touchdown in finite time takes place (such events present significant computational challenges).

In order to generate an interfacial travelling wave, along with its associated directional net flux properties, it was necessary to carry out a large number of numerical experiments before identifying a method that works and is robust. Our computational efforts and results in achieving this are outlined in §3*c*. The procedure that seems to work best for the problem at hand is as follows: (i) switch on a uniform background electric field of sufficient strength to drive the interface towards one of the walls (the lower wall in our case), (ii) when the interface is close to the wall, reduce the constant applied voltage to a value denoted as *C* in equation (3.7) and at the same time introduce a voltage travelling wave perturbation at the boundary as illustrated in equation (3.7) and figure 3, and (iii) select the value of *A*_{r}, the size of the applied voltage perturbation humps (see (3.7)), and *U*_{r}, the speed of the imposed voltage travelling wave, so as to prevent wall touchdown. If these steps are completed successfully, it was found that an interfacial travelling wave forms, and at later times the horizontal velocity can develop into a plug flow over roughly two-thirds of the channel height (figure 5, rightmost panel). The associated flux in the computations presented is negative (the induced interfacial wave travels in the negative direction) and increases as the voltage perturbation *A*_{r} increases (figure 6). On the other hand, the dependence of *U*_{r} is rather weak, and we conclude that it does not play a central role (it needs to be an order one quantity, however, in order to induce travelling wave phenomena).

We surmise that the procedure identified here can be useful in producing pumping on the microscale without the need of any moving mechanical parts. Several directions towards optimizing such operations suggest themselves and will be investigated in the future; foremost among them is the addition of a second voltage travelling wave on the upper wall with adjusted phase relative to the travelling wave on the lower electrode. Another area of practical importance is that of mixing on the microscale, and the multi-layer configuration presented here can use the spatio-temporal dynamics of the interface to provide mixing in either phase.

## Funding statement

The work of D.T.P. was partly supported by EPSRC grant no. EP/K041134/1.

## Acknowledgements

R.C. acknowledges a Roth Doctoral Fellowship from the Department of Mathematics, Imperial College London.

## Footnotes

One contribution of 15 to a Theme Issue ‘Stability, separation and close body interactions’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.