2.1. System description

Our system consists of two superimposed immiscible fluids, confined between two parallel plates, located at $y=-H$ and $y=H$ , as depicted schematically in figure 1. The $x$ axis runs along the channel. We denote the fluid on top as fluid 1 and that at the bottom as fluid 2. The viscosity, density and permittivity of the fluids are given by ${\it\eta}_{i}$ , ${\it\rho}_{i}$ and ${\it\varepsilon}_{i}$ respectively, where $i=1$ denotes fluid 1 and $i=2$ denotes fluid 2. We define the following parameters for the fluid property ratios: viscosity ratio ${\it\lambda}={\it\eta}_{1}/{\it\eta}_{2}$ , permittivity ratio ${\it\varepsilon}={\it\varepsilon}_{1}/{\it\varepsilon}_{2}$ and density ratio $\bar{{\it\rho}}={\it\rho}_{1}/{\it\rho}_{2}$ . We assume both the fluids to be perfect dielectrics. Initially, we consider the interface between the fluids to be flat and located along $y=0$ . The surface tension at the interface is assumed to be ${\it\sigma}$ .

Figure 1.

Schematic of the system under consideration: two fluids sandwiched between two parallel plates. An electric field is applied in the axial direction, which actuates the flow. The interface originally lies flat at $y=0$ . The potentials on the walls are modulated. The potentials on the two walls are, in general, not equal, unless $m=1$ , which creates an asymmetry in the domain.

Both the walls bear a potential, which is given for the bottom wall as ${\it\zeta}_{b}(x)={\it\zeta}_{0}(n+\cos (qx))$ and for the top wall as ${\it\zeta}_{t}(x)=m{\it\zeta}_{0}(n+\cos (qx+{\it\theta}))$ , where ${\it\theta}$ is the phase shift angle between the potentials at the two walls. We assume both $n$ and $m$ to be $O(1)$ constants. In the present work we consider the existence of EDLs near both the walls, which can be generated by mixing electrolytes with the liquids. We assume the channel height to be large enough so that the thin-EDL approximation remains valid. An axial electric field, with potential drop $V$ per modulation wavelength ( $d=2{\rm\pi}/q$ ), is applied along the channel, to actuate electro-osmotic flow. In general, the flow thus created will tend to deform the interface due to the presence of axially periodic components in the potential, which, in turn, will influence the flow. In the present study, we intend to find the deformed shape of the interface as well as the resulting flow patterns for the aforementioned system configurations. Since we are going to consider slow motion of liquids through narrow confinements, the deformation in the present case will be a strong function of the capillary number, defined as $\mathit{Ca}={\it\eta}_{2}u_{c}/{\it\sigma}$ , where $u_{c}$ is the characteristic velocity and ${\it\sigma}$ is the coefficient of surface tension.

Before we move on to the detailed formulation, it is important to mention the associated assumptions made to execute the analysis. First, we assume the thin-EDL limit to be valid (as mentioned in the previous paragraph), so that the Smoluchowski slip velocity conditions can be applied at the walls. In order for the slip velocity condition to be valid, we must meet several criteria, which include low to moderate value of the ionic Péclet number, low surface potential and low surface conductance (Squires & Bazant Reference Squires and Bazant2004), and we assume that all these conditions are satisfied for the present analysis. Second, we assume the Reynolds number (Re) to be small enough that the inertial terms in the Navier–Stokes equations can be safely neglected for the analytical solution. Third, we assume that the time scale of deformation for the interface is much larger than the time scale for the flow to reach steady state. As we show subsequently, the ratio of these two scales is given by the Ohnesorge number (Mandal, Chaudhury & Chakraborty Reference Mandal, Chaudhury and Chakraborty2014) (defined as $\mathit{Oh}={\it\eta}_{2}/\sqrt{{\it\rho}_{2}{\it\sigma}H}$ ). Therefore, in the present study we assume $\mathit{Oh}^{2}\gg 1$ . This assumption enables us to resolve the evolution of the interface in time, although the transience in the flow cannot be captured under these conditions, which, in other words, primarily indicates that the fluid instantaneously adjusts to any changes in the interface shape with time. Finally, we will assume that the gravitational effects are small enough to be neglected from the analysis.

We first present the general governing equations and boundary conditions for the system as described above. We then obtain approximate analytical solutions from these equations for the deformation and the flow structure using asymptotic methods (domain perturbation) in the limit of small deformation $(\mathit{Ca}\ll 1)$ . We further assume ${\it\theta}=0$ for the analytical solutions. In the next section we attempt to solve for the flow structure and the interface shape numerically over a wide range of governing parameters ( $\mathit{Ca},{\it\lambda},{\it\varepsilon}$ and $q$ ), albeit within the framework of low Re and the thin-EDL limit. Towards this end, we employ a diffuse interface model (Badalassi et al. Reference Badalassi, Ceniceros and Banerjee2003) (phase field model), where an order parameter is defined to denote the two fluids. We then solve the coupled Cahn–Hilliard–Navier–Stokes equations (Badalassi et al. Reference Badalassi, Ceniceros and Banerjee2003) to obtain the distribution of the order parameter and flow profile, which subsequently gives us the interface shape. In § 5, we compare our approximate analytic solutions with the numerical solutions for small values of the capillary number and execute a parametric study of the parameters responsible for changes in the flow dynamics. Finally, we conclude this work in § 6 and discuss possible future extensions.

2.2. Governing equations

The equations of motion for the two fluids in the limit of thin EDL are given by the Navier–Stokes equation in the following form:

where $\boldsymbol{v}_{\boldsymbol{i}}=(u_{i},v_{i})$ represents the velocity field. The subscript ‘ $i$ ’ denotes the $i$ th fluid. The fluid velocities also have to satisfy the continuity equation in the following form:

The potential for the axially applied electric field satisfies the following equation (Ajdari Reference Ajdari1995):

In (2.3), ${\it\varphi}_{i}$ denotes the potential in the $i$ th fluid. The axial electric field is thus given by

The fluid velocities and the potential are subject to the following boundary conditions. At the walls the applied potential will satisfy (Ajdari Reference Ajdari1995)

At the walls the $x$ velocities satisfy the Smoluchowski slip velocity condition (Ajdari Reference Ajdari1996), as attributable to the thin-EDL regime we consider here. Therefore,

At the interface, the tangential components of the axial electric field and the normal component of the displacement field will be continuous, since we assume the fluids to be perfect dielectrics. On the other hand, the velocities have to satisfy the no-slip and kinematic conditions, as well as the tangential and normal stress balance criterion (Leal Reference Leal2007). Since the interface shape is a priori unknown and it will be deformed in general due to the imbalance of normal stresses, we assume that the interface shape is given by the following equation:

Since we consider the transient interface deformation, the interface height is taken as a function of both the $x$ coordinate and time. In the present case, the aforementioned boundary conditions take the following form at the interface.

Continuity of tangential electric field:

Continuity of normal displacement field:

No-slip condition:

Kinematic condition:

Shear stress balance:

Normal stress balance:

Along the axial direction, the applied potential ${\it\varphi}_{i}$ will satisfy the following conditions:

