A hydrodynamic bifurcation in electroosmotically-driven periodic flows
Alexander Morozov Davide Marenduzzo Ronald G. Larson





alexander.morozov@ph.ed.ac.uk
SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom


SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom


Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109, USA




In this paper we report a novel inertial instability that occurs in electro-osmotically driven channel flows. We 
assume that the charge motion under the influence of an externally applied electric field is confined to a small vicinity of the channel walls that, effectively, drives a bulk flow through a prescribed slip velocity at the boundaries. Here, we study spatially-periodic wall velocity modulations in a two-dimensional straight channel numerically. At low slip velocities, the bulk flow consists of a set of vortices along each wall that are left-right symmetric, while at sufficiently high slip velocities, this flow loses its stability though a supercritical bifurcation. Surprisingly, the new flow state that bifurcates from a left-right symmetric base flow has a rather strong mean component along the channel, which is similar to pressure-driven velocity profiles. The instability sets in at rather small Reynolds numbers of about 20-30, and we discuss its potential applications in microfluidic devices.


A hydrodynamic bifurcation in electroosmotically-driven periodic flows
    Ronald G. Larson
    December 30, 2023
======================================================================



In microfluidic devices, the use of electric fields as a means of driving flow via electro-osmosis is an intriguing alternative to using pressure drops or moving surfaces <cit.>.
Electro-osmosis occurs when the ions in a double layer next to a charged surface are set in motion by an electric field, and the ions drag the solvent with them, producing bulk flow. Such flows are especially suitable for microfluidic applications, in which microfabrication techniques allow for control and patterning of electric and dielectric properties of channel surfaces.  In this way, not only can bulk flow be generated to transport analytes, but patterned flow fields can be imposed, allowing, for example, for creation of microfluidic mixers <cit.>. Such flows may also assist in separating particles or cells, possibly both modulating and augmenting the inertial forces that produce size-depending cross-stream drift <cit.>.
Even if uniform charge density is intended for a surface, some variation in charge is unavoidable, especially given the difficulty in controlling precisely  the surface chemistry producing the charge, and this will give rise to non-uniform electro-osmotic flows even in a straight channel.  One of the conceptually simpler non-uniform electroosmotic flows that can be produced is generated by a sinusoidally periodic surface charge on each side of a straight channel <cit.>; see Fig.<ref>.  This charge pattern leads to a spatially periodic charge in the double layer adjacent to the boundaries. When a voltage is applied along the channel, the velocity near the surface varies periodically as well, and acts like a periodic 'slip' velocity along the surface, generating a complex cellular flow in the fluid in the channel.  This flow is attractive as a simple boundary condition (straight walls, periodic charge) that nevertheless generates a complex flow, and that, moreover, has an analytic solution in the limit of creeping flow and thin double layers <cit.>. It is easy to add a uniform surface charge density to the periodic charge, mimicking, for example, an imperfectly treated surface with charge non-uniformity. A periodic deviation from a uniform charge might produce a deviation in the mean flow rate in the channel owing to nonlinear coupling between the flow produced by the uniform charge, and that produced by the periodic charge variations.  If the wall charge varies sinusoidally around zero, the electro-osmotic flow that is generated is periodic, and in the Stokes flow limit has no net flow direction. 