In (2.4k ) $d$ is the modulation wavelength, given by $d=2{\rm\pi}/q$ . This condition basically indicates the potential drop along a single modulation. In conditions (2.4e )–(2.4j ), $\hat{\boldsymbol{n}}$ (directed from fluid 2 to fluid 1) and $\hat{\boldsymbol{t}}$ denote normal and tangent unit vectors to the interface respectively. Additionally, ${\bf\tau}_{i}$ denotes the stress tensor in the $i$ th fluid, ${\it\sigma}$ is the surface tension and $\boldsymbol{{\rm\nabla}}_{\boldsymbol{s}}$ denotes the surface gradient operator, which has the following expression: $\boldsymbol{{\rm\nabla}}_{\boldsymbol{s}}=\boldsymbol{{\rm\nabla}}-\hat{\boldsymbol{n}}(\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})$ . The rest of the symbols bear their usual meanings. It should be noted that the stress tensors include contributions from Maxwell as well as viscous stresses, expressed as

In (2.5), $\boldsymbol{E}$ is the applied electric field given by (2.3b ). We shall denote the hydrodynamic part of the stress by ${\bf\tau}^{(\boldsymbol{h})}$ and the Maxwell stress by ${\bf\tau}^{(\boldsymbol{e})}$ . In view of the boundary conditions at the walls and taking into account the fact that the surface potentials are periodic in nature along the axial directions, we consider all the quantities to be periodic along the $x$ axis, which remains a standard assumption in EOF in the presence of charge-modulated surface conditions (Ajdari Reference Ajdari1996; Ghosh & Chakraborty Reference Ghosh and Chakraborty2012).

We now non-dimensionalize the governing equations, where the characteristic scale for each variable is used as the reference value so as to make the non-dimensional variables $O(1)$ . Towards this end, we define, for any variable ${\it\chi}$ , $\bar{{\it\chi}}={\it\chi}/{\it\chi}_{c}$ , where ${\it\chi}_{c}$ is the characteristic value. As characteristic values, we chose the characteristic electric field magnitude as $E_{0}=V/d$ , the characteristic velocity as $u_{c}={\it\varepsilon}_{2}E_{0}{\it\zeta}_{0}/{\it\eta}_{2}$ , the characteristic lengths as $x_{c}=y_{c}=H$ and the characteristic potential as ${\it\varphi}_{c}=V$ . The characteristic hydrodynamic stress is chosen as ${\it\tau}_{c}^{(h)}={\it\eta}_{2}u_{c}/H$ and $p_{c}={\it\tau}_{c}^{(h)}$ , while the characteristic Maxwell stress is chosen as ${\it\tau}_{c}^{(e)}={\it\varepsilon}_{2}E_{0}^{2}$ . The ratio of the two stresses is taken as ${\it\tau}_{c}^{(e)}/{\it\tau}_{c}^{(e)}=(E_{0}H)/{\it\zeta}_{0}={\it\beta}$ . With the above choices of characteristic variables, the non-dimensional electric field is defined as

Previous studies (see Mandal et al. Reference Mandal, Chaudhury and Chakraborty2014 and the references therein) have demonstrated that the time scale for deformation of the interface is of the order of $t_{c}\sim \mathit{Ca}(H/u_{c})$ ( $\mathit{Ca}$ is the capillary number, as defined earlier). Therefore, we chose this scale as the characteristic value of the time. It can be easily understood that for $\mathit{Ca}\ll 1$ , this scale is much smaller than the advective time scale for fluid flows. We also note that we have chosen the length scales in both directions to be the same, which indicates that the modulation wavelength ( $d=2{\rm\pi}/q$ ) must be of the same order as the channel height (or half-height). Enforcing the aforementioned non-dimensionalization scheme, we obtain the following non-dimensional equations:

along with the continuity equation in the form

The potentials satisfy

It should be noted that we have used the same symbols for the non-dimensional variables as the dimensional ones, since (unless otherwise mentioned) we will be using the non-dimensional variables throughout the present study. The boundary conditions take the following forms:

In (2.8e ) and (2.8f ), ${\bf\tau}_{i}={\bf\tau}_{i}^{(\boldsymbol{h})}+{\it\beta}{\bf\tau}_{\boldsymbol{i}}^{(\boldsymbol{e})}$ , for $i=1,2$ . As mentioned earlier, in the present study we consider the cases where $\mathit{Oh}\gg 1$ and $\mathit{Re}\ll 1$ . Enforcing these approximations, (2.7a ) and (2.7b ) become The rest of the equations (i.e. (2.7c )) and the boundary conditions (i.e. (2.8d )–(2.8f )) remain unchanged. It should be noted that the deformation $h(x,t)$ is not known a priori and has to be determined as part of the solution.

3.1. Leading-order solutions

The leading-order velocities can be obtained assuming a flat interface, and the stresses arising from these velocities will give the deformation at $O(\mathit{Ca})$ . Enforcing (3.1) in (2.9a )–(2.9b ) and (2.7c,d ) one obtains the leading-order equations of motion and electric field, which are given by

The boundary conditions are given by (following (2.4b ), (2.4c ) and (2.8a–f ))

The potential additionally satisfies

In (3.4b ), $\boldsymbol{e}_{\boldsymbol{x}}$ is the unit vector along the $x$ axis. ${\it\tau}_{yy,1}^{(0)}=2{\it\lambda}(\partial v_{1}^{(0)}/\partial y)+{\it\varepsilon}{\it\beta}$   $((E_{y,1}^{(0)})^{2}-(E_{x,1}^{(0)})^{2})/2$ and ${\it\tau}_{yy,2}^{(0)}=2(\partial v_{2}^{(0)}/\partial y)+{\it\beta}((E_{y,2}^{(0)})^{2}-(E_{x,2}^{(0)})^{2})/2$ , including contributions from hydrodynamic as well as Maxwell stresses. Likewise, in the shear stress balance equation (3.4g ), we have also included the contribution from the Maxwell as well as viscous stresses. From (3.3d ), (3.4a ), (3.4c ) and (3.4d ) we note that the potentials ${\it\varphi}_{i}^{(0)}$ can be solved for independently without the fluid velocities. The solution is simply

The electric field in the leading order is thus given by

In view of the boundary conditions (3.4b ) and the linearity of the governing equations, one can assume a solution for the velocity profiles in the following form:

In (3.7a ), ${\it\psi}$ is the stream function. Then, the equations of motion (3.3a )–(3.3c ) will give the biharmonic ( ${\rm\nabla}^{4}$ ) of the stream function as zero, which can be solved very easily with separation of variables assuming periodic solution along $x$ . Solution for $\hat{u} _{k}$ is even simpler, and, with the aid of conditions (3.4b,e–g ), can be obtained in the following form:

The solution for the stream function takes the following form:

where

and

The constants $a_{1}^{(0)},b_{1}^{(0)},\dots ,$  etc. can be found by enforcing the boundary conditions (3.4b,e–g ) as functions of ${\it\lambda}$ , $q$ , $m$ , ${\it\varepsilon}$ and $\partial h_{1}/\partial t$ . In (3.4g ), the expressions for the electric fields can be obtained from (3.6). Since, at the leading order, the electric field in the $y$ direction is zero, it is evident from (3.4g ) that the Maxwell stresses will have no contribution in the shear stress balance. The constants ( $a_{1}^{(0)},\,b_{1}^{(0)},\dots ,$  etc.) thus obtained have algebraically complicated expressions, and hence we do not mention them explicitly here. The MAPLE and MATLAB files containing the detailed expressions can be made available upon request. The solutions (3.7a )–(3.7e ) and the condition (3.4h ) suggest that the first-order deformation must have the following form:

Enforcing (3.8) in the solutions and in (3.4h ) results in a first-order initial value problem (IVP) for $c_{1}(t)$ , with the initial condition $c_{1}(0)=0$ . The final solution for $c_{1}(t)$ then takes the following form:

The explicit functional forms of $S_{2}$ and $S_{1}$ have been given in appendix A. It should be noted that, in steady state, $\partial h_{1}/\partial t=(\text{d}c_{1}/\text{d}t)\sin (qx)=0$ and $c_{1}(t)=S_{2}(q,{\it\lambda},m,{\it\varepsilon})$ . It is interesting to note that the Maxwell stresses at $O(1)$ contribute nothing to the deformation, through (3.4h ). Their sole effect at $O(1)$ is to cause a constant static pressure jump across the interface, which is given by

In (3.10), ${\rm\Delta}p_{s}^{(0)}=p_{s,1}^{(0)}-p_{s,2}^{(0)}$ is the static pressure jump across the interface at $O(1)$ . It should be noted that, for ${\it\varepsilon}>1$ , this jump causes higher pressure in fluid 2, whereas for ${\it\varepsilon}<1$ , fluid 1 has higher static pressure. For ${\it\varepsilon}=1$ , there is no jump in the static pressure. However, the permittivity ratio does alter the steady-state deformation, through the viscous stresses, as is evident from conditions (3.4b ). Equation (3.9) indicates that the deformation varies exponentially with time. Previous studies on interface deformation, for example deformation of droplets (Leal Reference Leal2007; Ramachandran & Leal Reference Ramachandran and Leal2012), have also demonstrated similar variation in time. Equation (3.9) further reveals that the deformation of the interface occurs at a 90° phase difference from the applied modulation in the surface potential. As expected, the deformation is a strong function of the viscosity ratio ( ${\it\lambda}$ ), the modulation frequency ( $q$ ), the permittivity ratio ( ${\it\varepsilon}$ ) and the asymmetry factor ( $m$ ). However, the time taken by the interface to reach steady state is only a function of ${\it\lambda}$ and $q$ ; the other parameters do not affect this time scale. A final point to note from (3.9) is that the $x$ -independent part of the applied potential ( $n$ , the axially constant portion) and the velocity resulting from it $(\hat{u} _{i})$ contribute nothing to the deformation at $O(\mathit{Ca})$ . Therefore, in cases of flows of superimposed fluids with uniform surface potential, no deformations occur. However, we subsequently show that the $O(\mathit{Ca}^{2})$ deformation does depend on the axially invariant part of the applied surface potential.

3.2. The $O(\mathit{Ca})$ solutions

The $O(\mathit{Ca})$ equations of motions can be easily obtained from (2.7c )–(2.7d ) and (2.9a )–(2.9b ), and read

The electric field at $O(\mathit{Ca})$ is given by

The boundary conditions for $O(\mathit{Ca})$ at the interface can be obtained from (2.8a )–(2.8f ) by expanding them in Taylor series around $y=0$ (which is the norm for domain perturbation). The conditions at $O(\mathit{Ca})$ hence obtained read

The normal stress balance at $O(\mathit{Ca})$ can be greatly simplified using the condition (3.4g ) (i.e. shear stress balance at $O(1)$ ) and the momentum equation at leading order. The simplified equation reads

We again note that, in (3.12f ) and (3.12g ), the stresses are the total stresses, including contributions from viscous as well as Maxwell stresses. The Maxwell stresses at $O(\mathit{Ca})$ are expressed as follows:

In (3.13), the electric fields can be evaluated from (3.11e ) and ${\it\delta}$ is the Kronecker delta. It should be noted that we have used the fact that $E_{y,i}^{(0)}=0$ , to write (3.13). The boundary conditions at the walls are given by (as obtained from the conditions (2.4a )–(2.4c ))

We note that (3.11a )–(3.11d ) along with the boundary conditions (3.12a )–(3.12f ) and (3.14a )–(3.14b ) are sufficient to give the full solutions for the velocity components in terms of $\partial h_{2}/\partial t$ , which can then be used in the normal stress balance (3.12g ) to obtain the $O(\mathit{Ca}^{2})$ deformation, i.e. $h_{2}$ . Similarly to the solutions at $O(1)$ , we find that the potential at $O(\mathit{Ca})$ can be calculated independently without the velocity field, from (3.11d ), subject to (3.12a,b ) and (3.14a ). The $O(\mathit{Ca})$ correction arises solely because of the deformation of the interface and the permittivity difference of the two fluids. Accordingly, the slip velocity condition at the wall also gets modified as shown in (3.14b ). Therefore, the deformed interface induces additional slip velocity at the walls, modifying the $O(1)$ Smoluchowski velocity in the process. The potentials at $O(\mathit{Ca})$ have the following solutions: For solving the potential, the $O(\mathit{Ca})$ corrections have been assumed to be periodic along the $x$ direction. The resulting electric fields can be calculated from (3.11e ). Interestingly, solutions (3.15) show that the $O(\mathit{Ca})$ potential and the corresponding electric field vary with time. Therefore, quite naturally, the slip velocity at the walls at $O(\mathit{Ca})$ , as given by (3.14b ), will also vary exponentially with time. Since the initial deformation is zero, at $t=0$ , this correction at $O(\mathit{Ca})$ remains zero and grows to steady state within a time $t\sim O(S_{1})$ . It should be noted that, for ${\it\varepsilon}=1$ , the expressions in (3.15) become exactly zero. Therefore, the $O(\mathit{Ca})$ corrections to the axial potential only exist when the permittivities of the two fluid are different.

On substituting the leading-order solutions and corresponding expressions for stresses in the condition (3.12f ), one can easily note that the contributions from the Maxwell stresses vanish exactly at $O(\mathit{Ca})$ as well. The conditions (3.12c–f ) further show that the $x$ velocities at $O(\mathit{Ca})$ must have the $\sin (2qx)$ and $\sin (qx)$ components. However, simultaneously, putting the Maxwell stresses at $O(\mathit{Ca})$ in the normal stress balance condition (3.12g ) demonstrates that the deformation must have a component of the form $g(t)\text{sin}(qx)$ . Therefore, from conditions (3.12c,d ) we can predict that along with the $\sin (2qx)$ and $\sin (qx)$ components, the $x$ velocities would also have a $\cos (qx)$ component. The solution for velocities at $O(\mathit{Ca})$ can be expressed in the following form (the solution can be obtained using a stream function as explained in the previous subsection, § 3.1):