Here we consider the effect of inertia on this simple flow. We employ a spectral method to solve the two-dimensional Navier-Stokes equations, and find, surprisingly, for the case of zero average surface charge, so that the Stokes flow is periodic with no net flow, that at a modest Reynolds number Re=v_0 L/ν of around 20, there is a bifurcation to a secondary flow with a non-zero mean flow, even when there is no mean flow induced by the boundary conditions themselves.  Here, v_0 is the characteristic velocity, L is the half-width of the channel, and ν is the kinematic viscosity of the fluid. The presence of this bifurcation means that, even for a boundary condition with no mean surface charge, and hence no mean current, a rectified mean flow can be produced through a purely oscillatory boundary condition.  The direction of the mean flow, to the right or the left, is arbitrary, but could be imposed by adding some small bias to the initial oscillatory flow, either electrically, geometrically, or in some other way. We believe that this is the first report of this hydrodynamic instability in a periodically patterned electroosmotic flow (although this discovery was alluded to in an earlier co-authored by one of the present authors <cit.>).  This bifurcation is of interest in its own right, but might also be a means of generating rectified flow in a channel with no net imposed current.  In fact, since the base flow is completely periodic, the applied voltage along the channel could in principle also be alternating, without changing either the zero net current, or the direction of the resulting flow. To reach the bifurcation condition with a wall charge that varies with position sinusoidally around zero, the Reynolds number must reach a value of close to 20, which, for water with ν = 10^-6m^2/s, requires a flow velocity and channel width 2L that are relatively large. The flow velocity v_0 is given by around σ_0 E/μκ, where σ_0 is the amplitude of the surface charge density, E the electric field imposed parallel to the walls, μ the fluid viscosity and κ the inverse Debye length at the wall <cit.>, where typical values are σ_0 ∼ 1 charge/nm^2, μ = 10^-3Pa s, and κ^-1 = 10nm. Under these conditions, a field of 10^4V/m would yield a velocity of 10^-2m/s, and so a channel of width 2L=2mm would suffice to yield an instability.  Note that the channel depth would need also to be comparable, or larger than, this scale, to prevent viscous suppression of the instability. 

This instability is distinct from well-known electrokinetic flow instabilities that result from coupling of electric fields and ionic conductivity gradients <cit.>, since there are no ionic conductivity gradients considered in what we report here. Such gradients can arise when a core fluid flow is focused by a 'sheath' fluid introduced at the walls of a microfluidic device, if the two fluids have differing ionic strengths. An instability occurs in such flows at a critical electric current Rayleigh number of Ra_e = (ϵ E_a^2 d^2/Dμ)((γ-1)/γ) |grad^*σ^*|_max = 205, as reported by Posner and Santiago <cit.> and Posner et al. <cit.>.  Here ϵ is the fluid permittivity, E_a is the applied electric field, d is the channel depth, D is the ion diffusivity, μ is the fluid viscosity, γ is the ratio of core-to-sheath conductivities, and |grad^*σ^*|_max is the maximum dimensionless conductivity gradient.  Note that the Rayleigh number representing the driving force for this instability disappears when the conductivities of the two fluids are equal to each other (i.e., γ = 1), since then there is no conductivity gradient to which the electric field can couple and the above Rayleigh number is zero. If a conductivity gradient is present,  the instability produced by it could in many cases occur at lower field strength than that needed to produce the inertial instability to be discussed here.  

We consider the flow of a Newtonian fluid in a 2D-channel forced by a prescribed slip velocity at the channel walls, similar to Fig.<ref>. We introduce a Cartesian coordinate system with the x-axis pointing along the length of the channel, and the y-axis – in the wall-normal direction. The walls are located at y=± L, so that the total channel width is 2L. The velocity components are v=(u(x,y,t),v(x,y,t)), and the flow is driven by the following slip velocity at the walls

    u(x,L,t) = v_+ coskx, 
        u(x,-L,t) = v_- cos(kx+ϕ), 
        v(x,± L,t) = 0.

Here, v_+ and v_- are the maximum slip velocities at the corresponding boundary, k is the wavenumber of the slip velocity modulation, and ϕ is the phase difference between the velocity at the upper and lower walls. The equations of motion are given by the Navier-Stokes equation

    ρ[ ∂ v/∂ t +  v·∇ v] = -∇ p + μ∇^2  v,

and the incompressibility condition

    ∇· v = 0.

Here, ρ and μ are the density and viscosity of the fluid, respectively, and p is the pressure. 
The problem is rendered dimensionless by the rescaling of all the variables, where 
we use the half-width of the channel L as the unit of length, max(v_+,v_-) as the unit of velocity, and L/max(v_+,v_-) as the unit of time. We also introduce the Reynolds number 

    Re=max(v_0,v_1) L/ν,

and the dimensionless wave-vector k̃ = k L. In what follows, all variables are dimensionless unless stated otherwise. 

To reduce the number of degrees of freedom, we introduce the streamfunction Ψ=Ψ(x,y,t), such that

    u = ∂Ψ/∂ y,      v = -∂Ψ/∂ x.