It should be noted that in (3.16b ), the first part of the solution (i.e. $G_{k}^{(1)}(y)\sin (qx)$ ) is the contribution of the axially invariant part of the surface potential, while the $F_{k}^{(1)}(y)\sin (2qx)$ part is the contribution from the axially varying periodic component of the surface potential. The remaining part of the stream function is the contribution from the Maxwell stresses in the normal stress balance at $O(\mathit{Ca})$ . The solution at $O(\mathit{Ca})$ is completely periodic and has no net flow component. The various constant coefficients of the hyperbolic terms (i.e. $k_{1}^{(1)},\,l_{1}^{(1)},a_{1}^{(1)},b_{1}^{(1)},\tilde{a}_{1}^{(1)},\tilde{b}_{1}^{(1)},\dots ,$  etc.) can be found as functions of $\partial h_{2}/\partial t$ by enforcing the boundary conditions (3.12c )–(3.12f ). Although the equations thus obtained are linear in these coefficients, their solutions take algebraically complicated forms (as in the previous section), and therefore we do not explicitly mention those solutions here. Again, the MAPLE and MATLAB files containing the detailed expressions for all of these constants can be made available upon request. The boundary conditions (3.12c )–(3.12f ) suggest that the deformation $h_{2}$ must be of the form

It is important to note that $c_{4}(t)$ is the contribution of the axially invariant component of the surface potential, while $c_{2}(t)$ is the contribution from the axially periodic part of the surface potential. The $c_{3}(t)$ part is caused by the action of the normal component of the Maxwell stresses at $O(\mathit{Ca})$ . It is important to mention here that the constants $\tilde{a}_{1}^{(1)},\tilde{b}_{1}^{(1)},\dots ,$  etc., which are associated with the velocity component resulting from $c_{3}(t)$ , are of the following form:

Since at steady state $\text{d}c_{3}/\text{d}t=0$ , all of these constants in (3.18) become zero. Therefore, at steady state the velocity arising from the $c_{3}(t)\text{sin}(qx)$ part of the deformation vanishes. Noting that $c_{3}(t)$ comes from the Maxwell stresses, we infer that the electrical stresses only induce a transient flow, which vanishes at steady state, inducing a permanent deformation in the process. In fact, previously Yeoh et al. (Reference Yeoh, Xu and Basaran2007) have shown that the Maxwell stresses can lead to a deformed equilibrium shape of a liquid–gas interface, in the presence of a modulated potential. We thus conclude that at $O(\mathit{Ca}^{2})$ of deformation (note that $O(\mathit{Ca})$ velocities lead to the $O(\mathit{Ca}^{2})$ deformation), both the components of the surface potential as well as the Maxwell stresses contribute, as claimed in the previous section. Following the procedure outlined in the previous section, enforcing (3.17) in the solutions (3.12c ), (3.12d ) and (3.12g ) results in first-order IVPs for $c_{2}(t)$ , $c_{3}(t)$ and $c_{4}(t)$ , with the initial conditions $c_{2}(0)=c_{3}(0)=c_{4}(t)=0$ . The respective solutions take the following form:

In (3.19b ), $K=((S_{2}{\it\beta})/q)(({\it\varepsilon}-1)^{2}/({\it\varepsilon}+1))\coth (q)$ . It thus follows that the part of the deformation corresponding to (3.19b ) is always out of phase with the $O(\mathit{Ca})$ deformation (3.8) for any choice of permittivity ratio and hence tends to suppress the leading-order deformation. However, the overall effect of the permittivity ratio is more involved, since the functions $S_{6}$ , $E_{2}$ , etc. also depend strongly on the permittivity ratio. We further note that the deformation from (3.19b ) vanishes when ${\it\varepsilon}=1$ . Equations (3.19a )–(3.19c ) indicate that the steady-state deformation at $O(\mathit{Ca}^{2})$ becomes

Interestingly, the time variation of the deformation at $O(\mathit{Ca}^{2})$ is not of simple exponential type, as seen in the previous order. Rather, there are two distinct scales involved, $S_{1}$ and $S_{3}$ , which in general are not equal. These two time scales govern the relaxation of the interface to final steady state. Depending on the relative magnitude of the coefficients present before exponential terms, it is possible to have non-monotonic behaviour at $O(\mathit{Ca}^{2})$ , although the total deformation is expected to vary monotonically in time, since the leading-order deformation (i.e. $h_{1}$ ) is the most dominant term and is monotonic in time. One can thus combine the results of the two subsections and write the steady-state deformation in the following way:

One can proceed further and find higher-order velocity fields and the corresponding deformations using similar iterative techniques to those outlined above. However, the calculations become increasingly tedious as we attempt to evaluate the $O(\mathit{Ca}^{2})$ velocity fields and $O(\mathit{Ca}^{3})$ deformations, and so on. Therefore, in the present analysis we restrict ourselves to $O(\mathit{Ca})$ velocity fields and $O(\mathit{Ca}^{2})$ deformations.

4.1. Governing equations with the phase field formalism

4.2. Numerical solution procedure

5.1. Deformation at steady state

Our first focus is on the nature of deformation in steady state, for a specified range of the aforementioned parameters. To this end, we first begin our discussion with figure 3(a–d), where the variation in the interface deformation ( $h(x,t)$ ) at steady state has been depicted as a function of $x$ , for different values of Ca (0.1, 0.2 and 0.3) in 3(a), different values of ${\it\lambda}$ (0.1, 1 and 10) in 3(b), different values of $q$ ( $2{\rm\pi}/3$ , $2{\rm\pi}/4$ and $2{\rm\pi}/5$ ) in 3(c) and different values of ${\it\alpha}$ (0.7/1.3, 1, 1.3/0.7) in 3(d); the other relevant parameters are given in the caption. We have plotted analytical results in figure 3(a–c) as obtained from (3.21) along with the numerical solutions for comparison. We note that in figure 3(a–c), reasonably good agreement is observed between the numerical and the analytical results, the agreement, quite intuitively, being better for lower values of Ca. Figure 3(a) demonstrates that for small capillary numbers (indicated by $\mathit{Ca}=0.1$ ), the deformation largely resembles a simple sinusoidal curve, while increasing the capillary number distorts the shape, increasing the magnitude of the deformation in the process.

Figure 3.

(a) Plot of the steady-state deformation ( $h_{\infty }$ ) versus $x$ , for three different values of $\mathit{Ca}=0.1$ , 0.2 and 0.3, while $q=2{\rm\pi}/3$ and ${\it\lambda}=1$ . Both the analytical (as obtained from (3.21)) and the numerical solutions (contour of ${\it\phi}=0$ ) have been plotted for comparison. (b) Plot (both analytical solution and numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for three different values of ${\it\lambda}=0.1$ , 1 and 10, while $\mathit{Ca}=0.1$ and $q=2{\rm\pi}/3$ . (c) Plot (both analytical solution and numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for three different values of $q=2{\rm\pi}/3,2{\rm\pi}/4,2{\rm\pi}/5$ , while $\mathit{Ca}=0.05$ and ${\it\lambda}=1$ . (d) Plot (only numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for three different values of ${\it\alpha}=0.7/1.3,1,1.3/0.7$ , while $\mathit{Ca}=0.1$ , $q=2{\rm\pi}/3$ , ${\it\lambda}=1$ . The other parameters are ${\it\varepsilon}=1$ , ${\it\beta}=1$ , $m=n=0.2$ and ${\it\theta}=0$ .

This can be explained very easily by considering (3.21), which indicates that in the leading order (i.e. the $O(\mathit{Ca})$ deformation) the deformation is indeed given by a simple sine curve. Since, for small values of Ca, the leading-order terms dominate, the deformation follows the shape of the curve $h_{1}(x,t)$ . However, as the value of Ca is augmented, the higher-order terms (i.e. the $O(\mathit{Ca}^{2})$ deformation) start to contribute to the deformation, which distorts the shape of the interface.

Figure 3(b) reveals the effects of the viscosity ratio on the deformation of the interface. It can be noted that increasing the viscosity ratio ( ${\it\lambda}$ ) results in an increase in the magnitude of the deformation. In a strict sense, the asymptotic approximations were derived with the assumption that ${\it\lambda}\sim O(1)$ , which is violated for both ${\it\lambda}=10$ and 0.1. However, the present figure indicates reasonably good agreement between the analytical and numerical results for ${\it\lambda}=0.1$ and 1, while the matching is relatively weaker for ${\it\lambda}=10$ . Nevertheless, the analytical results do predict that the deformation increases with ${\it\lambda}$ , which is in tandem with the numerical results. This particular behaviour of the deformation can be explained based on the fact that ${\it\lambda}>1$ means that the fluid at the top is more viscous. Noting that $m=0.2$ for the present figure, which decreases the electro-osmotic forcing on the top fluid, an augmented viscosity results in further decrease in the velocity, as described in (2.4b ). This finally results in a higher degree of asymmetry between the top and bottom halves of the channel, which in turn leads to increased deformation. However, for ${\it\lambda}<1$ , the decreased forcing on the top fluid (for $m<1$ ) is compensated to some extent by a decrease in the viscosity, which arrests the asymmetry, thus leading to lower deformation. We further note that a reverse trend in the deformation is expected for $m>1$ . Therefore, the effect of ${\it\lambda}$ on the extent of deformation is closely related to the value of $m$ . Figure 3(c) shows that the deformation increases when the wavelength of modulation is increased. This increase is quite drastic for a relatively small change in the wavelength: from 3 to 5. The increase can be very easily explained based on the fact that the curvature of the interface to $O(\mathit{Ca}^{2})$ is given by the relation $C=\boldsymbol{{\rm\nabla}}_{\boldsymbol{s}}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=-\mathit{Ca}(\partial ^{2}h_{1}/\partial x^{2})-\mathit{Ca}^{2}(\partial ^{2}h_{2}/\partial x^{2})+O(\mathit{Ca}^{3})$ . Thus, the curvature decreases as we decrease $q$ or increase the wavelength of modulation. Therefore, it follows that in order to balance the jump in the normal stresses across the interface through surface tension, the interface must have a larger deformation. This is exactly what is observed in the figure under consideration.

In figure 3(d) the effect of the relative heights of the liquid layers in terms of ${\it\alpha}$ is presented using numerical simulations only. Interestingly, it is seen that with a decrease in height of fluid 1 (top layer) there is a significant decrease in the steady-state deformation of the interface. This is due to the fact that with the decrease in height of fluid 1, there is an increase in velocity and subsequently an increase in the stresses acting on the interface. As the electro-osmotic forcing is less on fluid 1, the asymmetry between the top and bottom halves of the channel reduces and the deformation of the interface also reduces as a result of this. A reverse trend is observed when there is an increase in the height of fluid 1. We note that the effect of the parameter ${\it\alpha}$ on the deformation is also dependent on the magnitude of $m$ .

Figure 4.

Plot (both analytical solution and numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for two different values of ${\it\varepsilon}=0.5$ and ${\it\varepsilon}=3$ . The other parameters have the following values: $m=n=0.2$ , ${\it\lambda}=1$ , $q=2{\rm\pi}/3$ , ${\it\beta}=1$ , ${\it\theta}=0$ , $\mathit{Ca}=0.1$ and $\mathit{Cn}=0.02$ .

Figure 4 depicts the variation in the steady-state deformation, as a function of $x$ , for two different values of the permittivity ratio, ${\it\varepsilon}=0.5$ and 3; the other relevant parameters are given in the caption. Both analytical and numerical solutions have been plotted, and good agreement is observed between them. We note that a larger value of the permittivity ratio results in smaller deformation, for the particular choice of $m=0.2$ , as evident from numerical as well as analytical results. A simple explanation for this can be provided by the fact that the slip velocity for the top fluid at the wall is ${\sim}{\it\varepsilon}m/{\it\lambda}$ . Therefore, for $m<1$ , the choice of ${\it\varepsilon}<1$ further increases the asymmetry between the top and bottom halves of the channel. Since the deformation in the present study is in large part driven by asymmetry between the bottom and top halves of the channel, the choice of ${\it\varepsilon}<1$ would naturally lead to larger deformation. This same reasoning also explains the decrement in $h_{\infty }$ when we chose ${\it\varepsilon}=3>1$ . It is important to note here that the effect of ${\it\varepsilon}$ is threefold in the present study. First, it affects the deformation by controlling the Smoluchowski slip velocity, as depicted in the boundary conditions (2.4b ) and (2.4c ). Second, it changes the magnitude of the axial electric field across the two fluids through (2.3a ), thereby effecting the deformation. Finally, it also contributes to the deformation through the Maxwell stresses at the interface.

We next move to figure 5, where we show that, even for very small values of Ca, fairly large deformations can be obtained. To this end, in figure 5(a), we plot the deformation for $m=2$ and 5 with $\mathit{Ca}=0.01$ and ${\it\varepsilon}=2$ , while in figure 5(b), we depict the deformation for $\mathit{Ca}=0.005$ and ${\it\varepsilon}=3$ . All other relevant parameters are given in the caption. This figure clearly depicts that, even at very low values of $\mathit{Ca}\sim O(10^{-2}{-}10^{-3})$ , it is possible to have large deformations, to the tune of $h_{\infty }\sim 0.3$ . We further observe that both the numerical and analytical solutions predict similar results and they agree well in both the cases. A number of other case studies performed at very small values of Ca (not reproduced here for the sake of brevity), for a wide range of the rest of the parameters, indeed support the fact that large deformation, such as in figure 5, can be obtained for those cases as well. Therefore, when analysing the flow dynamics in such cases, one can no longer assume that interface remains flat; rather, the flow dynamics has to be analysed by taking the deformation into account.

Figure 5.

(a) Plot (both analytical solution and numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for two different values of $m=2$ and 5, ${\it\varepsilon}=2$ , $q=2{\rm\pi}/4$ and $\mathit{Ca}=0.01$ . (b) Plot (both analytical solution and numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for two different values of $m=2$ and 5, ${\it\varepsilon}=3$ , $q=2{\rm\pi}/5$ and $\mathit{Ca}=0.005$ . The other parameters have the following values: $n=0.2$ , ${\it\lambda}=1$ , ${\it\beta}=1$ , ${\it\theta}=0$ and $\mathit{Cn}=0.02$ . This figure shows that, even for very low values of Ca, one can have a fairly large amount of deformation.

Figure 6 demonstrates the maximum deformation at steady state ( $h_{max}$ ) as a function of Ca; the other parameters are given in the caption. Here, we plot the numerical results along with the deformation until $O(\mathit{Ca})$ (i.e. $\mathit{Ca}h_{1}$ ) as well as $O(\mathit{Ca}^{2})$ (i.e. $\mathit{Ca}h_{1}+\mathit{Ca}^{2}h_{2}$ ). We note that for low values of Ca ( $\text{until}\approx 0.1$ ), the agreement between the numerical and analytical results is near perfect. As we move towards higher values of Ca, the solution until $O(\mathit{Ca}^{2})$ is seen to offer greater accuracy. Reasonably good agreement between $h_{max}$ from analytical (until $O(\mathit{Ca}^{2})$ ) and numerical results is observed for $\mathit{Ca}\sim 0.2{-}0.25$ , beyond which the error becomes increasingly large. It should be noted that, for all values of Ca, the analytical solutions underestimate the deformation, albeit this is only noticeable for $\mathit{Ca}>0.1$ .

Figure 6.

The maximum deformation at steady state ( $h_{max}$ ) versus Ca is plotted for ${\it\lambda}=1$ , $q=2{\rm\pi}/3$ , ${\it\theta}=0$ , ${\it\varepsilon}=1$ , $m=0.2$ and $n=0.2$ . We plot the numerical solutions (dots) along with analytical solutions until $O(\mathit{Ca})$ (dotted line) and $O(\mathit{Ca}^{2})$ (solid line) for comparison. Quite intuitively, for low Ca, the agreement is near perfect, while for higher Ca, the deformation until $O(\mathit{Ca}^{2})$ offers a better match with the numerical results.

Figure 7.

(a) Plot of the maximum deformation at steady state ( $h_{max}$ , both analytical and numerical) versus Ca, for $q=2{\rm\pi}/3$ and $2{\rm\pi}/4$ , with ${\it\lambda}=1$ and ${\it\varepsilon}=1$ . Good agreement is observed between the two solutions for low Ca. (b) Plot of the maximum deformation at steady state ( $h_{max}$ , only numerical) versus Ca, for ${\it\lambda}=0.1$ , 1 and 10, with $q=2{\rm\pi}/3$ and ${\it\varepsilon}=1$ . (c) Plot (both analytical solution and numerical results) of the maximum deformation at steady state ( $h_{max}$ ) versus Ca, for ${\it\varepsilon}=0.2$ and 2, with $q=2{\rm\pi}/3$ and ${\it\lambda}=1$ . The other parameters have the following values: $m=n=0.2$ , ${\it\beta}=1$ , ${\it\theta}=0$ and $\mathit{Cn}=0.02$ (except for $\mathit{Ca}=0.025$ , 0.05 and 0.075, where we have used $\mathit{Cn}=0.01$ ).

Figure 7(a) depicts the variation in the maximum deformation at steady state as a function of Ca (both analytical and numerical) for two different values of $q=2{\rm\pi}/3$ and $2{\rm\pi}/4$ ; the values of other relevant parameters are given in the caption. The agreement of the analytical results with the numerical ones is good for lower values of Ca ( ${\sim}0.2$ for $q=2{\rm\pi}/3$ and $\mathit{Ca}\sim 0.15$ for $q=2{\rm\pi}/4$ ). However, in the case of $q=2{\rm\pi}/4$ , the deformation becomes large ( ${\sim}0.7$ ) for higher capillary numbers, and does not remain in the limit of small deformation, as assumed for the analytical solutions. This figure shows that the magnitude of the deformation increases as the wavelength of modulation is increased, which is in tandem with the conclusions drawn from figure 3(c). As expected, this increase is more pronounced for higher values of Ca. We next move to figure 7(b), where we again plot the maximum deformation at steady state as a function of Ca, this time for three different values of ${\it\lambda}=0.1$ , 1 and 10; the values of other relevant parameters are given in the caption. This figure generalizes the results obtained from figure 3(b). An increased viscosity ratio increases the deformation for $m<1$ , and as expected this increment is more pronounced for higher values of the capillary number. We observe very large deformations ( $h_{max}\sim 0.6$ ) for ${\it\lambda}=10$ at high capillary numbers in this figure. Finally, figure 7(c) depicts the variation in the maximum steady-state deformation with Ca, for two values of permittivity ratio, ${\it\varepsilon}=0.2$ and 2; the values of other parameters are given in the caption. This figure is in agreement with figure 4, where we showed that an increment in the permittivity ratio leads to a decrease in the deformation. The same trend is also observed in the present figure. As expected, the change in the deformation is more pronounced for higher values of the capillary number. The analytical and numerical solutions show good agreement for Ca values up to ${\sim}0.2$ . The agreement is better for ${\it\varepsilon}=2$ , since the deformation is smaller in those cases. This figure (7 a–c) shows that the interface deformation can be augmented as well as suppressed without changing the external parameters, just by adjusting the fluid properties. It is important to note here that although our numerical simulations, in principle, should be valid for large values of $h$ , when $h$ is very close to 1, our formulation breaks down. It should be noted that here we have assumed the Smoluchowski slip velocity condition at the wall to be valid. The main condition for its validity is that the EDL should be thin compared with the channel height, a condition that is met in most of the microscale flows. However, we must also note the fact that in spite of being thin, the EDL has some finite non-zero thickness. Therefore, for very large values of $h$ (for very high deformation), the interface will come very close to the wall and thus can eventually end up inside the EDL. In such cases, the simple Smoluchowski condition will not be valid and one has to analyse the details of charge dynamics and fluid flow inside the EDL as well, starting from the more basic Poisson–Nernst–Planck equations. Although potentially interesting, in the present study we do not actively pursue the aforementioned detailed analysis of the EDL and stick to a simpler Smoluchowski model.

Figure 8 demonstrates the phase plot of the analytical solution for the maximum deformation, as a function of ${\it\varepsilon}$ and ${\it\lambda}$ ; the other parameters are given in the caption. This figure also shows the line of zero or negligible deformation on the ${\it\varepsilon}{-}{\it\lambda}$ plane. The variation in $h_{max}$ here contains two distinct regions, shown as regions I and II in the figure. In region I, the maximum deformation decreases as we increase ${\it\varepsilon}$ , whereas with increase in ${\it\lambda}$ it increases. This region corresponds to the variations already shown in figure 7(b,c). However, in region II, this trend reverses, since the deformation (maximum) increases with increase in ${\it\varepsilon}$ and decreases when ${\it\lambda}$ is increased. These two regions are separated by the line of negligible or zero deformation, where the deformation almost vanishes. Therefore, we conclude that for every choice of viscosity ratio, there is a certain permittivity ratio for which the deformation is completely suppressed. In fact, upon crossing the ‘curve of zero deformation’ the phase of the deformation shifts by 180°, which results in the cusp shown in figure 8. Our numerical results (although not shown here for sake of brevity) also show that for a particular choice of ${\it\lambda}$ , one can find a certain permittivity ratio for which the deformation almost vanishes.

Figure 8.

Phase plot of the maximum deformation at steady state (from analytical solutions) with the permittivity ( ${\it\varepsilon}$ ) on the $x$ axis and the viscosity ratio ( ${\it\lambda}$ ) on the $y$ axis. We have also shown the line of minimum $h_{max}$ , denoted as ‘curve of zero deformation’. On this line the deformation is negligible, even for higher values of Ca. The rest of the parameters have been chosen as follows: $m=n=0.2$ , ${\it\beta}=1$ , $\mathit{Ca}=0.2$ .

5.2. Transient deformation

In this section, we focus on the variation of the interface deformation with time. It is important to note here that although the interface changes shape with time, the fluid flow always remains quasisteady, as attributable to the approximation $\mathit{Oh}^{2}\gg 1$ . The transient evolution of the interface is governed by the equations (3.8), (3.9), (3.17) and (3.19). From these equations it is apparent that the relaxation time (time required to reach steady state) is solely determined by the terms $S_{1}$ and $S_{3}$ , which are functions of ${\it\lambda}$ and $q$ . To study the relaxation of the deformed interface, we present the influence of these parameters on the relaxation. Since one of our primary aims is to compare the analytical results for the transience with the numerical results, it is of paramount importance to rescale the time, since the characteristic values used for non-dimensionalizing time in the governing equations are different for the analytical and numerical solutions. To this end, we note that the ratio of the numerical and analytical time scales is Ca. Therefore, comparison between the transience of the two solutions will be possible if the analytical time value ( $t$ ) is multiplied by Ca. In this section all the analytical solutions have been presented by multiplying the corresponding time by Ca.

Figure 9.

Variation in the deformation ( $h$ ) with time at an axial location $x_{0}=2.2$ , for $\mathit{Ca}=0.1$ , 0.2 and 0.3, while the other parameters are as follows: ${\it\lambda}=1$ , $q=2{\rm\pi}/3$ , ${\it\varepsilon}=1$ , ${\it\beta}=1$ , ${\it\theta}=0$ , $m=0.2$ and $n=0.2$ . We have plotted both the analytical and numerical results; a reasonable agreement is observed between them. The deformation clearly shows monotonic increment in time, as predicted in § 3.

Figure 9 demonstrates the variation in the deformation at an arbitrarily chosen axial location, $x_{0}=2.2$ , with time for three different capillary numbers (0.1, 0.2 and 0.3); the other relevant parameters are given in the caption. We have plotted both analytical results (as given in (3.8)–(3.9) and (3.17)–(3.20)) and the numerical solutions for comparison. As in the previous results, the agreement is reasonably good for lower capillary numbers ( $\text{until}\sim 0.2$ ), while the error is relatively higher for $\mathit{Ca}=0.3$ . Both the analytical and numerical results depict similar monotonic variation in time, with roughly equal time scales of variation. The monotonic variation observed in the present figure can be explained from (3.9) and (3.19), which shows that the leading-order deformation, i.e. the most dominant component of the deformation, as well as one of the components from $O(\mathit{Ca}^{2})$ (the $c_{4}(t)$ term) varies monotonically with time. Therefore, the overall deformation also turns out to be monotonic in time.

Figure 10 depicts the variation in the normalized deformation (analytical as well as numerical) as a function of time, for three different ${\it\lambda}=0.1$ , 1 and 10 in (a) and for three different $q=2{\rm\pi}/3,2{\rm\pi}/4$ and $2{\rm\pi}/5$ in (b); the other parameters are given in the caption. We have normalized the deformation with the corresponding steady-state deformation, $h_{\infty }$ . Figure 10(a) is intended to highlight the effects of the viscosity ratio on the time taken by the deformed interface to reach the steady state. To this end, good agreement between numerical and analytical predictions is observed, since the times taken for both are practically equal in figure 10(a). The major point to note from this figure is that the time to reach steady state increases as the viscosity ratio is increased, and this increase in quite drastic for a change of ${\it\lambda}$ from 1 to 10. The augmented relaxation time can be attributed to the fact that higher viscosity will result in slower motion of both the fluids, which will eventually lead to a slower deformation rate. This is manifested as a larger relaxation time for interface deformation. In addition to this, increased ${\it\lambda}$ also leads to increased asymmetry in the flow for $m<1$ . This leads to larger deformation, which eventually causes the total deformation time to increase. It is interesting to note that both the magnitude of the deformation as well as the time taken to reach the steady state are larger for larger  ${\it\lambda}$ . On the other hand, with increase in the modulation wavelength (decrease in $q$ ), there is a significant increase in the time required to reach steady state. This can be attributed to the fact that as we increase the modulation wavelength, the deformation tends to increase, as shown in figures 3(c) and 7(a). Thus, with all other parameters (especially the viscosity) kept constant, the interface is naturally expected to take more time to reach steady state, simply because it has to deform more. In this case also the asymptotic theory agrees well with the numerical simulation results.

Figure 10.

The normalized deformation (normalized with the steady-state deformation, $h_{\infty }$ ) as a function of time is plotted at an axial location of $x_{0}=2.2$ (a) and $x_{0}=2.5$ (b): (a) for three different ${\it\lambda}=0.1$ , 1 and 10, with $\mathit{Ca}=0.1$ , $q=2{\rm\pi}/3$ and $m=n=0.2$ ; (b) for three different $q=2{\rm\pi}/3$ , $2{\rm\pi}/4$ and $2{\rm\pi}/5$ , with $\mathit{Ca}=0.05$ , ${\it\lambda}=1$ and $m=n=0.2$ . The other parameters were taken as ${\it\varepsilon}=1$ , ${\it\theta}=0$ , ${\it\beta}=1$ . We plot both analytical and numerical solutions and good agreement is observed between the two. This plot is intended to demonstrate the effects of viscosity and modulation wavelength on the time taken to reach the steady state. It is important to note that the analytical and numerical values have been normalized with their respective steady-state deformations.

Figure 11.

(a1,a2) Contours of the steady-state streamlines depicting the flow path, for $\mathit{Ca}=0.1$ , ${\it\lambda}=1$ , $m=n=0.2$ , $q=2{\rm\pi}/3$ , and the interface, as shown by the thick line. (b1,b2) Contours of the steady-state streamlines depicting the flow path, for $\mathit{Ca}=0.2$ , ${\it\lambda}=2$ , $m=n=0.2$ , $q=2{\rm\pi}/3$ , and the interface, as shown by the thick line. In both the figures, (1) is from the numerical solution, whereas (2) shows the analytical solution. Reasonably good agreement is observed between the analytical and numerical solutions in both the cases. The other parameters are: ${\it\beta}=1$ , ${\it\epsilon}=1$ and ${\it\theta}=0$ .

5.3. Flow dynamics

In this section we discuss the nature of the flow for the combination of parameters presented in the previous two sections. We start with figure 11(a,b), where the contours of the steady-state streamlines are shown as obtained from numerical (figure 11 a1,b1) and analytical (figure 11 a2,b2) results; the relevant parameters are given in the caption. We also show the location of the deformed interface (thick line) to differentiate the streamlines of the two fluids. Analytical solutions for the streamlines are given in (3.7) and (3.16), while the streamlines from the numerical solutions are calculated from the relation $\text{d}{\it\psi}=u\,\text{d}y-v\,\text{d}x$ , where it has been assumed that for the numerical solutions ${\it\psi}=0$ at the bottom wall. In both the cases the agreement between the numerical and analytical solutions is reasonable. Several interesting points can be noted from the present figure. First, recirculation rolls are present in both the fluids, which are attributable to the axially modulated potential on the wall. The axially invariant part of the potential contributes to the wavy but non-recirculating streamlines, leading to a net throughput. The waviness in these streamlines is, of course, an outcome of the presence of an axially modulated surface potential. A second and the most interesting point to note is the effect of the deformation on the shape of the streamlines. It can be observed that the deformation distorts all the streamlines, including the recirculation rolls, since they tend to be parallel to the interface at steady state. This distortion is, of course, more pronounced for the case of $\mathit{Ca}=0.2$ and ${\it\lambda}=2$ in figure 11(b1,b2), since the deformation is higher for this combination of parameters.

Figure 12 (values of relevant parameters are given in the caption) embodies an extreme case of the distortion of streamlines, as reported in the preceding figures. It can be observed that the deformation is very large (maximum deformation ${\sim}0.7$ ), which distorts all the streamlines extremely to make them parallel to the interface at steady state. Here, we have only shown the numerically obtained streamlines, since the approximate analytical solutions are not expected to hold for such high deformation magnitudes.

Figure 12.

Contours of the steady-state streamlines (only numerical solutions) for large interface deformation, with $\mathit{Ca}=0.3$ , ${\it\lambda}=1$ , $m=n=0.2$ , ${\it\varepsilon}=1$ , ${\it\beta}=1$ , ${\it\theta}=0$ and $q=2{\rm\pi}/4$ . The location of the interface has also been plotted (thick line). This figure demonstrates an extreme case of deformation and streamline distortion.

Figure 13.

(a,b) Contours of the steady-state streamlines for $\mathit{Ca}=0.1$ , ${\it\lambda}=1$ , $q=2{\rm\pi}/3$ , ${\it\varepsilon}=1$ , ${\it\beta}=1$ , ${\it\theta}=0$ , $m=-0.2$ and $n=0$ . The location of the interface has been shown with a thick line. Panel (a) demonstrates the results obtained from numerical solutions, whereas (b) is from the analytical solution. (c) Plot of the steady-state deformation for two different values of $m=0.2$ (curve with lower amplitude) and $-0.2$ (curve with larger amplitude), while the other parameters are the same as in (a) and (b). We have plotted both numerical (circles) and analytical (lines) solutions, and the agreement between them is reasonably good.

We next investigate a special case in figure 13(a–c), where we depict the contours of the streamlines for $m=-0.2$ and $n=0$ (both analytical (figure 13 a) and numerical (figure 13 b)), which indicates that the potentials at the top and the bottom walls are out of phase (values of other relevant parameters are given in the caption). Good agreement is observed between the analytical and numerical results. It is interesting to note that in the present case two recirculation rolls appear in the top fluid, whereas for cases with $m>0$ , reported in the preceding figures, only one roll exists in the top fluid. This can be attributed to the fact that for $m<0$ the fluid on top tries to move in the opposite direction to the bottom fluid, whereas the bottom fluid tries to drag the top fluid with itself. Since $m=0.2$ , the slip velocity in the top fluid is weak, and hence the part of this fluid residing close to the interface feels more drag due to the bottom fluid (note that ${\it\lambda}=1$ ) compared with the electro-osmotic forcing due to Smoluchowski slip. Therefore, this part of the fluid tends to move along the direction of the bottom fluid, whereas the part near the top wall moves in the opposite direction. This mismatch in the direction of movement gives rise to multiple rolls in the fluid on top. Finally, we note that the deformation of the interface (at steady state) for $m=-0.2$ has a greater magnitude than that for $m=0.2$ (i.e. $m>0$ ), as shown in figure 13(c), where the steady-state deformation has been plotted keeping the parameter values the same as in (a) and (b), only altering $m$ between $-0.2$ and 0.2. This can be explained very easily based on the argument of asymmetry, mentioned in the explanation of figure 3(b). Clearly, $m=-0.2$ signifies more asymmetry, and since the deformation is primarily caused here by asymmetry, it leads to a higher amount of deformation.

Following figure 13, we now attempt to show the deformation and the resulting flow structures for the intermediate phase shift angles ( ${\it\theta}$ ) as depicted in figure 14(a–c) (only numerical results). Figure 14(a) demonstrates the variation in the interface deformation $(h_{\infty })$ at steady state, for five different choices of ${\it\theta}=0$ , ${\rm\pi}/4$ , ${\rm\pi}/2$ , $3{\rm\pi}/4$ and ${\rm\pi}$ ; the values of the rest of the parameters are given in the caption. We note that as the phase shift angle increases, the deformation increases, reaching a maximum at ${\it\theta}={\rm\pi}$ . This can again be explained from the arguments of asymmetry, since increasing ${\it\theta}$ leads to increased asymmetry between the two halves of the channel, which in turn causes the interface to deform more. Finally, in figure 14(b,c) we show the streamlines for two choices of ${\it\theta}={\rm\pi}/2$ and ${\rm\pi}$ . We can see that the only notable effect of ${\it\theta}$ here is that the position of the recirculation rolls in the upper fluid (fluid 1) changes and tends to shift to the left as ${\it\theta}$ is increased. At the same time, the roll in fluid 2 tends to get distorted as ${\it\theta}$ is increased.

Figure 14.

(a) Plot (only numerical results) of the steady-state deformation $(h_{\infty })$ versus $x$ , for five different values of the phase shift angle, ${\it\theta}=0$ , ${\rm\pi}/4$ , ${\rm\pi}/2$ , $3{\rm\pi}/4$ and ${\rm\pi}$ . (b,c) Contours of the steady-state streamlines (only numerical results) depicting the flow path, for two different values of the phase shift angle, (b) ${\it\theta}={\rm\pi}/2$ and (c) ${\it\theta}={\rm\pi}$ . The other parameters have the following values: $m=2$ , $n=0.2$ , ${\it\lambda}=0.5$ , $q=2{\rm\pi}/3$ , ${\it\beta}=1$ , ${\it\varepsilon}=1$ , $\mathit{Ca}=0.1$ and $\mathit{Cn}=0.02$ . The interface is represented by the red line.

We acknowledge that more parametric studies, apart from those presented here, can be executed depicting the resulting streamline and deformation patterns. However, we strongly feel that the parametric studies reported herein are sufficient to highlight the essential features of flow dynamics and deformation patterns, within the assumptions of the present framework. Therefore, we restrain ourselves from performing any more such studies and conclude the present analysis.