In terms of the streamfunction, the equation of motion is given by

    ( ∂/∂ t + ∂Ψ/∂ y∂/∂ x - ∂Ψ/∂ x∂/∂ y)∇^2 Ψ = 1/Re∇^4 Ψ,

with the following boundary conditions

    ∂Ψ/∂ y(x,1,t) = ṽ_+ cosk̃ x, 
       ∂Ψ/∂ y(x,-1,t) = ṽ_-  cos(k̃ x+ϕ),
       ∂Ψ/∂ x(x,± 1,t) = 0,

where ṽ_+=v_+/max(v_+,v_-), and ṽ_-=v_-/max(v_+,v_-), and ∇^4 is the biharmonic operator.



It is instructive to consider the limit of zero inertia and steady velocity field. In this case, Eq.(<ref>) reduces to the Stokes equation, ∇^4 Ψ = 0, that has a simple analytical solution satisfying the boundary conditions, Eqs.(<ref>),

    Ψ_0 = 2k̃/sinh^22k̃-4k̃^2[
    A(x)(1+y)sinhk̃(1-y)
                                - B(x)(1-y)sinhk̃(1+y)],

where

    A(x) = ṽ_+ cosk̃x + ṽ_- sinh2k̃/2k̃cos(k̃x+ϕ), 
    
    B(x) = ṽ_- cos(k̃x+ϕ) + ṽ_+ sinh2k̃/2k̃cosk̃x.

This solution is similar to the one obtained by Ajdari <cit.>. Although Ψ_0 differs from the true solution to Eq.(<ref>) for any finite amount of inertia, it is nevertheless useful for gaining an insight into the structure of the flow.  In Fig.<ref> we plot the velocity profile given by Ψ_0 for ṽ_+=ṽ_-=1 (top row), and ṽ_+=1, ṽ_-=1/3 (bottom row) for three values of the phase difference ϕ: 0 (left column), π/2 (middle column), and π (right column). As can be seen from the figure, the flow consists of two arrays of vortices aligned along each wall with their relative position and strength set by the phase-difference ϕ and the velocity magnitudes ṽ_+ and ṽ_-, respectively.

To assess the effect of inertia on this solution, we solve Eqs.(<ref>) and (<ref>) numerically using a Fourier-Chebyshev pseudo-spectral method <cit.>. We express the streamfunction as a Fourier series

    Ψ(x,y,t) = ∑_n=-N^Nψ_n(y,t)e^i n k̃ x,

where ψ_n(y,t) = ψ_-n^*(y,t) to ensure that Ψ(x,y,t) is real, and * denotes the complex conjugate. At any time t, ψ_n(y,t) is represented by its values at M Gauss-Lobatto points <cit.> in the wall-normal direction, and the y-derivatives are taken by multiplying these values with the Chebyshev pseudo-spectral differentiation matrix <cit.>. The non-linear terms are calculated by performing a discrete Fourier transform of the streamfunction to real space, evaluating the non-linear terms there, and performing an inverse discrete Fourier transform back to spectral space; the 3/2-rule is used to avoid aliasing errors and the boundary conditions are implemented using the tau-method <cit.>. For each set of parameters, we check convergence of the velocity field by comparing it at several resolutions (N,M); convergence was always reached for N=5 (before de-aliasing) and M=80. Most of the results presented below are obtained by using the Newton-Raphson algorithm <cit.> to solve the time-independent version of Eq.(<ref>). We also performed direct numerical simulations of Eq.(<ref>) using a fully-implicit Crank-Nicolson method <cit.>; for all parameters studied, convergence was reached for the dimensionless time-step of 10^-2.

First, we study how the presence of inertia modifies the Stokes solution, Eq.(<ref>), at relatively low Reynolds numbers. Using the Newton-Raphson method, we find steady solutions of Eq.(<ref>), and compare them to the Stokes profile Ψ_0. The difference is quantified by calculating the kinetic energy of the flow, defined as

    E = k̃/2π∫_0^2π/k̃dx  1/2∫_-1^1dy 1/2[ (∂Ψ/∂ y)^2+(∂Ψ/∂ x)^2 ],

for the inertial, E_i, and Stokes solutions, E_s. In Fig.<ref> we plot the ratio E_i/E_s for ṽ_+=ṽ_-=1, ϕ=0, and k̃=π. The data demonstrate that the inertial contribution to the kinetic energy is only about 4% of the total kinetic energy at Re=30, and that that contribution decreases for smaller values of Re.  The difference E_i/E_s-1  scales quadratically with Re (fit not shown), implying that the small inertial correction to the Stokes profile can be obtained from the leading-order term of the perturbation theory in Re, even for Re∼ 30. Visual inspection of the inertial velocity profiles together with the data in Fig.<ref> suggests that the Stokes solution, Eq.(<ref>), is a very good approximation to the actual inertial solution even at moderate Reynolds numbers. 



The situation changes significantly at higher Reynolds numbers. In Fig.<ref>(left) we plot the velocity profile for the symmetric boundary conditions at Re=40 and k̃=π, and observe that it no longer posseses a translation-reflection symmetry along the x-axis, c.f. Fig.<ref>(top, left). This is associated with the emergence of the zeroth Fourier mode U(y) of the horizontal velocity component u(x,y), see Fig.<ref>(right), absent at lower Reynolds numbers. This x-independent, mean flow along the x-direction reaches significant amplitudes of about 23% of the maximum slip velocity at the wall.



To characterise this new flow state, we introduce a dimensionless order parameter

    χ = Re k̃/2π∫_0^2π/k̃ dx ∫_-1^1 dy u(x,y) ≡
    Re ∫_-1^1 U(y) dy,

which is a two-dimensional flow rate along the channel (in physical units) scaled by the kinematic viscosity ν of the fluid. In Fig.<ref>(left) we plot χ as a function of the Reynolds number for ṽ_+=ṽ_-=1, ϕ=0, and k̃=π (black line). For low values of Re the flow is left-right-symmetric, there is no mean flow, and χ=0, while at larger Re, χ acquires non-zero values indicating the presence of a mean flow. The direction of the mean flow is selected by a spontaneous symmetry breaking, and can be in either direction along the channel. The state diagram, Fig.<ref>(left), therefore has two symmetric branches, ±χ, typical of a super-critical (pitchfork) bifurcation. By combining the Newton-Raphson and time-iteration techniques, we have verified that the left-right symmetric solution with χ=0 is also present for higher values of Re but is linearly unstable. The final flow state with χ 0 is stationary and stable with respect to small perturbations. Therefore, we conclude that the new flow state is a result of a linear instability that sets in at Re_crit≈ 33.3, for this set of parameters. In Fig.<ref>(left) we also show the bifurcation diagrams for other values of the phase-difference ϕ, and observe that the lowest Re_crit is achieved for ϕ=π; the corresponding base profile is shown in Fig.<ref>(top,right). 


The instability thresholds presented above were calculated by imposing a fixed value of k̃, i.e. assuming a particular spatial period of the solution. To find the critical condition in an infinitely long channel, we now study how Re_crit depends on k̃. In Fig.<ref>(right) we plot the non-linear stability thresholds for two values of ϕ, and observe that Re_crit^min = 21.2 for k̃=2.6 and ϕ=π. The stability thresholds for other values of ϕ lie in-between the two cases presented in Fig.<ref>(right), similar to Fig.<ref>(left).

We also studied the effect of the asymmetry in the wall slip velocity (not shown), with either ṽ_+ or ṽ_- smaller than unity. For every set of ϕ and k̃ considered, the corresponding Re_crit was found to be larger than Re_crit for ṽ_+=ṽ_-=1.

As mentioned in the Introduction, this instability can potentially be utilised as a means of creating a unidirectional flow in a microfluidic device, although relatively high transitional Reynolds numbers and the a priory unknown direction of the flow could make it impractical. We now attempt to assess whether a modification of the slip boundary condition, Eq.(<ref>), can alleviate both problems. To this end, we consider the following (dimensional) velocity profile prescribed at the walls

    u(x,± L,t) = v_c ± v coskx,

where the spatially-oscillatory component is the same as in Eq.(<ref>) for the most unstable parameters (v_+=v_-≡ v, ϕ=π, and k L=2.6), and we have introduced v_c – the amplitude of a constant slip velocity in the positive x-direction. Here, we study whether a small value of v_c can produce a significant mean flow at Re<Re_crit^min. The problem is made dimensionless as before, and we define additionally another Reynolds number, Re_c=v_c L/ν, based on the constant slip velocity. 



In the absence of the spatially-oscillatory component, equations of motion are trivially solved by a plug-like flow, U(y)=v_c (in physical units),which to the order parameter given by χ_c=2Re_c. For Re>0, we expect that the interaction between the plug-like and spatially-oscillatory components will generate flow rates enhanced beyond χ_c. In Fig.<ref> we present the bifurcation diagram for the modified boundary conditions, Eq.(<ref>), varying Re but keeping Re_c fixed to a particular value. The Re_c=0 data is the same as in Fig.<ref>(left) for ϕ=π. In the presence of a constant bias, the bifurcation diagram loses its ±χ symmetry, and we only plot the dimensionless flow rate  in the same direction as the bias. For Re_c=0.1 and Re_c=1, the flow rate is dominated by the plug-like profile at low values of Re, while at larger Re there is an enhancement of the flow rate due to the instability. The bifurcation diagram now looks like an imperfect pitch-fork bifurcation. For yet larger Re_c, the effect of the underlying instability is masked by the presence of a strong bias and only a mild enhancement is observed. While the presence of the bias clearly enhances the mean flow rate and breaks the left-right symmetry, the enhancement is mild and it remains to be seen whether there are practical advantages of generating a steady flow in a microfluidic device by a slip velocity Eq.(<ref>) instead of a stronger steady component alone.

The results developed here have several implications.  First, we find that a periodic variation in wall charge will have only a small effect on average velocity in a microfluidic device with an otherwise uniform wall charge, if the periodic component is small (or even modest) in magnitude compared to the uniform component. This conclusion seems likely to hold if the non-uniform component is irregular or non-periodic, as long as it is significantly smaller than the uniform wall charge.  Thus, surface charge in a microfluidic device does not need to be nearly perfectly uniform to achieve a uniform flow rate, whose magnitude is set by the average surface charge, a conclusion of importance in practical applications where wall charging is unlikely to be exquisitely uniform.  Secondly, if a rectified flow with a sharp onset is desired in a microfluidic device using electric fields to drive the flow,  this can be accomplished by exploiting the bifurcation described here, albeit only for rather large channel widths and heights (i.e., millimeters) and strong fields.  In addition, there may be benefit in using periodic, or nearly periodic flows for separation of particles or cells based on size or other characteristics, including separations based on inertial forces.  These inertial forces are already being exploited in pressure-driven flows to separate rare circulating tumor cells from white blood cells <cit.>.  Electroosmotic flow driven by a periodic wall charge, along with fluid inertial forces, may expand the options for improving the efficiency of such devices. We note that inertial fluid forces in pressure-driven microfluidic devices are strong enough to induced circulating Dean flows, which are of great significance for separating particles and cells. Thus, the addition of electroosmotically driven flow, combined with inertial effects, opens multiple new opportunities for separations. Thirdly, the flows generated by periodic charges may provide a good experimental test of one's ability to control elecroosmotic flow fields, and of the ability to created controlled charge at walls.  Since the flow field is readily predicted, including the effect of surface charge amplitude and other parameters,  a measurement of the flow (even without the bifurcation) could be used to validate methods of controlling surface charge, for example. Fourthly, both the circulating primary flow and the secondary bifurcation flow described here occurs in a geometry of trivial simplicity  (a straight channel), which allows it to be used as a test flow field for exploring various advanced simulating methods, such as mesoscopic flow simulations <cit.>, and for exploring the behaviour of complex fluids in complex flows, but with simple geometry and boundary conditions <cit.>. Finally, the flow is essentially completely viscous prior to the bifurcation and described by an analytical solution to the Stokes equation, and thus represents a particularly simple and elegant example of a classical forward bifurcation at a very modest Reynolds number, and is the simplest bifurcation so far presented for electroosmotic flow.  

R.G.L. acknowledges support and hospitality of the Higgs Centre for Theoretical Physics, University of Edinburgh where a part of this work was performed. A.M. acknowledges support from the UK Engineering and Physical Sciences Research Council (EP/I004262/1). Research outputs generated through the EPSRC grant EP/I004262/1 can be found at http://dx.doi.org/xxx-xxx.



