Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium

Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium //// Advanced Search home journals subjects indexing guidelines submit resources about Journal Menu About this Journal Abstracting and Indexing Aims and Scope Article Processing Charges Articles in Press Bibliographic Information Contact Information Editorial Board Editorial Workflow Free eTOC Alerts Table of Contents Abstract Full-Text PDF Full-Text HTML Full-Text ePUB Linked References How to Cite this Article ISRN Applied Mathematics Volume 2011 (2011), Article ID 498718, 35 pages doi:10.5402/2011/498718 Review Article <h2>Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium</h2> M. F. El-Sayed , 1,2 M. H. M. Moussa , 1 A. A. A. Hassan , 1 and N. M. Hafez 1 1 Department of Mathematics, Faculty of Education, Ain Shams University, Heliopolis (Roxy), Cairo 11757, Egypt 2 Department of Mathematics, College of Science, Qassim University, P.O. Box 6644, Buraidah 51452, Saudi Arabia Received 23 March 2011; Accepted 4 May 2011 Academic Editor: R. Cardoso Copyright © 2011 M. F. El-Sayed et al. This is an open access article distributed under the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract The effect of an applied electric field on the stability of the interface between two thin viscous leaky dielectric fluid films in porous medium is analyzed in the long-wave limit. A systematic asymptotic expansion is employed to derive coupled nonlinear evolution equations of the interface and interfacial free charge distribution. The linearized stability of these equations is determined and the effects of various parameters are examined in detail. For perfect-perfect dielectrics, the various parameters affect only for small wavenumber values. For dielectrics, the various parameters affect only for small wavenumber values. For effect for small wavenumbers, and a stabilizing effect afterwards, and for high wavenumber values for the other physical parameters, new regions of stability or instability appear. For leaky-leaky dielectrics, the conductivity of upper fluid has a destabilizing effect for small or high wavenumbers, while it has a dual role on the stability of the system in a wavenumber range between them. The effects of all other physical parameters behave in the same manner as in the case of perfect-leaky dielectrics, except that in the later case, the stability or instability regions occur more faster than the corresponding case of leaky-leaky dielectrics. 1. Introduction The effect of electric fields on the stability and dynamics of fluid-fluid interfaces has been an area of extensive research, beginning from the classic works of Taylor and McEwan [ 1 ] and Melcher and Smith [ 2 ]. These works and subsequent studies, see, for example, the reviews of Saville [ 3 ] and Griffiths [ 4 ], have amply demonstrated the role of electrical stresses on fluid interfaces and the associated electrohydrodynamic instabilities in such systems. One of the basic problems here is to understand the stability of the interface between two fluid layers bounded on the top and bottom by rigid plates, and this has been the subject of many previous studies. Mohamed et al. [ 5 ] concentrated on two superposed viscous fluids in a channel subjected to a normal electric field, where the upper fluid is highly conducting, while the lower fluid is dielectric, and they performed the long-wave linear stability analysis, [ 6 , 7 ], and showed that the electric field always has a destabilizing effect on the flow. Abdella and Rasmussen [ 8 ] studied Couette flow of two viscous fluids with different viscosities, densities, conductivities and permittivities, in an unbounded domain subjected to a normal electric field. They studied, following Melcher [ 9 ], two special cases in detail: the electrohydrodynamic free-charge configuration (EH-If) and the electrohydrodynamic polarization charge configuration (EH-If). These studies have largely considered systems in which gravitational effects are important, and therefore, a critical applied voltage is required to cause the instability, very long waves are stabilized by interfacial tension, and waves of intermediate lengths become unstable. These earlier studies have focused on how the critical voltage required for instability is affected by the nature of the fluids, namely. whether they are perfect dielectrics, or whether they are leaky dielectrics in which there is the possibility of free charge conduction in the fluids, and also the possibility of accumulation and redistribution of charges on the interface between two fluids. Recently, there has been a renewed interest in this area, in part due to the relevance of such phenomena in the formation of well-controlled patterns using the application of electric fields to thin liquid films [ 10 – 12 ], which has demonstrated that the application of an external electric field to polymer-air or polymer-polymer interfaces enhances the spontaneous fluctuations at the interface leading to an instability. However, for leaky dielectrics with surface charges at the interface and two fluids that are not perfectly conducting, we must bear in mind that there is an electrical tangential shear stress at the interface, induced by the electric field, and hence, it changes the stability of the two-fluid layer system as previously investigated by Ozen et al. [ 13 ]. Wu and Chou [ 14 ] performed linear stability analysis of a leaky dielectric viscoelastic fluid whose constitutive behavior was described using the Jeffrey's model [ 15 ]. The surface instability of a Newtonian fluid (modeled as both leaky and perfect dielectrics) under the effect of electric field is now well understood [ 16 ]. For recent developments of this topic, see the investigations of Papageorgiou and Petropoulos [ 17 ], Shankar and Sharma [ 18 ], Craster and Matar [ 19 ], Li et al. [ 20 ], Tomar et al. [ 21 ], and Supeene et al. [ 22 ]. Theoretical studies that have considered the linear stability characteristics of thin fluid films subjected to electric fields were restricted to the following configurations: (i) the interface between a perfect dielectric liquid and air [ 10 ], (ii) the interface between two perfect dielectric liquids [ 11 ], (iii) the interface between a dielectric fluid of finite thickness and a conducting fluid of much larger thickness [ 23 ], and (iv) the interface between a leaky dielectric liquid and air [ 24 ]. Except the study of Pease III and Russel [ 24 ], all these studies have focused only on perfect dielectric systems. The presence of conductivity in one or both of the liquids could have a significant impact on the length scales and growth rates. It should be noted that in all the above-mentioned studies, the medium has been considered to be nonporous. Porous media theories, on one hand, play an important role in many branches of engineering, including material science, petroleum industry, chemical engineering, and soil mechanics as well as biomechanics. The flow through porous media has gained considerable interest in recent years, particularly among geophysical fluid dynamicists. It is well known that in Darcy's law, which relates the pressure gradient, the bulk viscous resistance, and the gravitational force in a porous medium, the usual viscous term in the equation of motion is replaced by the resistive force − ( 𝜇 / 𝑘 1 ) 𝐯 , where 𝜇 is the fluid viscosity, 𝑘 1 is the medium permeability, and 𝐯 is the Darcian velocity of the fluid. Much of the recent studies on this topic are given by Ingham and Pop [ 25 ], Vafai [ 26 ], Del Rio and Whitaker [ 27 ], Pop and Ingham [ 28 ], and Nield and Bejan [ 29 ]. On the other hand, electrohydrodynamic instability studies for flows in porous media has attracted little attention in the scientific literature [ 30 – 36 ] despsite their applications in various diverse fields with great interest. Thus, there is a growing need for original research in the updated electrohydrodynamic phenomena which have some physical and engineering applications. In this study, we consider the most general case of the effect of an external applied electric field on the stability of the interface between two thin leaky dielectric fluids with arbitrary viscosities and conductivities. We first use a systematic long-wave asymptotic analysis to derive the nonlinear evolution equations for the interface position and interfacial charge distribution, and we subsequently study the linearized stability of the nonlinear differential equations. This paper is organized in the following manner: Section 2 discusses the relevant governing equations and boundary conditions, and in Section 3 , we nondimensionalize these equations and conditions. In Sections 4 and 5 , we outline the long-wave asymptotic analysis used to derive the nonlinear evolution equations for the interface position and charge, and in Section 6 , we develop the linear stability analysis of the nonlinear equations, and discuss the important representative studies results obtained from the general case of two leaky dielectric interface and the limiting cases of both perfect-leaky dielectric, and perfect-perfect dielectric interfaces. Finally, the salient conclusions of the present study are discussed in Section 7 . 2. Problem Formulation and Governing Equations The system of interest consists of two leaky dielectric fluids in porous medium of arbitrary viscosities occupying the regions − 𝐻 < 𝑦 < 0 (fluid 2) and 0 < 𝑦 < 𝛽 𝐻 (fluid 1) in the initial unperturbed state, see Figure 1 , where 𝛽 is the ratio of the thicknesses of top and bottom fluids. The perturbed interface between the two fluids is denoted by 𝑦 = ℎ ( 𝑥 ) . The two fluids are stationary in the initial state, with viscosities 𝜇 𝑖 , dielectric constants 𝜖 𝑖 , conductivities 𝜎 𝑖 , porosity of porous medium 𝜀 and medium permeability 𝑘 1 , where 𝑖 = 1 , 2 . Fluid 2 is bounded at the bottom 𝑦 = − 𝐻 by a rigid plate which is maintained at an electric potential 𝜓 = 𝜓 𝑏 , while fluid 1 is bounded at the top 𝑦 = 𝛽 𝐻 by a rigid plate maintained at an electric potential 𝜓 = 0 . In the ensuing analysis, we assume that the material properties of the fluid such as viscosities 𝜇 𝑖 , dielectric constants 𝜖 𝑖 , conductivities 𝜎 𝑖 , porosity of porous medium 𝜀 , and medium permeability 𝑘 1 are constants and independent of spatial position. Following the leaky dielectric model formulation of Saville [ 3 ], we assume that electroneutrality is valid in the bulk, while free charge is assumed to accumulate at the fluids interface. We also neglect the diffusion of free charge within the interface. Figure 1: Schematic diagram showing the configuration and coordinate system. For leaky dielectric fluids of constant conductivities and with zero net charge in the bulk, the following governing equations are appropriate for the electric field 𝐄 in the two fluids 1 and 2 [ 3 ] ∇ ⋅ 𝐄 𝑖 = 0 , 𝑖 = 1 , 2 . ( 2 . 1 ) Since the electric fields are irrotational, 𝐄 𝑖 = − ∇ 𝜓 𝑖 , where 𝜓 𝑖 is the electric potential: in the fluid 𝑖 . Substituting this into ( 2.1 ) gives the following Laplace's equation for the electric potential: 𝜓 𝑖 ∇ 2 𝜓 𝑖 = 0 . ( 2 . 2 ) These governing equations are supplemented by the following boundary conditions. (1) The normal component of the electric field, at the interface 𝑦 = ℎ ( 𝑥 ) , satisfies 𝜖 2 𝜖 0  ∇ 𝜓 2  ⋅ 𝐧 − 𝜖 1 𝜖 0  ∇ 𝜓 1  ⋅ 𝐧 = 𝑞 ( 𝑥 , 𝑡 ) a t 𝑦 = ℎ ( 𝑥 ) , ( 2 . 3 ) where 𝐧 is the unit normal to the interface at 𝑦 = ℎ ( 𝑥 ) (see Figure 1 ), 𝜖 𝑖 ( 𝑖 = 1 , 2 ) is the dielectric constant in fluid 𝑖 , 𝜖 0 is the permittivity of free space, and 𝑞 ( 𝑥 , 𝑡 ) is the surface charge density of free charges at the interface. (2) The continuity of the tangential component of the electric field at the interface 𝑦 = ℎ ( 𝑥 ) translates to the continuity of the electric potentials; that is, 𝜓 1 = 𝜓 2 a t 𝑦 = ℎ ( 𝑥 ) . ( 2 . 4 ) (3) The electric potentials satisfy the following conditions at the rigid boundaries: 𝜓 2 = 𝜓 𝑏 𝜓 a t 𝑦 = − 𝐻 , 1 = 0 a t 𝑦 = 𝛽 𝐻 . ( 2 . 5 ) We next turn to the equations governing the motion of the two fluids. Owing to the relatively small thicknesses of the fluids, we ignore inertial effects in both fluids, and hence, the governing equations are the Stokes equations for continuity and momentum balance ∇ ⋅ 𝐕 𝑖 = 0 , ∇ ⋅ 𝐓 𝑖 = 0 , ( 2 . 6 ) where 𝐕 𝑖 and 𝐓 𝑖 are the velocity field and the total stress tensor, respectively, in fluid 𝑖 . Both the fluids are assumed to be irrotational, then the fluid velocity 𝐕 𝑖 can be derived from a scalar velocity potential 𝜑 𝑖 such that 𝐕 𝑖 = − ∇ 𝜑 𝑖 . We have neglected the effects of gravity on the length scales of interest here. In addition, the effects of van der Waals dispersion forces are negligible for the films considered in the experimental studies [ 11 ]. The total stress tensor 𝐓 is given by a sum of isotropic pressure, deviatoric viscous stresses for the Newtonian fluid, the electrical Maxwell stress tensor, and a Darcy's law term describing the isotropic porous medium 𝐓 𝑖 = − 𝑝 𝑖 𝐈 + 𝜇 𝑖  ∇ 𝐕 𝑖 + ∇ 𝐕 𝑇 𝑖  + 𝐦 𝑖 + 𝜇 𝑖 𝑘 1 𝜑 𝑖 𝐈 , 𝑖 = 1 , 2 , ( 2 . 7 ) where 𝑝 𝑖 is the pressure in fluid 𝑖 , 𝐈 is the identity tensor, and the Maxwell stress tensor 𝐦 𝑖 is given by Saville [ 3 ] 𝐦 𝑖 = 𝜖 𝑖 𝜖 0  𝐄 𝑖 𝐄 𝑖 − 1 2  𝐄 𝑖 ⋅ 𝐄 𝑖  𝐈  . ( 2 . 8 ) The divergence of the Maxwell stress tensor ∇ ⋅ 𝐦 𝑖 = 0 , because the bulk of the fluid is free of net charge, and the dielectric constants are independent of spatial position in the two fluids; that is, ∇ ⋅ 𝐦 𝑖  𝜖 = ∇ ⋅ 𝑖 𝜖 0  𝐄 𝑖 𝐄 𝑖 − 1 2  𝐄 𝑖 ⋅ 𝐄 𝑖  𝐈 𝜖   = − 0 𝐄 𝑖 2 ⋅ 𝐄 𝑖 ∇ 𝜖 𝑖 + 𝜌 𝑓 𝐄 𝑖 = 0 , ( 2 . 9 ) with 𝜌 𝑓 is the bulk free charge. Thus, the Maxwell stress tensor will not appear in the momentum balance but will affect the flow only through the conditions at the interface. The governing momentum equations in the two fluids, therefore, become 𝜌 𝜀 𝐷 𝐕 𝑖 𝐷 𝑡 = − ∇ 𝑝 𝑖 + ∇ ⋅ 𝐦 𝑖 + 𝜇 𝑖 𝜀 ∇ 2 𝐕 𝑖 − 𝜇 𝑖 𝑘 1 𝐕 𝑖 , 𝑖 = 1 , 2 . ( 2 . 1 0 ) Since there is no time scale, then 𝐷 𝐯 𝑖 / 𝐷 𝑡 = 0 , and we get from the above two equations that ∇ 𝑝 𝑖 − 𝜇 𝑖 𝜀 ∇ 2 𝐕 𝑖 + 𝜇 𝑖 𝑘 1 𝐕 𝑖 = 0 , 𝑖 = 1 , 2 . ( 2 . 1 1 ) The fluid velocities satisfy no-slip and no-penetration conditions at the top and bottom plates 𝐕 1 ( 𝑦 = 𝛽 𝐻 ) = 0 , 𝐕 2 ( 𝑦 = − 𝐻 ) = 0 . ( 2 . 1 2 ) At the interface 𝑦 = ℎ ( 𝑥 ) between the two fluids, continuity of velocities and stresses apply ( 𝐕 ⋅ 𝐧 ) 1 = ( 𝐕 ⋅ 𝐧 ) 2 , ( ( 2 . 1 3 ) 𝐕 ⋅ 𝐭 ) 1 = ( 𝐕 ⋅ 𝐭 ) 2 , ( 2 . 1 4 ) ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 2 = ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 1 + 𝛾 𝜅 , ( 2 . 1 5 ) ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 2 = ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 1 , ( 2 . 1 6 ) where 𝛾 is the interfacial tension between the two fluids, 𝐭 is the unit tangent to the interface (see Figure 1 ), and 𝜅 = 𝜕 2 ℎ ( 𝑥 ) / 𝜕 𝑥 2 is the mean curvature of the interface. We restrict our attention to two-dimensional systems which are invariant in the 𝑧 direction and denote the velocity components in the 𝑥 and 𝑦 directions by 𝑢 and 𝜐 , respectively; that is, 𝐕 𝑖 = ( 𝑢 𝑖 , 𝜐 𝑖 ) . Upon substituting the expression for the Maxwell stress tensor ( 2.8 ) in ( 2.15 ) and ( 2.16 ), respectively, we get ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 𝑖 = − 𝑝 𝑖 + 2 𝜇 𝑖 𝜕 𝜐 𝑖 𝜕 𝑦 + 𝜖 𝑖 𝜖 0  𝐸 𝑁 2 𝑖 − 1 2 𝐸 2 𝑖  + 𝜇 𝑖 𝑘 1 𝜑 𝑖 , ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 𝑖 = 𝜇 𝑖  𝜕 𝑢 𝑖 + 𝜕 𝑦 𝜕 𝜐 𝑖  𝜕 𝑥 + 𝜖 𝑖 𝜖 0 𝐸 𝑇 𝑖 𝐸 𝑁 𝑖 , ( 2 . 1 7 ) where 𝐸 𝑇 𝑖 and 𝐸 𝑁 𝑖 are the tangential and normal components of the electric field 𝐄 𝑖 in the fluid 𝑖 , respectively. Hence, we obtain the following conditions to be applied at the interface 𝑦 = ℎ ( 𝑥 ) :  𝑝 1 − 2 𝜇 1 𝜕 𝜐 1  −  𝑝 𝜕 𝑦 2 − 2 𝜇 2 𝜕 𝜐 2  + 1 𝜕 𝑦 𝑘 1  𝜇 2 𝜑 2 − 𝜇 1 𝜑 1  𝜖 = 𝛾 𝜅 + 1 𝜖 0 2   ∇ 𝜓 1  ⋅ 𝐧 2 −  ∇ 𝜓 1  ⋅ 𝐭 2  − 𝜖 2 𝜖 0 2   ∇ 𝜓 2  ⋅ 𝐧 2 −  ∇ 𝜓 2  ⋅ 𝐭 2  , ( 2 . 1 8 ) for the normal stress continuity, and 𝜇 2  𝜕 𝑢 2 + 𝜕 𝑦 𝜕 𝜐 2  𝜕 𝑥 − 𝜇 1  𝜕 𝑢 1 + 𝜕 𝑦 𝜕 𝜐 1   𝜕 𝑥 = − ∇ 𝜓 1  ⋅ 𝐭 𝑞 ( 𝑥 , 𝑡 ) , ( 2 . 1 9 ) for the tangential stress continuity. In the last equation, we have used the normal electric field continuity condition, ( 2.3 ), to simplify the right-hand side. The kinematic condition at the interface prescribes the evolution of the interface position ℎ ( 𝑥 , 𝑡 ) , 𝜐 1 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜐 2 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 ℎ , ( 2 . 2 0 ) where ∇ 𝑠 is the gradient operator along the interface 𝑦 = ℎ ( 𝑥 ) . Finally, the interfacial charge is governed by a conservation equation 𝜀 𝜕 𝑞 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 𝜎 𝑞 − 𝑞 𝐧 ⋅ ( 𝐧 ⋅ ∇ ) 𝐕 = 𝜀   2 𝐄 2  −  𝜎 ⋅ 𝐧 1 𝐄 1 ⋅ 𝐧   , ( 2 . 2 1 ) where the terms on the left-hand side represent, respectively, the accumulation, convection, and variation of the charge due to dilation of the interface, while the right side represents the migration of charge to or from the interface due to ion conduction in the bulk [ 3 ]. 3. Nondimensional Forms It is useful at this point to nondimensionalize the governing equations and boundary conditions by setting 𝜓  𝑖 = 𝜓 𝑖 𝜓 𝑏 , 𝑝  𝑖 = 𝑝 𝑖 𝐻 2 𝜖 0 𝜓 2 𝑏 , 𝑇  𝑖 = 𝑇 𝑖 𝐻 2 𝜖 0 𝜓 2 𝑏 , 𝑥  = 𝑥 𝐻 , 𝑦  = 𝑦 𝐻 , 𝐄  𝑖 = 𝐄 𝑖 𝐻 𝜓 𝑏 , 𝐕  𝑖 = 𝐕 𝑖 𝜇 2 𝐻 𝜖 0 𝜓 2 𝑏 , 𝑡  = 𝑡 𝜖 0 𝜓 2 𝑏 𝜇 2 𝐻 2 , 𝜑  𝑖 = 𝜇 2 𝜑 𝑖 𝜖 0 𝜓 2 𝑏 , 𝑞  = 𝑞 𝐻 𝜖 0 𝜓 𝑏 , 𝜅  = 𝜅 𝐻 , 𝑘  1 = 𝑘 1 𝐻 2 , ℎ  = ℎ 𝐻 . ( 3 . 1 ) Upon using these scales to nondimensionalize the above governing equations and boundary conditions, we end up with the following nondimensional set of equations. Without loss of clarity, and for the sake of brevity, we represent nondimensional variables with the same notation in the ensuing discussion by dropping dashes. The nondimensional governing equations for the electric potentials 𝜓 𝑖 are ∇ 2 𝜓 𝑖 = 0 , 𝑖 = 1 , 2 , ( 3 . 2 ) with the following boundary conditions at the interface 𝑦 = ℎ ( 𝑥 ) ,  𝜖 2 ∇ 𝜓 2  −  𝜖 ⋅ 𝐧 1 ∇ 𝜓 1  ⋅ 𝐧 = 𝑞 , 𝜓 1 = 𝜓 2 , ( 3 . 3 ) and the following boundary conditions at the top and bottom boundaries 𝜓 1 ( 𝑦 = 𝛽 ) = 0 , 𝜓 2 ( 𝑦 = − 1 ) = 1 . ( 3 . 4 ) Similarly, the nondimensional equations governing the fluid motion are ∇ ⋅ 𝐕 1 = 0 , ∇ ⋅ 𝐕 2 = 0 , ( 3 . 5 ) ∇ 𝑝 1 − 𝜇 𝑟 𝜀 ∇ 2 𝐕 1 + 𝜇 𝑟 𝑘 1 𝐕 1 = 0 , ∇ 𝑝 2 − 1 𝜀 ∇ 2 𝐕 2 + 1 𝑘 1 𝐕 2 = 0 . ( 3 . 6 ) with 𝜇 𝑟 = 𝜇 1 / 𝜇 2 being the ratio of viscosities of the two fluids. The nondimensional normal and tangential stress continuity conditions at the interface become  𝑝 1 − 2 𝜇 𝑟 𝜕 𝜐 1  −  𝑝 𝜕 𝑦 2 − 2 𝜕 𝜐 2  + 1 𝜕 𝑦 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝜖 𝛾 𝜅 + 1 2   ∇ 𝜓 1  ⋅ 𝐧 2 −  ∇ 𝜓 1  ⋅ 𝐭 2  − 𝜖 2 2   ∇ 𝜓 2  ⋅ 𝐧 2 −  ∇ 𝜓 2  ⋅ 𝐭 2  ,  ( 3 . 7 ) 𝜕 𝑢 2 + 𝜕 𝑦 𝜕 𝜐 2  𝜕 𝑥 − 𝜇 𝑟  𝜕 𝑢 1 + 𝜕 𝑦 𝜕 𝜐 1   𝜕 𝑥 = − ∇ 𝜓 1  ⋅ 𝐭 𝑞 ( 𝑥 , 𝑡 ) , ( 3 . 8 ) where 𝛾 = 𝛾 𝐻 / 𝜖 0 𝜓 2 𝑏 is the nondimensional interfacial tension. The boundary conditions for the velocities at the top and bottom plates become 𝑢 1 𝜐 ( 𝑦 = 𝛽 ) = 0 , 1 𝑢 ( 𝑦 = 𝛽 ) = 0 , 2 ( 𝑦 = − 1 ) = 0 , 𝜐 2 ( 𝑦 = − 1 ) = 0 . ( 3 . 9 ) The nondimensional kinematic condition at the interface is 𝜐 1 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜐 2 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 ℎ , ( 3 . 1 0 ) and the nondimensional charge conservation equation at the interface is 𝜀 𝜕 𝑞 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 𝑆 𝑞 − 𝑞 𝐧 ⋅ ( 𝐧 ⋅ ∇ ) 𝐕 = 𝜀   1 ∇ 𝜓 1  −  𝑆 ⋅ 𝐧 2 ∇ 𝜓 2 ⋅ 𝐧   , ( 3 . 1 1 ) where 𝑆 𝑖 = 𝜎 𝑖 𝜇 2 𝐻 2 / 𝜖 2 0 𝜓 2 𝑏 , 𝑖 = 1 , 2 are the nondimensional conductivities in the two fluids. This completes the specification of the governing equations and boundary conditions, which are highly coupled. Due to the negligible effect of gravity at length scales of interest here, the above system of equations undergoes a long-wave instability. We now carry out a long-wave asymptotic analysis to make the above system of equations tractable, and thereby derive coupled nonlinear evolution equations for the interface position ℎ ( 𝑥 , 𝑡 ) and charge 𝑞 ( 𝑥 , 𝑡 ) . While the main focus of this paper is to analyze the stability of the linearized equations, it is nonetheless useful to first derive the nonlinear evolution equations, since these equations can be used in future studies to understand (by numerical simulations) the nonlinear evolution processes that occur after the linear instability. 4. Long-Wave Asymptotic Analysis In the long-wave limit, the wavelength 𝐿 of the fastest growing modes is much larger than the transverse length scale 𝐻 in the system, and it is useful to define a small parameter 𝛿 = 𝐻 / 𝐿 ≪ 1 . The lateral length scale 𝐿 is determined self-consistently in the following analysis to be  𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 , and this is further estimated below to be much larger than 𝐻 . Similarly, a slow time scale is necessary to describe the dynamics of the interface motion at such large length scales, and this is introduced a little later in ( 4.15 ). In the limit 𝐻 ≪ 𝐿 , the derivatives in the 𝑥 direction should be scaled with 𝐿 . To this end, we define the slowly varying scale 𝜒 in the following manner: 𝜕 𝜕 𝜕 𝑥 = 𝛿 𝜕 𝜒 , ( 4 . 1 ) and 𝜕 / 𝜕 𝜒 ∼ 𝑂 ( 1 ) . When we apply the above scalings, ( 4.1 ), to the Laplace equation, ( 3.2 ), for the electric potential 𝜓 𝑖 , it simplifies in the limit 𝛿 ≪ 1 to 𝜕 2 𝜓 𝑖 𝜕 𝑦 2 = 0 , 𝑖 = 1 , 2 . ( 4 . 2 ) The continuity condition for the normal component of the electric field at the interface 𝑦 = ℎ ( 𝜒 ) , ( 3.3 ), is similarly simplified in the long-wave limit as 𝜖 2 𝜕 𝜓 2 𝜕 𝑦 − 𝜖 1 𝜕 𝜓 1 𝜕 𝑦 = 𝑞 , ( 4 . 3 ) while the other interface condition (second equation in ( 3.3 )) and the boundary conditions, ( 3.4 ), remain unchanged in the long-wave limit. We now turn to the simplification of the momentum equations ( 3.6 ) for the fluid motion in the long-wave limit. It is useful to define the variable 𝜇 𝑟 , 𝑖 such that 𝜇 𝑟 , 𝑖 = 𝜇 𝑟 for 𝑖 = 1 and 𝜇 𝑟 , 𝑖 = 1 for 𝑖 = 2 . The 𝑥 -momentum equation can be simplified in the long-wave limit as 𝛿 𝜕 𝑝 𝑖 − 𝜇 𝜕 𝜒 𝑟 , 𝑖 𝜀 𝜕 2 𝑢 𝑖 𝜕 𝑦 2 + 𝜇 𝑟 , 𝑖 𝑘 1 𝑢 𝑖 = 0 . ( 4 . 4 ) This suggests that 𝑢 𝑖 ∼ 𝑂 ( 𝛿 ) 𝑝 𝑖 . In order to make this explicit and to make ordering of various quantities simpler, we represent the pressure 𝑝 𝑖 and the 𝑥 -component velocity 𝑢 𝑖 in the following manner: 𝑝 𝑖 = 𝑝 𝑖 ( 0 ) , 𝑢 𝑖 = 𝛿 𝑢 𝑖 ( 0 ) . ( 4 . 5 ) The above variables are the leading order quantities in an asymptotic expansion in 𝛿 , and we will be concerned only with the leading order variables in this paper. The continuity equation in both fluids ( 3.5 ) becomes, upon using 𝜕 / 𝜕 𝑥 ∼ 𝛿 𝜕 / 𝜕 𝜒 and using the above expansion for 𝑢 𝑖 , 𝛿 2 𝜕 𝑢 𝑖 ( 0 ) + 𝜕 𝜒 𝜕 𝜐 𝑖 𝜕 𝑦 = 0 , ( 4 . 6 ) which suggests the following expansion for 𝜐 𝑖 : 𝜐 𝑖 = 𝛿 2 𝜐 𝑖 ( 0 ) . ( 4 . 7 ) Upon using this expansion, the nondimensional 𝑦 component of the momentum equation yields to leading order in 𝛿 , 𝜕 𝑝 𝑖 ( 0 ) / 𝜕 𝑦 = 0 , implying that the pressure is constant in both films across the 𝑦 direction, and so 𝑝 𝑖 = 𝑝 𝑖 ( 𝜒 , 𝑡 ) . The simplified 𝑥 -momentum equation, therefore, is given by 𝑑 𝑝 𝑖 ( 0 ) − 𝜇 𝑑 𝜒 𝑟 , 𝑖 𝜀 𝜕 2 𝑢 𝑖 ( 0 ) 𝜕 𝑦 2 + 𝜇 𝑟 , 𝑖 𝑘 1 𝑢 𝑖 ( 0 ) = 0 . ( 4 . 8 ) The normal stress condition at the interface, ( 3.7 ), simplifies to give the following equation in the long-wave limit: 𝑝 1 ( 0 ) − 𝑝 2 ( 0 ) + 1 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝛾 𝛿 2 𝜕 2 ℎ 𝜕 𝜒 2 + 𝜖 1 2  𝜕 𝜓 1  𝜕 𝑦 2 − 𝜖 2 2  𝜕 𝜓 2  𝜕 𝑦 2 . ( 4 . 9 ) In order for the interfacial tension to be of the same order as the other terms in the above equation, we require 𝛾 𝛿 2 ∼ 𝑂 ( 1 ) , where 𝛿 = 𝐻 / 𝐿 . We set 𝛾 ( 𝐻 / 𝐿 ) 2 = 1 , and from this relation, we determine the lateral length scale 𝐿 to be √ 𝐿 = 𝛾 𝐻 2 =  𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 . Upon using the relation 𝛾 𝛿 2 = 1 , ( 4.9 ) becomes 𝑝 1 ( 0 ) − 𝑝 2 ( 0 ) + 1 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝜕 2 ℎ 𝜕 𝜒 2 + 𝜖 1 2  𝜕 𝜓 1  𝜕 𝑦 2 − 𝜖 2 2  𝜕 𝜓 2  𝜕 𝑦 2 . ( 4 . 1 0 ) The tangential stress continuity, ( 3.8 ), simplifies to the following condition in the long-wave limit: 𝜕 𝑢 2 ( 0 ) 𝜕 𝑦 − 𝜇 𝑟 𝜕 𝑢 1 ( 0 ) 𝜕 𝑦 = − 𝑞 𝜕 𝜓 1 𝜕 𝜒 . ( 4 . 1 1 ) The nondimensional kinematic condition at the interface, ( 3.10 ), after using the asymptotic expansion for 𝜐 𝑖 ( 4.7 ), yields 𝛿 2 𝜐 1 ( 0 ) ( 𝑦 = ℎ ( 𝜒 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝛿 2 𝑢 1 ( 0 ) 𝜕 ℎ 𝜕 𝜒 . ( 4 . 1 2 ) In order for the time derivative term in the above equation to be of the same order as the other two terms, it is necessary to stipulate a slow time scale in the long-wave limit such that 𝜕 𝜕 𝑡 = 𝛿 2 𝜕 𝜕 𝜏 , ( 4 . 1 3 ) where 𝜕 / 𝜕 𝜏 is 𝑂 ( 1 ) . The kinematic condition thus gives 𝜐 1 ( 0 ) ( 𝑦 = ℎ ( 𝜒 ) ) = 𝜀 𝜕 ℎ 𝜕 𝜏 + 𝑢 1 ( 0 ) 𝜕 ℎ 𝜕 𝜒 . ( 4 . 1 4 ) The dimensional slow time scale is obtained as follows: 𝜕 = 𝜇 𝜕 𝑡 2 𝐻 2 𝜖 0 𝜓 2 𝑏 𝜕 𝜕 𝑡 d i m = 𝛿 2 𝜕 , 𝜕 𝜏 ( 4 . 1 5 ) where 𝑡 d i m is the dimensional time. Thus, in order to nondimensionalize the dimensional time 𝑡 d i m in long-wave limit, the appropriate time scale is 𝜇 2 𝐻 2 / ( 𝜖 0 𝜓 2 𝑏 𝛿 2 ) . After using 𝛾 𝛿 2 = 1 to eliminate 𝛿 2 , the time scale becomes 𝜇 2 𝐻 3 𝛾 / ( 𝜖 0 𝜓 2 𝑏 ) 2 . Finally, the nondimensional interfacial charge balance, ( 3.11 ) is simplified in the long-wave limit as 𝜀 𝛿 2 𝜕 𝑞 𝜕 𝜏 + 𝛿 2 𝑢 1 ( 0 ) 𝜕 𝑞 𝜕 𝜒 − 𝛿 2 𝑞 𝜕 𝜐 1 ( 0 )  𝑆 𝜕 𝑦 = 𝜀 1 𝜕 𝜓 1 𝜕 𝑦 − 𝑆 2 𝜕 𝜓 2  . 𝜕 𝑦 ( 4 . 1 6 ) The above equation suggests that the nondimensional conductivities 𝑆 1 and 𝑆 2 both should scale as 𝛿 2 in order to balance the left side of the equation. So, we let 𝑆 𝑖 = 𝛿 2 𝑆 𝑖 ( 0 ) , 𝑖 = 1 , 2 , where 𝑆 𝑖 ( 0 ) = 𝜎 𝑖 𝜇 2 𝛾 𝐻 3 / ( 𝜖 3 0 𝜓 4 𝑏 ) . Upon using these rescaled conductivities, the charge conservation equation becomes, after using the continuity equation, ( 4.6 ), 𝜀 𝜕 𝑞 + 𝜕  𝑢 𝜕 𝜏 1 ( 0 ) 𝑞   𝑆 𝜕 𝜒 = 𝜀 1 ( 0 ) 𝜕 𝜓 1 𝜕 𝑦 − 𝑆 2 ( 0 ) 𝜕 𝜓 2  . 𝜕 𝑦 ( 4 . 1 7 ) This completes the derivation of the simplified governing equations in the long-wave limit. 5. Nonlinear Evolution Equations We now outline the derivation of the nonlinear evolution equations for the interfacial position ℎ ( 𝜒 , 𝜏 ) and surface charge density 𝑞 ( 𝜒 , 𝜏 ) . The simplified Laplacian for the potential 𝜓 𝑖 , ( 4.2 ), is easily solved along with the boundary conditions, ( 4.3 ), to give the following expressions for the potentials 𝜓 𝑖 ( 𝑖 = 1 , 2 ): 𝜓 1  𝜖 ( 𝜒 , 𝑦 ) = ( 𝛽 − 𝑦 ) 2 +  ( 1 + ℎ ( 𝜒 ) ) 𝑞 ( 𝜒 ) 𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  , 𝜓 ℎ ( 𝜒 ) 2 ( 𝜒 , 𝑦 ) = 𝛽 𝜖 2 − 𝜖 1  𝜖 𝑦 + 𝛽 ( 1 + 𝑦 ) 𝑞 ( 𝜒 ) + ℎ ( 𝜒 ) 1 − 𝜖 2 −  ( 1 + 𝑦 ) 𝑞 ( 𝜒 ) 𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  , ℎ ( 𝜒 ) ( 5 . 1 ) where the interfacial charge density 𝑞 ( 𝜒 , 𝜏 ) is determined below by the interface charge conservation equation ( 4.17 ). The simplified 𝑥 -momentum equation ( 4.8 ) can be integrated with respect to 𝑦 , since 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 is independent of 𝑦 and the two constants of integration that arise are determined by the boundary conditions 𝑢 1 ( 0 ) ( 𝛽 ) = 0 , 𝑢 1 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝑢 ( 0 ) i n t 𝑢 ( 𝜒 ) , 2 ( 0 ) ( − 1 ) = 0 , 𝑢 2 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝑢 ( 0 ) i n t ( 𝜒 ) . ( 5 . 2 ) Here, 𝑢 ( 0 ) i n t ( 𝜒 ) is the 𝑥 component of the velocity at the interface 𝑦 = ℎ ( 𝜒 ) , and this quantity will eventually be determined by using the tangential stress continuity condition ( 4.11 ). For the purposes of keeping the algebra tractable, it is found convenient to keep 𝑢 ( 0 ) i n t ( 𝜒 ) undetermined at present. The solutions for the 𝑥 -component velocities 𝑢 𝑖 ( 0 ) ( 𝑖 = 1 , 2 ) thus obtained are 𝑢 1 ( 0 ) 𝑘 ( 𝜒 , 𝑦 ) = 1 𝜇 𝑟 𝑑 𝑝 1 ( 0 ) ⎡ ⎢ ⎢ ⎣  √ 𝑑 𝜒 s i n h 𝜀 / 𝑘 1   √ ( 𝑦 − 𝛽 ) − s i n h 𝜀 / 𝑘 1  ( 𝑦 − ℎ ( 𝜒 ) )  √ s i n h 𝜀 / 𝑘 1 [ ]  ⎤ ⎥ ⎥ ⎦ ℎ ( 𝜒 ) − 𝛽 − 1 + 𝑢 ( 0 ) i n t  √ ( 𝜒 ) s i n h 𝜀 / 𝑘 1  ( 𝑦 − 𝛽 )  √ s i n h 𝜀 / 𝑘 1 [ ]  , 𝑢 ℎ ( 𝜒 ) − 𝛽 ( 5 . 3 ) 2 ( 0 ) ( 𝜒 , 𝑦 ) = 𝑘 1 𝑑 𝑝 2 ( 0 ) ⎡ ⎢ ⎢ ⎣  √ 𝑑 𝜒 s i n h 𝜀 / 𝑘 1   √ ( 𝑦 + 1 ) − s i n h 𝜀 / 𝑘 1  ( 𝑦 − ℎ ( 𝜒 ) )  √ s i n h 𝜀 / 𝑘 1  ⎤ ⎥ ⎥ ⎦ ( ℎ ( 𝜒 ) + 1 ) − 1 + 𝑢 ( 0 ) i n t  √ ( 𝜒 ) s i n h 𝜀 / 𝑘 1  ( 𝑦 + 1 )  √ s i n h 𝜀 / 𝑘 1  , ( ℎ ( 𝜒 ) + 1 ) ( 5 . 4 ) where the pressure gradients 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 ( 𝑖 = 1 , 2 ) are determined below. The continuity equation ( 4.6 ), after substituting the asymptotic expansion for 𝜐 𝑖 , ( 4.7 ), simplifies to 𝜕 𝑢 𝑖 ( 0 ) + 𝜕 𝜒 𝜕 𝜐 𝑖 ( 0 ) 𝜕 𝑦 = 0 . ( 5 . 5 ) The above equation is integrated with respect to 𝑦 from 𝑦 = 𝛽 to 𝑦 = ℎ ( 𝜒 ) for fluid 1 and 𝑦 = − 1 to 𝑦 = ℎ ( 𝜒 ) for fluid 2, to yield the following expressions for the normal velocities at the interface 𝜐 𝑖 ( 0 ) [ 𝑦 = ℎ ( 𝜒 ) ] for 𝑖 = 1 , 2 , where the integration constant is set to zero in order to satisfy the no-penetration condition at 𝑦 = 𝛽 and 𝑦 = − 1 : 𝜐 1 ( 0 ) [ ]  𝑦 = ℎ ( 𝜒 ) = − 𝛽 ℎ ( 𝜒 ) 𝜕 𝑢 1 ( 0 ) 𝜐 𝜕 𝜒 𝑑 𝑦 2 ( 0 ) [ ]  𝑦 = ℎ ( 𝜒 ) = − ℎ ( 𝜒 ) − 1 𝜕 𝑢 2 ( 0 ) 𝜕 𝜒 𝑑 𝑦 . ( 5 . 6 ) Note that the normal velocities of the two fluids are equal at the interface (normal velocity continuity condition), and so, 𝜐 1 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝜐 2 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝜐 ( 0 ) i n t . ( 5 . 7 ) Therefore, we equate the two integrals in the above equation and apply Leibnitz rule 𝜕  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 − 𝜕 ℎ 𝑢 𝜕 𝜒 1 ( 0 ) ℎ 𝜕 ( 𝜒 ) =  𝜕 𝜒 ℎ ( 𝜒 ) − 1 𝑢 2 ( 0 ) 𝑑 𝑦 − 𝜕 ℎ 𝑢 𝜕 𝜒 2 ( 0 ) ℎ ( 𝜒 ) . ( 5 . 8 ) Noting that the 𝑥 -component velocities are equal at the interface, 𝑢 1 ( 0 ) [ ℎ ( 𝜒 ) ] = 𝑢 2 ( 0 ) [ ℎ ( 𝜒 ) ] , the above equation is simplified to 𝜕   𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 )  𝑑 𝑦 − ℎ ( 𝜒 ) − 1 𝑢 2 ( 0 )  𝑑 𝑦 = 0 . ( 5 . 9 ) We next integrate the above equation with respect to 𝜒 , and set the integration constant (which is at most a function of time) to zero. The constant of integration is zero, because the pressure gradients 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 in the two fluids should be zero in the absence of electric field, and when the interfaces are flat. We then substitute the expressions for 𝑢 𝑖 ( 0 ) ( 𝑦 ) , ( 5.3 ) and ( 5.4 ), in the above equation, carry out the integrations with respect to 𝑦 , and substitute the simplified normal stress continuity condition ( 4.10 ) to eliminate 𝑑 𝑝 2 ( 0 ) / 𝑑 𝜒 in terms of 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 . Prior to determining 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 , it is useful to determine the 𝑥 -component of the fluid velocity at the interface 𝑢 ( 0 ) i n t from the simplified tangential stress continuity condition, ( 4.11 ). Once 𝑢 ( 0 ) i n t is determined, the pressure gradient 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 is determined from the integrated version of ( 5.9 ), and thus the velocity profile (( 5.3 ) and ( 5.4 )), is known completely. Therefore, we get 𝑑 𝑝 1 ( 0 ) =  𝑘 𝑑 𝜒 1  ℎ ( 1 + ℎ ( 𝑥 ) )      ( 𝜒 ) − 𝑀 ( 𝜒 ) − 𝐺 𝐹 𝑘 1 𝜀 − 𝑘 1  ℎ    (    𝜒 ) − 𝑀 ( 𝜒 ) 𝐺 + 4 𝑘 1 𝜀    t a n h 𝜀 4 𝑘 1 [ ]    𝑘 1 + ℎ ( 𝜒 ) 1  1 1 + ℎ ( 𝜒 ) + 𝜇 𝑟 [ ]    𝛽 − ℎ ( 𝜒 ) − 𝐺 t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) − t a n h 𝜀 4 𝑘 1 [ ] −  ℎ ( 𝜒 ) − 𝛽   4 𝑘 1 𝜀    t a n h 𝜀 4 𝑘 1 [ ]  − 1 1 + ℎ ( 𝜒 ) 𝜇 𝑟   t a n h 𝜀 4 𝑘 1 [ ] ℎ ( 𝜒 ) − 𝛽     − 1 , ( 5 . 1 0 ) 𝑑 𝑝 2 ( 0 ) = 𝑑 𝜒 𝑑 𝑝 1 ( 0 ) ( 𝜒 ) − 𝜕 𝑑 𝜒 3 ℎ ( 𝜒 ) 𝜕 𝜒 3 + 𝑢 ( 0 ) i n t ( 𝜒 ) 𝑘 1  𝜇 𝑟  𝑢 − 1 + 𝑀 ( 𝜒 ) , ( 5 . 1 1 ) ( 0 ) i n t   ( 𝜒 ) = 𝑘 1 𝜀 𝐹 − 𝑘 1 𝑑 𝑝 1 ( 0 )    𝑑 𝜒 t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) − t a n h 𝜀 4 𝑘 1 [ ] ℎ ( 𝜒 ) − 𝛽   + 𝑘 1   t a n h 𝜀 4 𝑘 1 [ ] 𝜕 1 + ℎ ( 𝜒 )   3 ℎ ( 𝜒 ) 𝜕 𝜒 3 ×    − 𝑀   c o t h 𝜀 𝑘 1 [ ]  1 + ℎ ( 𝜒 ) − 𝜇 𝑟   c o t h 𝜀 𝑘 1 [ ]  +  𝜇 ℎ ( 𝜒 ) − 𝛽 𝑟    − 1 t a n h 𝜀 4 𝑘 1 [ ] 1 + ℎ ( 𝜒 )   − 1 , ( 5 . 1 2 ) where 𝐺 , 𝑀 , and 𝐹 are functions of 𝜒 , and they are defined as    𝐺 ( 𝜒 ) = c o t h 𝜀 𝑘 1 [ ]  1 + ℎ ( 𝜒 ) − 𝜇 𝑟   c o t h 𝜀 𝑘 1 [ ]  +  𝜇 ℎ ( 𝜒 ) − 𝛽 𝑟    − 1 t a n h 𝜀 4 𝑘 1 [ ] 1 + ℎ ( 𝜒 )   − 1 ×   𝜇 𝑟  [ ] +  − 1 1 + ℎ ( 𝜒 ) 𝑘 1 𝜀   1 − 2 𝜇 𝑟    t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) + t a n h 𝜀 4 𝑘 1 [ ]  𝜖 ℎ ( 𝜒 ) − 𝛽    𝑀 ( 𝜒 ) = 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) − 3  𝜖 1 𝜖 2 ℎ   ( 𝜒 ) 𝑞 ( 𝜒 ) ( 1 + 𝛽 ) + 𝜖 2 − 𝜖 1 ×  𝜖   2 + 𝜖 1 [ ]  +  𝜖 + 𝑞 ( 𝜒 ) 1 + 2 ℎ ( 𝜒 ) − 𝛽 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) − 2 ×  𝑞 ( 𝜒 ) 𝑞   𝜖 ( 𝜒 ) 1 [ ] 1 + ℎ ( 𝜒 ) 2 + 𝜖 2 [ ] 𝛽 − ℎ ( 𝜒 ) 2  + 𝜖 1 𝜖 2 𝑞  [ ]  , 𝐹 ( 𝜒 ) 1 + 2 ℎ ( 𝜒 ) − 𝛽 ( 𝜒 ) = − 𝑞 ( 𝜒 ) ℎ  ( 𝜒 ) 𝜖 2  ( 𝛽 − ℎ ( 𝜒 ) ) 𝑞 ( 𝜒 ) ( 𝛽 + 1 ) − 𝜖 1 + 𝜖 2   𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) 2 − 𝑞 ( 𝜒 ) 𝑞  ( 𝜒 ) ( 𝛽 − ℎ ( 𝜒 ) ) ( 1 + ℎ ( 𝜒 ) )  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   . ℎ ( 𝜒 ) ( 5 . 1 3 ) Finally, the kinematic condition at the interface, ( 4.14 ), is used to derive the evolution equation for ℎ ( 𝜒 , 𝜏 ) , as follows. We first substitute the expressions for the normal fluid velocities at the interface, ( 5.6 ), after using Leibnitz rule on the integral 𝜀 𝜕 ℎ 𝜕 𝜏 + 𝑢 1 ( 0 ) [ ] ℎ ( 𝜒 ) 𝜕 ℎ 𝜕 𝜕 𝜒 = −  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 + 𝑢 1 ( 0 ) [ ] ℎ ( 𝜒 ) 𝜕 ℎ , 𝜕 𝜒 ( 5 . 1 4 ) which is simplified to give 𝜀 𝜕 ℎ + 𝜕 𝜕 𝜏  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 = 0 , ( 5 . 1 5 ) where 𝑢 1 ( 0 ) is substituted from ( 5.3 ), after using the expressions for 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 and 𝑢 ( 0 ) i n t determined using the procedure outlined above. The evolution equation for the interfacial charge density 𝑞 ( 𝜒 , 𝜏 ) is obtained from ( 4.17 ) after substituting the expressions for 𝑢 ( 0 ) i n t and the gradients of the potential (from ( 5.1 )) in that equation. The nonlinear evolution equations for ℎ ( 𝜒 , 𝜏 ) and 𝑞 ( 𝜒 , 𝜏 ) , respectively, take the forms 𝜀 𝜕 ℎ + 1 𝜕 𝜏 𝜇 𝑟 𝑘   1 𝑑 𝑝 1 ( 0 ) + 𝜇 𝑑 𝜒 𝑟 2 𝑢 ( 0 ) i n t  ( 𝜒 ) 𝜕 ℎ 𝜕 𝜒 s e c h 2  1 2  𝜀 𝑘 1  ( 𝛽 − ℎ ) + 𝑘 1  𝑑 ( 𝛽 − ℎ ) 2 𝑝 1 ( 0 ) 𝑑 𝜒 2 − 𝑑 𝑝 1 ( 0 ) 𝑑 𝜒 𝜕 ℎ  −  𝜕 𝜒 𝑘 1 𝜀  2 𝑘 1 𝑑 2 𝑝 1 ( 0 ) 𝑑 𝜒 2 + 𝜇 𝑟 𝑢  ( 0 ) i n t   1 ( 𝜒 ) × t a n h 2  𝜀 𝑘 1   𝜀 ( 𝛽 − ℎ ) = 0 , 𝜕 𝑞 𝜕 𝜏 + 𝑢 ( 0 ) i n t ( 𝜒 ) 𝜕 𝑞 𝜕 𝜒 − 𝑞 𝜕 ℎ  √ 𝜕 𝜒 𝜀 𝑘 1 𝜇 𝑟 𝑑 𝑝 1 ( 0 )  1 𝑑 𝜒 t a n h 2  𝜀 𝑘 1  +  ( ℎ − 𝛽 ) 𝜀 𝑘 1 𝑢 ( 0 ) i n t   ( 𝜒 ) c o t h 𝜀 𝑘 1   ( ℎ − 𝛽 ) + 𝑞 𝑢  ( 0 ) i n t + ( 𝜒 ) 𝜀 𝑆 1 ( 0 )  ( 1 + ℎ ) 𝑞 + 𝜖 2  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  ℎ + 𝜀 𝑆 2 ( 0 )  ( 𝛽 − ℎ ) 𝑞 − 𝜖 1  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  ℎ = 0 . ( 5 . 1 6 ) These coupled nonlinear equations can be solved numerically with appropriate initial conditions to determine the evolution of the interface in the presence of electric fields and porous medium. In the present work, however, we restrict ourselves to studying the linear stability properties of these equations, an issue we turn to next. 6. Stability Analysis and Discussion Before linearizing the coupled nonlinear equations, it is necessary to first determine the base state about which we perturb. The steady base state we consider is that of stationary fluids with a flat interface ℎ ( 𝜒 , 𝜏 ) = 0 , and with a constant interfacial charge density 𝑞 0 which is independent of 𝜒 and 𝜏 . This base state interfacial charge 𝑞 0 is determined from ( 4.17 ) with the left side set to zero, since 𝜕 / 𝜕 𝜏 = 0 and 𝑢 𝑖 ( 0 ) = 0 in the base state 𝑆 1 ( 0 )  𝜕 𝜓 1  𝜕 𝑦 𝑦 = 0 = 𝑆 2 ( 0 )  𝜕 𝜓 2  𝜕 𝑦 𝑦 = 0 . ( 6 . 1 ) The derivatives of the potentials in the above equation are calculated from ( 5.1 ), with ℎ ( 𝜒 ) set to zero. This yields the following expression for the steady interfacial charge density: 𝑞 0 = 𝜖 1 𝑆 2 ( 0 ) − 𝜖 2 𝑆 1 ( 0 ) 𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 ) . ( 6 . 2 ) The variables ℎ and 𝑞 are now perturbed about their base state values ℎ ( 𝜒 , 𝜏 ) = ℎ 1 [ ] , 𝑞 e x p 𝑖 𝑘 𝜒 + 𝜔 𝜏 ( 𝜒 , 𝜏 ) = 𝑞 0 + 𝑞 1 [ ] , e x p 𝑖 𝑘 𝜒 + 𝜔 𝜏 ( 6 . 3 ) where ℎ 1 and 𝑞 1 are the amplitudes of the perturbations which are independent of 𝜒 and 𝜏 , 𝑘 is the nondimensional wavenumber based on the lateral length scale  𝐿 = 𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 and 𝜔 is the nondimensional growth rate based on the time scale 𝜇 2 𝐻 3 𝛾 / ( 𝜖 0 𝜓 2 𝑏 ) 2 . We substitute the expressions for 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 and 𝑢 ( 0 ) i n t from ( 5.10 ), ( 5.12 ), and ( 6.3 ) in the nonlinear evolution equations ( 5.16 ), and apply Taylor expansion on the hyperbolic functions about ℎ ( 𝜒 ) = 0 . We then linearize the resulting coupled nonlinear evolution equations with respect to ℎ 1 and 𝑞 1 , to obtain a set linear homogeneous equations for ℎ 1 and 𝑞 1 of the form  𝐻 1 ℎ 1 +  𝑄 1 𝑞 1  𝐻 = 0 , 2 ℎ 1 +  𝑄 2 𝑞 1 = 0 , ( 6 . 4 ) where  𝐻 1  𝑅 = − 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  𝛽 4 𝜖 4 2 + 𝜖 4 1  𝑅 1 4 + 𝑘 4 𝑅 1 3 𝑅 1 − 𝑘 2 𝑅 1 3 𝑅 1 𝜖 2  − 𝜖 2 1 𝜖 2 2  6  𝑅 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  𝛽 2 + 𝑘 2 𝑅 1 3 𝑅 1 𝜖 2  + 𝜖 3 1 𝜖 2  𝑅 1 5  𝑅 − 4 𝛽 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  + 𝑘 2 𝑅 1 3 𝑅 1 𝛽 𝜖 2  − 𝜖 1 𝜖 3 2  𝑅 1 5 + 4 𝛽 3  𝑅 1 4 + 𝑘 4 𝑅 1 3 R 1  + 𝑘 2 𝑅 1 3 𝑅 1 𝛽 𝜖 2  + 𝑞 2 0 ( 1 + 𝛽 ) 𝜖 2  − 𝑅 6 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  × 𝑅   1 𝐶 ( 𝛽 − 1 ) − 1 𝐵 1 𝛽  𝑘 1 𝜀  𝜖 1 − 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2   − 𝑞 0 𝜖 2  𝑅 8 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  ×  𝛽  2 𝑅 1 − 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 1 +  2 𝑅 1 − 𝐶 1 𝐵 1  𝛽 ( 𝛽 − 1 ) 𝑘 1 𝜀  𝜖 1 𝜖 2 + 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 2   + 𝐵 1 𝑅 5 𝜔 𝜀 3 / 2 𝜇 𝑟 ,  𝑄 1 = ( 𝛽 − 1 ) 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2 𝑅   1 5 + 𝑘 2 𝑅 1 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2   − 𝑞 0  𝑅 7 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  2 × 𝑅   1 𝐶 + 𝛽 1 𝐵 1  𝑘 1 𝜀  𝜖 1 + 𝛽 2  𝑅 1 + 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 ,  𝐻   2 =  𝜖 1 + 𝛽 𝜖 2   𝜀  𝜀 𝑅 2  𝑅 1 8 − 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 4 1 𝜇 r + 𝜖 3 1  𝑅 1 6 + 𝑅 3 𝑅 4 + 𝐵 1 𝑅 2 𝑅 4 𝑅 3 4  + 𝜖 2 1 𝜖 2  𝑅 1 7 + 3 𝛽 𝑅 3 𝑅 4 + 3 𝐵 1 𝑅 2 𝑅 4 𝑅 3 5  + 𝛽 𝜖 2  𝑅 3 0 𝑅 2 2  k 1 𝜀 q 0 𝜇 r + 𝛽 𝜖 2 2  𝑅 1 9 + 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 R 2 ×  𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑅 1 8 + 𝐵 1 𝑅 4 𝑆 1 ( 0 )  𝜖 2  𝜇 𝑟   + 𝜖 1  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2  𝑅 2 0 + 3 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  3 𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  4 𝛽 𝑅 1 8 + 3 𝐵 1 𝑅 4 𝑆 1 ( 0 ) − 𝛽 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 2  𝜇 𝑟    − 𝑅 2 𝑞 0 𝜇 𝑟  𝑅 9  𝑅 2 1 + 𝑅 2 3 + 𝑅 2 5 + 𝑅 2 8 − 𝑅 2 9  + 𝑅 3 2 𝑅 3 6   ,  𝑄 2 = 𝑅 2  𝜖 1 + 𝛽 𝜖 2  𝜇 𝑟  𝐵 1 𝑅 2 𝜀 2  𝜖 1 + 𝛽 𝜖 2  3  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  𝜖 + 𝜔 1 + 𝛽 𝜖 2   − 𝑞 0  𝑅 9  𝑅 2 2 + 𝑅 2 4 + 𝑅 2 6 + 𝑅 2 7  − 𝑅 3 1 𝑅 2 √ 𝜀 𝑘 1  𝜖 1 + 𝛽 𝜖 2  + 𝑅 3 3 𝑅 3 6 ,   ( 6 . 5 ) in which 𝐵 1 , 𝐵 2 , 𝐶 1 , 𝐶 2 , and 𝑅 1 − 𝑅 4 9 are given in the appendix. The set of linear homogeneous equations ( 6.4 ) for ℎ 1 and 𝑞 1 is written in the matrix form 𝑀 ⋅ 𝐶 𝑇 = 0 , where the vector 𝐶 = [ ℎ 1 , 𝑞 1 ] and the determinant of this matrix 𝑀 is set to zero for nontrivial solutions in order to obtain the characteristic equation for 𝜔 . This characteristic equation, which gives the growth rate 𝜔 as a function of 𝑘 , 𝑘 1 , 𝛽 , 𝜇 𝑟 , 𝑆 1 ( 0 ) , 𝑆 2 ( 0 ) , 𝜀 , 𝜖 1 , and 𝜖 2 , is a quadratic equation for 𝜔 . The roots of the characteristic equation for 𝜔 can be written as 𝜔 1 𝑅 = − 4 8 𝑅 4 9 −  𝑅 2 4 8  𝐻 + 2 2  𝑄 1 𝑅 4 9 − 2 𝑅 4 1 𝑅 4 7 𝑅 4 9 𝑅 4 9 , 𝜔 ( 6 . 6 ) 2 𝑅 = − 4 8 𝑅 4 9 +  𝑅 2 4 8  𝐻 + 2 2  𝑄 1 𝑅 4 9 − 2 𝑅 4 1 𝑅 4 7 𝑅 4 9 𝑅 4 9 . ( 6 . 7 ) The roots 𝜔 1 and 𝜔 2 are always real, with one of them 𝜔 1 is always negative, and the other one 𝜔 2 can be positive or negative depending on the choice of the system parameters. We can take 𝜔 = 𝜔 2 as a function of 𝑘 only in ( 6.7 ) if the values of the other parameters are known. Now, to see the effects of various parameters on the stability of the considered system, we calculate the growth rate 𝜔 given by ( 6.7 ) as a function of the wavenumber 𝑘 for different values of all physical parameters included in the analysis. These calculations are presented in Figures 2 – 7 for the general case of leaky-leaky dielectric fluids, Figures 8 and 9 for the limiting case of perfect-leaky dielectric fluids in which 𝑆 1 ( 0 ) = 0 , and Figures 10 – 13 for another limiting case of perfect-perfect dielectric fluids in which 𝑆 1 ( 0 ) = 𝑆 2 ( 0 ) = 0 , where we have given the growth rate 𝜔 against the wavenumber 𝑘 for the porosity of porous medium 𝜀 , medium permeability 𝑘 1 , nondimensional conductivities 𝑆 1 ( 0 ) , 𝑆 2 ( 0 ) , dielectric constants 𝜖 1 , 𝜖 2 , the ratio of the thicknesses of top and bottom fluids 𝛽 , and the ratio of viscosities of the two fluids 𝜇 𝑟 = 𝜇 1 / 𝜇 2 , respectively. Figure 2: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 2 0 and (b) 0 ≤ 𝑘 ≤ 5 . Figure 3: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑘 1 when 𝜀 = 0 . 2 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑘 1 = 5 , 1 0 , 1 2 , respectively. Figure 4: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑆 1 ( 0 ) when 𝜀 = 0 . 4 , 𝑘 1 = 1 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑆 1 ( 0 ) = 1 0 2 , 𝑆 1 ( 0 ) = 1 0 3 , 1 0 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 1 0 0 and (b) 0 ≤ 𝑘 ≤ 1 0 . Figure 5: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑆 2 ( 0 ) when 𝜀 = 0 . 4 , 𝑘 1 = 1 0 , 𝑆 1 ( 0 ) = 1 0 3 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑆 2 ( 0 ) = 1 0 2 , 1 0 4 , 1 0 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 1 0 0 and (b) 0 ≤ 𝑘 ≤ 3 . Figure 6: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 1 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 1 = 3 , 5 , 7 , respectively. Figure 7: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜇 𝑟 when 𝜀 = 0 . 4 , 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , and 𝜇 𝑟 = 1 , 2 , 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 3 0 , and (b) 0 ≤ 𝑘 ≤ 5 . Figure 8: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 2 0 , and (b) 0 ≤ 𝑘 ≤ 5 . Figure 9: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜇 𝑟 when 𝜀 = 0 . 4 , 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , and 𝜇 𝑟 = 1 , 2 , 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 3 0 and (b) 0 ≤ 𝑘 ≤ 1 0 . Figure 10: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively. Figure 11: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑘 1 when 𝜀 = 0 . 2 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑘 1 = 5 , 1 0 , 1 2 , respectively. Figure 12: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 1 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 1 = 3 , 5 , 7 , respectively. Figure 13: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 2 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 2 = 4 , 6 , 1 0 , respectively. 6.1. Leaky Dielectric-Leaky Dielectric Interface This is the general case in which 𝑆 1 ( 0 ) ≠ 0 and 𝑆 2 ( 0 ) ≠ 0 , which is discussed in Figures 2 – 7 , when the conductivities of the upper and lower fluids are present in the analysis. Figures 2(a) and 2(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the porosity of porous medium 𝜀 . It is clear from Figure 2(b) that for small values of the porosity (e.g., 𝜀 = 0 . 1 ), the growth rate 𝜔 increases by increasing the wavenumber 𝑘 till a maximum value of 𝜀 after which 𝜔 decreases by increasing 𝑘 ; that is, there exists a maximum mode of instability only for small values of the porosity 𝜀 , while for any other value of 𝜀 , we found that 𝜔 decreases by increasing 𝑘 , that is. the system is always stable in this case. It is clear also from Figure 2(a) , by increasing the porosity values and at any wavenumber value, that the porosity of porous medium has a stabilizing effect for the wavenumber range 0 ≤ 𝑘 ≤ 1 5 , and it has a destabilizing effect for higher wavenumber values 𝑘 > 1 5 ; that is, the porosity of porous medium has a dual role on the stability of the considered system (stabilizing and then destabilizing) depends on the wavenumber range (less than or hifher than a critical wavenumber value 𝑘 = 1 5 , resp.). Figure 3 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the medium permeability 𝑘 1 . We conclude from this figure that for any wavenumber value 𝑘 ≤ 2 0 , the growth rate 𝜔 increases by increasing the medium permeability 𝑘 1 values, while for wavenumber values 𝑘 > 2 0 , all the curves correspond to different values of medium permeability 𝑘 1 coincide. This means that the medium permeability 𝑘 1 has a destabilizing effect for the wavenumber range 0 < 𝑘 ≤ 2 0 and it has no effect on the stability of the considered system for higher wavenumber values 𝑘 > 2 0 . Figures 4(a) and 4(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the conductivity of upper fluid 𝑆 1 ( 0 ) . It is clear from this figure that the conductivity 𝑆 1 ( 0 ) has a stabilizing effect for small wavenumber values and a destabilizing effect for high wavenumber values, as the growth rates 𝜔 decrease and increase by increasing the conductivity 𝑆 1 ( 0 ) values, respectively. It is also seen that between these two wavenumber ranges, the conductivity of upper fluid 𝑆 1 ( 0 ) has a dual role on the stability of the considered system; that is, it has a destabilizing effect for 𝑆 1 ( 0 ) values greater than 1 0 2 , while it has a stabilizing effect for 𝑆 1 ( 0 ) values greater than 1 0 4 . Therefore, we conclude that the conductivity of upper fluid 𝑆 1 ( 0 ) has different effects on the stability of the system depending on the chosen wavenumber range. Figures 5(a) and 5(b) shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the conductivity of lower fluid 𝑆 2 ( 0 ) . It indicates that the conductivity 𝑆 1 ( 0 ) has a destabilizing effect for the wavenumber range 0 < 𝑘 < 1 0 , while it has a stabilizing effect for higher wavenumber values 𝑘 ≥ 1 0 , since the growth rate 𝜔 increases in the first wavenumber range and decreases in the second wavenumber range by increasing the conductivity of lower fluid 𝑆 2 ( 0 ) . We conclude also from Figure 5(b) that there exists a mode of maximum instability for small values of the conductivity 𝑆 2 ( 0 ) which disappears for high values of the conductivity of lower fluid 𝑆 2 ( 0 ) ≥ 1 0 4 . Figure 6 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the dielectric constant of upper fluid 𝜖 1 . It is seen from this figure that there exists a critical wavenumber value 𝑘 = 4 . 6 , before which the growth rates decrease by increasing the dielectric constant 𝜖 1 , and after which they increase by increasing 𝜖 1 values. Thus, the dielectric constant of upper fluid 𝜖 1 has a stabilizing as well as a destabilizing effects, for wavenumber ranges before and after this critical wavenumber value, respectively. Similarly, the effects of both the dielectric constant of lower fluid 𝜖 2 , and the ratio of thicknesses of upper and lower fluids 𝛽 on the stability of the considered system are found to have opposite effects to the effect of the dielectric constant 𝜖 1 given by Figure 6 , but the corresponding figures are not given here. In other words, we conclude that the dielectric constant of lower fluid 𝜖 2 has a destabilizing as well as a stabilizing effects for wavenumber ranges before and after the same critical wavenumber value 𝑘 = 4 . 6 , respectively, while the ratio of thicknesses of upper and lower fluids 𝛽 has also a destabilizing as well as a stabilizing effects for wavenumber ranges before and after a less critical wavenumber value 𝑘 = 2 . 5 , respectively. Figures 7(a) and 7(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the ratio of viscosity of upper and lower fluids 𝜇 𝑟 . It indicates that the viscosity ratio 𝜇 𝑟 has a stabilizing effect for small wavenumber values, and it has also a destabilizing effect for higher wavenumber values, since the growth rate 𝜔 decreases and increases by increasing the increase of 𝜇 𝑟 , respectively. It is clear also from Figure 7(b) that for 𝜇 𝑟 > 1 , that is, when the viscosity of upper fluid is larger than the viscosity of lower fluid, there exists a mode of maximum instability which disappear for 𝜇 𝑟 ≤ 1 . 6.2. Perfect Dielectric-Leaky Dielectric Interface This is the limiting case in which 𝑆 1 ( 0 ) = 0 and 𝑆 2 ( 0 ) ≠ 0 , which is discussed in Figures 8 and 9 , when the conductivity of the upper fluid is not included in the analysis. Figures 8(a) and 8(b) show the variation of the growth rate 𝜔 against the wavenumber 𝑘 for various values of the porosity of porous medium 𝜀 . In comparison with Figures 2(a) and 2(b) , we conclude that the porosity of porous medium 𝜀 behaves in the same manner as in the previous case of leaky-leaky dielectric fluids, except that the values of growth rates 𝜔 are higher than their values in the previous case and the corresponding curves intersect the 𝑘 -axis at larger wavenumber values than in the previous case. We conclude also that for small wavenumber values, there exists a mode of maximum instability for porosity values 𝜀 ≤ 0 . 3 . Similarly, the effects of medium permeability 𝑘 1 , the conductivity of lower fluid 𝑆 2 ( 0 ) , the dielectric constants 𝜖 1 , 𝜖 2 , and the ratio of thicknesses of upper and lower fluids 𝛽 on the stability of the considered system are found to behave in the same manner as their effects in the previous case of leaky-leaky dielectric fluids, but figures are excluded to save space, except that in this case: for the effect of 𝑘 1 , the growth rate values are higher than their values in the previous case; for the effect of 𝑆 2 ( 0 ) , the obtained curves intersect the 𝑘 -axis at larger wavenumber values than in the previous case; for the effect of 𝜖 1 , there exists a mode of maximum instability for all values of 𝜖 1 ; for the effect of 𝜖 2 , the obtained curves are exactly similar to those obtained in the previous case; finally, for the effect of 𝛽 , there is a mode of maximum instability and the obtained curves intersect the 𝑘 -axis at bigger wavenumber values than those obtained in the previous case. Figures 9(a) and 9(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the viscosity ratio of upper and lower fluids 𝜇 r . In comparison with Figures 7(a) and 7(b) , we conclude that the viscosity ratio 𝜇 𝑟 behaves in the same manner as in the previous case of leaky-leaky dielectric fluids with the only difference that in this case, for all values of 𝜇 𝑟 ⪋ 1 , there exists a mode of maximum instability and not only for 𝜇 𝑟 > 1 shown in the previous case. 6.3. Perfect Dielectric-Perfect Dielectric Interface This is the limiting case in which 𝑆 1 ( 0 ) = 𝑆 2 ( 0 ) = 0 , which is discussed in Figures 10 – 13 , when the conductivities of the upper and lower fluids are absent in the analysis. Figure 10 shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for different values of the porosity of porous medium 𝜀 . In view of the above discussion, we conclude from this figure that the porosity 𝜀 has a slightly stabilizing effect for small wavenumber values 𝑘 ≤ 2 . 5 and that it has no effect on the stability of the considered system afterwards for higher wavenumber values. Figure 11 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for various values of the medium permeability 𝑘 1 . It is clear from this figure that the permeability 𝑘 1 has a destabilizing effect, since the growth rate 𝜔 increases by increasing the medium permeability values at any fixed wavenumber value. Similarly, the effect of dielectric constant of the upper fluid 𝜖 1 on the stability of the system is illustrated in Figure 12 , and it shows that the dielectric constant 𝜖 1 has a stabilizing effect since the growth rate 𝜔 decreases by increasing the dielectric constant 𝜖 1 at any wavenumber value. Figure 13 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the dielectric constant of lower fluid 𝜖 2 , and from which we conclude thatthe dielectric constant 𝜖 1 , and after which they increase by increasing 𝜖 2 has a destabilizing effect in the wavenumber range 0 < 𝑘 ≤ 4 . 6 , and it has no effect on the stability of the considered system afterwards for higher wavenumber values 𝑘 > 4 . 6 . Similarly, the effects of both the ratio of thicknesses of upper and lower fluids 𝛽 , and the ratio of viscosities of upper and lower fluids 𝜇 𝑟 , respectively, on the stability of the system are found to has the same and the opposite effects as the effect of the dielectric constant 𝜖 2 shown in Figure 13 , but the corresponding figures are not given here to avoid any kind of repitation; that is, the parameters 𝛽 and 𝜇 𝑟 have destabilizing and stabilizing effects for small wavenumber values, respectively. 7. Concluding Remarks In conclusion, we have provided a general formulation for analyzing the effect of an externally applied electric field on the stability and dynamics of the interface between two leaky dielectric fluids of arbitrary viscosities and conductivities in porous medium. A systematic long-wave asymptotic analysis was used to derive coupled nonlinear evolution equations for the position of the interface and free charge density at the interface. Attention was restricted to linearized stability of the coupled nonlinear equations and the effect of a variety of system parameters on the stability of the considered system. Two limiting cases are also studied, that is, the case of perfect-leaky dielectric fluids and the case of two perfect dielectric fluids, and recovered the previous studies in absence of porous medium. The obtained results in these limiting cases and the general case of two leaky dielectric fluids can be summarized as follows: (I) For perfect-perfect dielectric fluids, we conclude for small wavenumbers that (i) the porosity of porous medium 𝜀 , the dielectric constant of upper fluid 𝜖 1 , and the ratio of viscosity of upper and lower fluids 𝜇 𝑟 have stabilizing effects, (ii) the medium permeability 𝑘 1 , the dielectric constant of lower fluid 𝜖 2 , and the ratio of thickness of upper and lower fluids 𝛽 have destabilizing effects, (iii) at any value of these physical parameters there are no modes of maximum instability, that is. the system is always stable, (iv) these physical parameters have no effect on the stability of the system for high wavenumber values. (II) For perfect-leaky dielectric fluids, we found that (i) the conductivity of lower fluid S 2 ( 0 ) has a destabilizing effect for small wavenumbers and a stabilizing effect for high wavenumbers, (ii) for small wavenumber values, the physical parameters 𝜀 , 𝜖 1 , and 𝜇 𝑟 have stabilizing effects, while the parameters, 𝑘 1 , 𝜖 2 , and 𝛽 have destabilizing effects, as in the case of perfect-perfect dielectrics, (iii) for high wavenumber values, new regions of stability or instability appear; in other words, the physical parameters 𝜀 , 𝜖 1 , and 𝜇 r have destabilizing effects for high wavenumbers, while the parameters 𝜖 2 , and 𝛽 have stabilizing effects, and k 1 has no effect on the stability of the system for high wavenumber values, (iv) there exists a mode of maximum instability for some of these physical parameters which do not appear in the previous case of perfect-perfect dielectrics. (III) For leaky-leaky dielectric fluids, we found that (i) the conductivity of upper fluid S 1 ( 0 ) has a destabilizing effect for small wavenumbers 0 < 𝑘 < 𝑘 1 and also a stabilizing effect for high wavenumbers 𝑘 > 𝑘 2 , while it has a dual role on the stability of the considered system between them in the wavenumber range 𝑘 1 < 𝑘 < 𝑘 2 . (ii) the effects of all other physical parameters on the stability of the considered system behave in the same manner as their effects in the case of perfect-leaky dielectrics, except that in the case of perfect-leaky dielectrics the stability or instability regions occur more faster than the corresponding case of leaky-leaky dielectrics, and the maximum instability holds for more values of the physical parameters included in the analysis. It should be mentioned that the problem investigated in this article can be generalized to study the linear electrohydrodynamic instabilities at the interface between two immiscible fluids, either perfect or leaky dielectrics, subjected to alternating electric fields and moving through a porous medium in the limit of the electrode spacing being large compared to the wavelength of the perturbation using the Floquet theory analysis [ 37 ], and this case is now in a current research. Appendix 𝐵 1   = c o t h 𝜀 𝑘 1  + 𝜇 𝑟  𝛽  c o t h 𝜀 𝑘 1  +  𝜇 𝑟   1 − 1 t a n h 2  𝜀 𝑘 1  , 𝐵 2 =  𝜀 𝑘 1  − c s c h 2   𝜀 𝑘 1  + 𝜇 𝑟 c s c h 2  𝛽  𝜀 𝑘 1  + 1 2  𝜇 𝑟  − 1 s e c h 2  1 2  𝜀 𝑘 1 , 𝐶   1 = 𝜇 𝑟  − 1 − 𝑘 1 𝜀   2 𝜇 𝑟   1 − 1 t a n h 2  𝜀 𝑘 1   1 + t a n h 2 𝛽  𝜀 𝑘 1 , 𝐶   2 = 𝜇 𝑟 1 − 1 + 2   1 + 1 − 2 𝜇 𝑟  s e c h 2  1 2  𝜀 𝑘 1  − t a n h 2  1 2 𝛽  𝜀 𝑘 1 , 𝑅   1 = 𝑘 1 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀   1 t a n h 2  𝜀 𝑘 1   , 𝑅 − 1 2 = 1 𝜇 𝑟   𝛽 − 2 𝑘 1 𝜀  𝛽 t a n h 2  𝜀 𝑘 1  − 𝜇 𝑟  𝑅 1 𝑘 1 + 𝐶 1 𝐵 1  𝛽 t a n h 2  𝜀 𝑘 1  , 𝑅   3 = 𝐵 1 𝑅 2 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  − 𝜖 1 𝑆 2 ( 0 ) + 𝜖 2 𝑆 1 ( 0 )  , 𝑅 4  = − − 𝐵 2 𝐶 1 + 𝐵 1 𝐶 2 𝐵 2 1    1 t a n h 2  𝜀 𝑘 1   1 + t a n h 2 𝛽  𝜀 𝑘 1 , 𝑅   5 = 𝑅 2 2  𝜖 1 + 𝛽 𝜖 2  4 , 𝑅 6 = 𝜀 𝜖 1  1 ( 1 − 𝛽 ) s i n h 2  𝜀 𝑘 1   𝜖 + 𝛽 1 + 𝛽 𝜖 2   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1  , 𝑅 7 =  𝜖 1 + 𝛽 𝜖 2   𝜀  𝜖 1 + 𝛽 2 𝜖 2   1 s i n h 2  𝜀 𝑘 1   𝜖 + 𝛽 1 + 𝛽 𝜖 2   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1 , 𝑅   8 = 2 𝜀 𝜖 1  𝜖 2 + 𝛽 𝜖 1   1 s i n h 2  𝜀 𝑘 1   + 𝛽 𝛽 𝜖 2 + 𝜖 1 𝜖   2 − 𝜖 1   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1  , 𝑅 9  1 = t a n h 2  𝜀 𝑘 1   𝛽 + t a n h 2  𝜀 𝑘 1  , 𝑅 1 0 = 𝑘 2 𝑘 1 𝜀 𝑅 2  𝜖 1 + 𝛽 𝜖 2   1 s e c h 2  𝜀 𝑘 1  , 𝑅 1 1 = 𝐵 1  𝛽 √ √ 𝜀 − 2 𝑘 1  𝛽 t a n h 2  𝜀 𝑘 1 , 𝑅   1 2 = 𝑅 2 𝜇 𝑟 √ 𝑘 1  𝛽 t a n h 2  𝜀 𝑘 1  , 𝑅 1 3 = 𝑅 9 𝑅 1 2 + 𝑅 1 1 𝑅 2 , 𝑅 1 4 = 𝑅 2 𝑅 1 2 𝑘 5 1  1 t a n h 2  𝜀 𝑘 1  , 𝑅 1 5 = 𝜀 𝑅 1 0 𝑅 1 2  1 s i n h 2  𝜀 𝑘 1  , 𝑅 1 6 = 𝑅 2  4 𝑅 3 + 𝑞 0 𝜇 𝑟 𝑅 2  𝜀 𝐵 2  𝑆 1 ( 0 ) − 𝑆 2 ( 0 )  + 4 𝛽 𝑘 1 𝑘 4 𝜖 2  1 t a n h 2  𝜀 𝑘 1 , 𝑅    1 7 = 𝑅 2  4 𝑅 3 ( 2 𝛽 − 1 ) + 3 𝛽 q 0 𝜇 𝑟 𝑅 2  𝜀 𝐵 1  𝑆 1 ( 0 ) − 𝑆 2 ( 0 )  + 2 𝛽 𝑘 1 𝑘 4 𝜖 2  1 t a n h 2  𝜀 𝑘 1 , 𝑅    1 8 = 𝛽 𝑞 0 𝑘 1 𝑘 4 𝑅 2 𝜀  1 t a n h 2  𝜀 𝑘 1  , 𝑅 1 9 = 𝑅 2  − 4 𝑅 3 + 𝛽 𝜀 𝑞 0 𝜇 𝑟 𝐵 1 𝑅 2  𝑆 1 ( 0 ) − 𝑆 2 ( 0 ) , 𝑅   2 0 = 𝑅 2  4 ( 𝛽 − 2 ) 𝑅 2 + 3 𝛽 𝜀 𝑞 0 𝜇 𝑟 𝐵 1 𝑅 3  𝑆 1 ( 0 ) − 𝑆 2 ( 0 ) , 𝑅   2 1 = − 𝑘 1 𝑘 2 (  𝜖 1 + 𝛽 ) 1 + 𝛽 𝜖 2  𝑞 2 0 𝜖 2  ( 𝛽 − 1 ) 𝜀 𝜖 1 +  𝜖 1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 2 = 𝑞 0 𝑘 1 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2  𝜀  𝜖 1 + 𝛽 2 𝜖 2  −  𝜖 1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 3 = 𝑞 0 𝜖 2  𝜖 1 + 𝛽 𝜖 2  𝑘 1 𝑘 2  2 𝜀 𝜖 1  𝛽 𝜖 1 + 𝜖 2  +  𝜖 1 − 𝜖 2 𝜖   1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 4 = 𝜀 ( 1 − 𝛽 ) 𝑘 1 𝑘 2 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2  2 , 𝑅 2 5 = 𝜀 𝑘 1 𝑘 2  𝜖 1 + 𝛽 𝜖 2 𝜖   3 1  𝑘 2 − 𝜖 2  + 3 𝛽 𝑘 2 𝜖 2 1 𝜖 2 + 𝛽 3 𝑘 2 𝜖 3 2 + 𝜖 1 𝜖 2 2  3 𝛽 2 𝑘 2 + 𝜖 2 , 𝑅   2 6 = − 𝜀 𝑞 0 𝑘 1 𝑘 2  𝐶 1 𝐵 1  + 2 𝑘 1 𝜀  𝜖 2 1  𝜖 1 + 𝛽 ( 2 + 𝛽 ) 𝜖 2   1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 7 = 𝜀 𝜖 2 𝑘 1 𝑘 2  𝐶 1 𝐵 1  + 2 𝑘 1 𝜀  ×  ( 𝛽 − 1 ) 𝜖 1  𝜖 1 + 𝛽 𝜖 2  2 − 𝛽 2 𝑞 0 𝜖 2  𝛽 2 𝜖 2 + ( 1 + 2 𝛽 ) 𝜖 1    1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 8 = 𝜀 𝑘 1 𝑘 2 𝜖 1 𝜖 2 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀  ×  𝛽  𝛽 2  𝑞 − 1 2 0 𝜖 2 +  𝜖 1 + 𝛽 𝜖 2 𝜖   2 1 − 𝜖 2 2  − 2 𝑞 0 𝜖 2   1 + 𝛽 2  𝜖 1 + 𝛽 𝜖 2   − 4 𝛽 𝑞 0 𝜖 2 1  𝑘 1 𝜀   1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 9 = 𝜀 𝑘 1 𝑘 2 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀  ×  𝑘 2  𝜖 4 1 + 4 𝛽 3 𝜖 1 𝜖 3 2 + 𝛽 4 𝜖 4 2  + 𝜖 2 1 𝜖 2   1 − 𝛽 2  𝑞 2 0 + 2 𝛽 𝑘 2  2 𝜖 1 + 3 𝛽 𝜖 2  𝐶    + 2 𝛽 1 𝐵 1  𝑞 0 𝜖 3 1 𝜖 2   1 t a n h 2  𝜀 𝑘 1  , 𝑅 3 0 = 𝛽 𝑞 0 𝑘 2 𝜖 2  𝜖 1 + 𝛽 𝜖 2   ( 1 + 𝛽 ) 𝑞 0 − 𝜖 1 + 𝜖 2  , 𝑅 3 1 = 𝛽 𝑞 0 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2 , 𝑅 3 2 = 𝑘 2 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2 (   𝛽 − 1 ) 𝑞 0 − 𝜖 1 − 𝜖 2 (   𝛽 + 1 ) 𝑞 0 − 𝜖 1 + 𝜖 2  , 𝑅 3 3 = 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2  ( 𝛽 − 1 ) 𝜖 1 𝜖 2 − 𝑞 0  𝜖 1 + 𝛽 2 𝜖 2 , 𝑅   3 4 = 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑆 1 ( 0 ) − 3 𝛽 𝑆 2 ( 0 )  𝜖 2  , 𝑅 3 5 = 𝛽 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑆 1 ( 0 ) − 𝛽 𝑆 2 ( 0 )  𝜖 2  , 𝑅 3 6 = 𝜀 𝑘 1 𝑅 2  1 t a n h 2  𝜀 𝑘 1  , 𝑅 3 7 = 𝑅 1 4 + 𝑘 4 𝑅 1 𝑅 1 3 , 𝑅 3 8 =  𝑅 1 𝐶 ( 𝛽 − 1 ) − 1 𝐵 1 𝛽  𝑘 1 𝜀  𝜖 1 − 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 , 𝑅 3 9  = 𝛽 2 𝑅 1 − 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 1 +  2 𝑅 1 − 𝐶 1 𝐵 1  𝛽 ( 𝛽 − 1 ) 𝑘 1 𝜀  𝜖 1 𝜖 2 + 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 2 , 𝑅 4 0 =  𝑅 1 𝐶 + 𝛽 1 𝐵 1  𝑘 1 𝜀  𝜖 1 + 𝛽 2  𝑅 1 + 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜀 2 , 𝑅 4 1 = − 𝑅 3 7 𝛽 4 𝜖 4 2 + 𝜖 4 1  𝑅 3 7 − 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  − 𝜖 2 1 𝜖 2 2  6 𝑅 3 7 𝛽 2 + 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  + 𝜖 3 1 𝜖 2  𝑅 1 5 − 4 𝛽 𝑅 3 7 + 𝛽 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  − 𝜖 1 𝜖 3 2  𝑅 1 5 + 4 𝛽 3 𝑅 3 7 + 𝛽 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  + ( 1 + 𝛽 ) 𝑞 2 0 𝜖 2  − 𝑅 6 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3 𝑅 3 8  𝜖 1 + 𝛽 𝜖 2   − 𝑞 0 𝜖 2  𝑅 8 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3 𝑅 3 9  𝜖 1 + 𝛽 𝜖 2 , 𝑅   4 2 = 𝛽 𝜖 2  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2 ×  𝑅 1 9 + 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑅 1 8 + 𝐵 1 𝑅 4 𝑆 1 ( 0 )  𝜖 2  𝜇 𝑟   , 𝑅 4 3 = 𝜖 1  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2 ×  𝑅 2 0 + 3 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  3 𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  4 𝛽 𝑅 1 8 + 3 𝐵 1 𝑅 4 𝑆 1 ( 0 ) − 𝛽 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 2  𝜇 𝑟   , 𝑅 4 4 = 𝜀 𝑅 2  𝑅 1 8 − 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 4 1 𝜇 𝑟 + 𝜖 3 1  𝑅 1 6 + 𝑅 3 𝑅 4 + 𝐵 1 𝑅 2 𝑅 4 𝑅 3 4  + 𝜖 2 1 𝜖 2  𝑅 1 7 + 3 𝛽 𝑅 3 𝑅 4 + 3 𝐵 1 𝑅 2 𝑅 4 𝑅 3 5  , 𝑅 4 5 = − 𝑅 2  𝜖 1 + 𝛽 𝜖 2  𝜇 𝑟 𝑞 0  𝑅 9  𝑅 2 2 + 𝑅 2 4 + 𝑅 2 6 + 𝑅 2 7  − 𝑅 3 1 𝑅 2 √ 𝜀 𝑘 1  𝜖 1 + 𝛽 𝜖 2  + 𝑅 3 3 𝑅 3 6  , 𝑅 4 6 = 𝐵 1 𝑅 2 2 𝜀 2 𝜇 𝑟  𝜖 1 + 𝛽 𝜖 2  4 , 𝑅 4 7 = 𝑅 4 6  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  + 𝑅 4 5 , 𝑅 4 8 = 𝑅 4 1 𝑅 4 6  𝜖 1 + 𝛽 𝜖 2  + 𝐵 1 𝑅 5 𝑅 4 7 𝜇 𝑟 𝜀 3 / 2 , 𝑅 4 9 = 2 𝐵 1 𝑅 5 𝑅 4 6 𝜇 𝑟 𝜀 3 / 2  𝜖 1 + 𝛽 𝜖 2  . ( A . 1 ) <h4>References</h4> G. I. Taylor and A. D. McEwan, “ The stability of a horizontal fluid interface in a vertical electric field ,” The Journal of Fluid Mechanics , vol. 22, pp. 1–15, 1965. J. R. Melcher and C. V. Smith, “ Electrohydrodynamic charge relaxation and interfacial perpendicular-field instability ,” Physics of Fluids , vol. 12, no. 4, pp. 778–790, 1969. D. A. Saville, “ Electrohydrodynamics: the Taylor-Melcher leaky dielectric model ,” Annual Review of Fluid Mechanics , vol. 29, pp. 27–64, 1997. D. J. Griffiths, Introduction to Electrohydrodynamics , Pearson Education, Delhi, India, 3rd edition, 2006. A. A. Mohamed, E. F. Elshehawey, and M. F. El-Sayed, “ Electrohydrodynamic stability of two superposed viscous fluids ,” Journal of Colloid And Interface Science , vol. 169, no. 1, pp. 65–78, 1995. N. T. M. El-Dabe, “ Electrohydrodynamic stability of two superposed elasticoviscous liquids in plane Couette flow ,” Journal of Mathematical Physics , vol. 28, no. 11, pp. 2791–2800, 1987. D. D. Joseph and Y. Y. Renardy, Fundamentals of Two-Fluid Dynamics—Part I: Mathematical Theory and Application , vol. 3 of Interdisciplinary Applied Mathematics , Springer, New York, NY, USA, 1993. K. Abdella and H. Rasmussen, “ Electrohydrodynamic instability of two superposed fluids in normal electric fields ,” Journal of Computational and Applied Mathematics , vol. 78, no. 1, pp. 33–61, 1997. J. R. Melcher, Continuum Electromechanics , MIT Press, Cambridge, UK, 1981. E. Schaffer, T. Thurn-Albrecht, T. P. Russell, and U. Steiner, “ Electrohydrodynamic instabilities in polymer films ,” Europhysics Letters , vol. 53, no. 4, pp. 518–523, 2001. Z. Lin, T. Kerle, S. M. Baker, et al., “ Electric field induced instabilities at liquid/liquid interfaces ,” Journal of Chemical Physics , vol. 114, no. 5, pp. 2377–2381, 2001. M. D. Morariu, N. E. Voicu, E. Schaffer, Z. Lin, T. P. Russell, and U. Steiner, “ Hierarchical structure formation and pattern replication induced by an electric field ,” Nature Materials , vol. 2, no. 1, pp. 48–52, 2003. O. Ozen, N. Aubry, D. T. Papageorgiou, and P. G. Petropoulos, “ Electrohydrodynamic linear stability of a two immiscible fluids in channel flow ,” Electrochimica Acta , vol. 51, pp. 5316–5323, 2006. L. Wu and S. Y. Chou, “ Electrohydrodynamic instability of a thin film of viscoelastic polymer underneath a lithographically manufactured mask ,” Journal of Non-Newtonian Fluid Mechanics , vol. 125, no. 2-3, pp. 91–99, 2005. R. G. Larson, Constitutive Equations for Polymer Melts and Solutions , Butterworth, Stoneham, Mass, USA, 1988. R. Verma, A. Sharma, K. Kargupta, and J. Bhaumik, “ Electric field induced instability and pattern formation in thin liquid films ,” Langmuir , vol. 21, no. 8, pp. 3710–3721, 2005. D. T. Papageorgiou and P. G. Petropoulos, “ Generation of interfacial instabilities in charged electrified viscous liquid films ,” Journal of Engineering Mathematics , vol. 50, no. 2-3, pp. 223–240, 2004. V. Shankar and A. Sharma, “ Instability of the interface between thin fluid films subjected to electric fields ,” Journal of Colloid and Interface Science , vol. 274, pp. 294–308, 2004. R. V. Craster and O. K. Matar, “ Electrically induced pattern formation in thin leaky dielectric films ,” Physics of Fluids , vol. 17, no. 3, Article ID 032104, 2005. F. Li, O. Ozen, N. Aubry, D. T. Papageorgiou, and P. G. Petropoulos, “ Linear stability of a two-fluid interface for electrohydrodynamic mixing in a channel ,” Journal of Fluid Mechanics , vol. 583, pp. 347–377, 2007. G. Tomar, V. Shankar, A. Sharma, and G. Biswas, “ Electrohydrodynamic instability of a confined viscoelastic liquid film ,” Journal of Non-Newtonian Fluid Mechanics , vol. 143, no. 2-3, pp. 120–130, 2007. G. Supeene, C. R. Koch, and S. Bhattacharjee, “ Deformation of a droplet in an electric field: nonlinear transient response in perfect and leaky dielectric media ,” Journal of Colloid and Interface Science , vol. 318, no. 2, pp. 463–476, 2008. S. Herminghaus, “Dynamical instability of thin liquid films between conducting media,” Physical Review Letters , vol. 83, pp. 2539–2542, 1999. L. F. Pease III and W. B. Russel, “ Linear stability analysis of thin leaky dielectric films subjected to electric fields ,” Journal of Non-Newtonian Fluid Mechanics , vol. 102, no. 2, pp. 233–250, 2002. D. B. Ingham and I. Pop, Transport Phenomena in Porous Media , Pergamon Press, Oxford UK, 1998. K. Vafai, Handbook of Porous Media , Marcel Dekker, New York, NY, USA, 2000. J. A. Del Rio and S. Whitaker, “ Electrohydrodynamics in porous media ,” Transport in Porous Media , vol. 44, no. 2, pp. 385–405, 2001. I. Pop and D. B. Ingham, Convevtive Heat Transfer: Mathematical and Computational Modeling of Viscous Fluids and Porous Media , Pergamon Press, Oxford, UK, 2001. D. A. Nield and A. Bejan, Convection in Porous Media , Springer, New York, NY, USA, 3rd edition, 2006. M. F. El-Sayed, “ Electrohydrodynamic instability of two superposed viscous streaming fluids through porous media ,” Canadian Journal of Physics , vol. 75, no. 7, pp. 499–508, 1997. M. F. El-Sayed, “ Effect of normal electric fields on Kelvin-Helmholtz instability for porous media with Darcian and Forchheimer flows ,” Physica A , vol. 255, no. 1-2, pp. 1–14, 1998. M. F. El-Sayed, “ Electrohydrodynamic instability of dielectric fluid layer between two semi-infinite identical conducting fluids in porous medium ,” Physica A , vol. 367, pp. 25–41, 2006. M. F. El-Sayed, “ Instability of two streaming conducting and dielectric bounded fluids in porous medium under time-varying electric field ,” Archive of Applied Mechanics , vol. 79, no. 1, pp. 19–39, 2009. M. F. El-Sayed, “ Thermohydrodynamic instability of viscous rotating dielectric fluid layer in porous medium with vertical ac electric field ,” Special Topics & Reviews in Porous Media , vol. 1, pp. 15–29, 2010. M. F. El-Sayed, “ Electrohydrodynamic instability conditions for two superposed dielectric bounded fluids streaming with fine dust in a porous medium ,” Journal of Porous Media , vol. 13, pp. 457–473, 2010. M. F. El-Sayed, G. M. Moatimid, and T. M. N. Metwaly, “Nonlinear electrohydrodynamic stability of two superposed streaming finite dielectric fluids in porous medium with interfacial surface charges,” Transport in Porous Media , vol. 86, pp. 559–578, 2011. P. Gambhire and R. M. Thaokar, “Electrohydrodynamic instabilities at interfaces subjected to alternating electric field,” Physics of Fluids , vol. 22, Article ID 064103, 2010. // http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png ISRN Applied Mathematics Hindawi Publishing Corporation

Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium

Loading next page...
 
/lp/hindawi-publishing-corporation/electrohydrodynamic-instability-of-two-thin-viscous-leaky-dielectric-lpYo2NP9mR

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2011 M. F. El-Sayed et al.
ISSN
2090-5564
eISSN
2090-5572
Publisher site
See Article on Publisher Site

Abstract

Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium //// Advanced Search home journals subjects indexing guidelines submit resources about Journal Menu About this Journal Abstracting and Indexing Aims and Scope Article Processing Charges Articles in Press Bibliographic Information Contact Information Editorial Board Editorial Workflow Free eTOC Alerts Table of Contents Abstract Full-Text PDF Full-Text HTML Full-Text ePUB Linked References How to Cite this Article ISRN Applied Mathematics Volume 2011 (2011), Article ID 498718, 35 pages doi:10.5402/2011/498718 Review Article <h2>Electrohydrodynamic Instability of Two Thin Viscous Leaky Dielectric Fluid Films in a Porous Medium</h2> M. F. El-Sayed , 1,2 M. H. M. Moussa , 1 A. A. A. Hassan , 1 and N. M. Hafez 1 1 Department of Mathematics, Faculty of Education, Ain Shams University, Heliopolis (Roxy), Cairo 11757, Egypt 2 Department of Mathematics, College of Science, Qassim University, P.O. Box 6644, Buraidah 51452, Saudi Arabia Received 23 March 2011; Accepted 4 May 2011 Academic Editor: R. Cardoso Copyright © 2011 M. F. El-Sayed et al. This is an open access article distributed under the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract The effect of an applied electric field on the stability of the interface between two thin viscous leaky dielectric fluid films in porous medium is analyzed in the long-wave limit. A systematic asymptotic expansion is employed to derive coupled nonlinear evolution equations of the interface and interfacial free charge distribution. The linearized stability of these equations is determined and the effects of various parameters are examined in detail. For perfect-perfect dielectrics, the various parameters affect only for small wavenumber values. For dielectrics, the various parameters affect only for small wavenumber values. For effect for small wavenumbers, and a stabilizing effect afterwards, and for high wavenumber values for the other physical parameters, new regions of stability or instability appear. For leaky-leaky dielectrics, the conductivity of upper fluid has a destabilizing effect for small or high wavenumbers, while it has a dual role on the stability of the system in a wavenumber range between them. The effects of all other physical parameters behave in the same manner as in the case of perfect-leaky dielectrics, except that in the later case, the stability or instability regions occur more faster than the corresponding case of leaky-leaky dielectrics. 1. Introduction The effect of electric fields on the stability and dynamics of fluid-fluid interfaces has been an area of extensive research, beginning from the classic works of Taylor and McEwan [ 1 ] and Melcher and Smith [ 2 ]. These works and subsequent studies, see, for example, the reviews of Saville [ 3 ] and Griffiths [ 4 ], have amply demonstrated the role of electrical stresses on fluid interfaces and the associated electrohydrodynamic instabilities in such systems. One of the basic problems here is to understand the stability of the interface between two fluid layers bounded on the top and bottom by rigid plates, and this has been the subject of many previous studies. Mohamed et al. [ 5 ] concentrated on two superposed viscous fluids in a channel subjected to a normal electric field, where the upper fluid is highly conducting, while the lower fluid is dielectric, and they performed the long-wave linear stability analysis, [ 6 , 7 ], and showed that the electric field always has a destabilizing effect on the flow. Abdella and Rasmussen [ 8 ] studied Couette flow of two viscous fluids with different viscosities, densities, conductivities and permittivities, in an unbounded domain subjected to a normal electric field. They studied, following Melcher [ 9 ], two special cases in detail: the electrohydrodynamic free-charge configuration (EH-If) and the electrohydrodynamic polarization charge configuration (EH-If). These studies have largely considered systems in which gravitational effects are important, and therefore, a critical applied voltage is required to cause the instability, very long waves are stabilized by interfacial tension, and waves of intermediate lengths become unstable. These earlier studies have focused on how the critical voltage required for instability is affected by the nature of the fluids, namely. whether they are perfect dielectrics, or whether they are leaky dielectrics in which there is the possibility of free charge conduction in the fluids, and also the possibility of accumulation and redistribution of charges on the interface between two fluids. Recently, there has been a renewed interest in this area, in part due to the relevance of such phenomena in the formation of well-controlled patterns using the application of electric fields to thin liquid films [ 10 – 12 ], which has demonstrated that the application of an external electric field to polymer-air or polymer-polymer interfaces enhances the spontaneous fluctuations at the interface leading to an instability. However, for leaky dielectrics with surface charges at the interface and two fluids that are not perfectly conducting, we must bear in mind that there is an electrical tangential shear stress at the interface, induced by the electric field, and hence, it changes the stability of the two-fluid layer system as previously investigated by Ozen et al. [ 13 ]. Wu and Chou [ 14 ] performed linear stability analysis of a leaky dielectric viscoelastic fluid whose constitutive behavior was described using the Jeffrey's model [ 15 ]. The surface instability of a Newtonian fluid (modeled as both leaky and perfect dielectrics) under the effect of electric field is now well understood [ 16 ]. For recent developments of this topic, see the investigations of Papageorgiou and Petropoulos [ 17 ], Shankar and Sharma [ 18 ], Craster and Matar [ 19 ], Li et al. [ 20 ], Tomar et al. [ 21 ], and Supeene et al. [ 22 ]. Theoretical studies that have considered the linear stability characteristics of thin fluid films subjected to electric fields were restricted to the following configurations: (i) the interface between a perfect dielectric liquid and air [ 10 ], (ii) the interface between two perfect dielectric liquids [ 11 ], (iii) the interface between a dielectric fluid of finite thickness and a conducting fluid of much larger thickness [ 23 ], and (iv) the interface between a leaky dielectric liquid and air [ 24 ]. Except the study of Pease III and Russel [ 24 ], all these studies have focused only on perfect dielectric systems. The presence of conductivity in one or both of the liquids could have a significant impact on the length scales and growth rates. It should be noted that in all the above-mentioned studies, the medium has been considered to be nonporous. Porous media theories, on one hand, play an important role in many branches of engineering, including material science, petroleum industry, chemical engineering, and soil mechanics as well as biomechanics. The flow through porous media has gained considerable interest in recent years, particularly among geophysical fluid dynamicists. It is well known that in Darcy's law, which relates the pressure gradient, the bulk viscous resistance, and the gravitational force in a porous medium, the usual viscous term in the equation of motion is replaced by the resistive force − ( 𝜇 / 𝑘 1 ) 𝐯 , where 𝜇 is the fluid viscosity, 𝑘 1 is the medium permeability, and 𝐯 is the Darcian velocity of the fluid. Much of the recent studies on this topic are given by Ingham and Pop [ 25 ], Vafai [ 26 ], Del Rio and Whitaker [ 27 ], Pop and Ingham [ 28 ], and Nield and Bejan [ 29 ]. On the other hand, electrohydrodynamic instability studies for flows in porous media has attracted little attention in the scientific literature [ 30 – 36 ] despsite their applications in various diverse fields with great interest. Thus, there is a growing need for original research in the updated electrohydrodynamic phenomena which have some physical and engineering applications. In this study, we consider the most general case of the effect of an external applied electric field on the stability of the interface between two thin leaky dielectric fluids with arbitrary viscosities and conductivities. We first use a systematic long-wave asymptotic analysis to derive the nonlinear evolution equations for the interface position and interfacial charge distribution, and we subsequently study the linearized stability of the nonlinear differential equations. This paper is organized in the following manner: Section 2 discusses the relevant governing equations and boundary conditions, and in Section 3 , we nondimensionalize these equations and conditions. In Sections 4 and 5 , we outline the long-wave asymptotic analysis used to derive the nonlinear evolution equations for the interface position and charge, and in Section 6 , we develop the linear stability analysis of the nonlinear equations, and discuss the important representative studies results obtained from the general case of two leaky dielectric interface and the limiting cases of both perfect-leaky dielectric, and perfect-perfect dielectric interfaces. Finally, the salient conclusions of the present study are discussed in Section 7 . 2. Problem Formulation and Governing Equations The system of interest consists of two leaky dielectric fluids in porous medium of arbitrary viscosities occupying the regions − 𝐻 < 𝑦 < 0 (fluid 2) and 0 < 𝑦 < 𝛽 𝐻 (fluid 1) in the initial unperturbed state, see Figure 1 , where 𝛽 is the ratio of the thicknesses of top and bottom fluids. The perturbed interface between the two fluids is denoted by 𝑦 = ℎ ( 𝑥 ) . The two fluids are stationary in the initial state, with viscosities 𝜇 𝑖 , dielectric constants 𝜖 𝑖 , conductivities 𝜎 𝑖 , porosity of porous medium 𝜀 and medium permeability 𝑘 1 , where 𝑖 = 1 , 2 . Fluid 2 is bounded at the bottom 𝑦 = − 𝐻 by a rigid plate which is maintained at an electric potential 𝜓 = 𝜓 𝑏 , while fluid 1 is bounded at the top 𝑦 = 𝛽 𝐻 by a rigid plate maintained at an electric potential 𝜓 = 0 . In the ensuing analysis, we assume that the material properties of the fluid such as viscosities 𝜇 𝑖 , dielectric constants 𝜖 𝑖 , conductivities 𝜎 𝑖 , porosity of porous medium 𝜀 , and medium permeability 𝑘 1 are constants and independent of spatial position. Following the leaky dielectric model formulation of Saville [ 3 ], we assume that electroneutrality is valid in the bulk, while free charge is assumed to accumulate at the fluids interface. We also neglect the diffusion of free charge within the interface. Figure 1: Schematic diagram showing the configuration and coordinate system. For leaky dielectric fluids of constant conductivities and with zero net charge in the bulk, the following governing equations are appropriate for the electric field 𝐄 in the two fluids 1 and 2 [ 3 ] ∇ ⋅ 𝐄 𝑖 = 0 , 𝑖 = 1 , 2 . ( 2 . 1 ) Since the electric fields are irrotational, 𝐄 𝑖 = − ∇ 𝜓 𝑖 , where 𝜓 𝑖 is the electric potential: in the fluid 𝑖 . Substituting this into ( 2.1 ) gives the following Laplace's equation for the electric potential: 𝜓 𝑖 ∇ 2 𝜓 𝑖 = 0 . ( 2 . 2 ) These governing equations are supplemented by the following boundary conditions. (1) The normal component of the electric field, at the interface 𝑦 = ℎ ( 𝑥 ) , satisfies 𝜖 2 𝜖 0  ∇ 𝜓 2  ⋅ 𝐧 − 𝜖 1 𝜖 0  ∇ 𝜓 1  ⋅ 𝐧 = 𝑞 ( 𝑥 , 𝑡 ) a t 𝑦 = ℎ ( 𝑥 ) , ( 2 . 3 ) where 𝐧 is the unit normal to the interface at 𝑦 = ℎ ( 𝑥 ) (see Figure 1 ), 𝜖 𝑖 ( 𝑖 = 1 , 2 ) is the dielectric constant in fluid 𝑖 , 𝜖 0 is the permittivity of free space, and 𝑞 ( 𝑥 , 𝑡 ) is the surface charge density of free charges at the interface. (2) The continuity of the tangential component of the electric field at the interface 𝑦 = ℎ ( 𝑥 ) translates to the continuity of the electric potentials; that is, 𝜓 1 = 𝜓 2 a t 𝑦 = ℎ ( 𝑥 ) . ( 2 . 4 ) (3) The electric potentials satisfy the following conditions at the rigid boundaries: 𝜓 2 = 𝜓 𝑏 𝜓 a t 𝑦 = − 𝐻 , 1 = 0 a t 𝑦 = 𝛽 𝐻 . ( 2 . 5 ) We next turn to the equations governing the motion of the two fluids. Owing to the relatively small thicknesses of the fluids, we ignore inertial effects in both fluids, and hence, the governing equations are the Stokes equations for continuity and momentum balance ∇ ⋅ 𝐕 𝑖 = 0 , ∇ ⋅ 𝐓 𝑖 = 0 , ( 2 . 6 ) where 𝐕 𝑖 and 𝐓 𝑖 are the velocity field and the total stress tensor, respectively, in fluid 𝑖 . Both the fluids are assumed to be irrotational, then the fluid velocity 𝐕 𝑖 can be derived from a scalar velocity potential 𝜑 𝑖 such that 𝐕 𝑖 = − ∇ 𝜑 𝑖 . We have neglected the effects of gravity on the length scales of interest here. In addition, the effects of van der Waals dispersion forces are negligible for the films considered in the experimental studies [ 11 ]. The total stress tensor 𝐓 is given by a sum of isotropic pressure, deviatoric viscous stresses for the Newtonian fluid, the electrical Maxwell stress tensor, and a Darcy's law term describing the isotropic porous medium 𝐓 𝑖 = − 𝑝 𝑖 𝐈 + 𝜇 𝑖  ∇ 𝐕 𝑖 + ∇ 𝐕 𝑇 𝑖  + 𝐦 𝑖 + 𝜇 𝑖 𝑘 1 𝜑 𝑖 𝐈 , 𝑖 = 1 , 2 , ( 2 . 7 ) where 𝑝 𝑖 is the pressure in fluid 𝑖 , 𝐈 is the identity tensor, and the Maxwell stress tensor 𝐦 𝑖 is given by Saville [ 3 ] 𝐦 𝑖 = 𝜖 𝑖 𝜖 0  𝐄 𝑖 𝐄 𝑖 − 1 2  𝐄 𝑖 ⋅ 𝐄 𝑖  𝐈  . ( 2 . 8 ) The divergence of the Maxwell stress tensor ∇ ⋅ 𝐦 𝑖 = 0 , because the bulk of the fluid is free of net charge, and the dielectric constants are independent of spatial position in the two fluids; that is, ∇ ⋅ 𝐦 𝑖  𝜖 = ∇ ⋅ 𝑖 𝜖 0  𝐄 𝑖 𝐄 𝑖 − 1 2  𝐄 𝑖 ⋅ 𝐄 𝑖  𝐈 𝜖   = − 0 𝐄 𝑖 2 ⋅ 𝐄 𝑖 ∇ 𝜖 𝑖 + 𝜌 𝑓 𝐄 𝑖 = 0 , ( 2 . 9 ) with 𝜌 𝑓 is the bulk free charge. Thus, the Maxwell stress tensor will not appear in the momentum balance but will affect the flow only through the conditions at the interface. The governing momentum equations in the two fluids, therefore, become 𝜌 𝜀 𝐷 𝐕 𝑖 𝐷 𝑡 = − ∇ 𝑝 𝑖 + ∇ ⋅ 𝐦 𝑖 + 𝜇 𝑖 𝜀 ∇ 2 𝐕 𝑖 − 𝜇 𝑖 𝑘 1 𝐕 𝑖 , 𝑖 = 1 , 2 . ( 2 . 1 0 ) Since there is no time scale, then 𝐷 𝐯 𝑖 / 𝐷 𝑡 = 0 , and we get from the above two equations that ∇ 𝑝 𝑖 − 𝜇 𝑖 𝜀 ∇ 2 𝐕 𝑖 + 𝜇 𝑖 𝑘 1 𝐕 𝑖 = 0 , 𝑖 = 1 , 2 . ( 2 . 1 1 ) The fluid velocities satisfy no-slip and no-penetration conditions at the top and bottom plates 𝐕 1 ( 𝑦 = 𝛽 𝐻 ) = 0 , 𝐕 2 ( 𝑦 = − 𝐻 ) = 0 . ( 2 . 1 2 ) At the interface 𝑦 = ℎ ( 𝑥 ) between the two fluids, continuity of velocities and stresses apply ( 𝐕 ⋅ 𝐧 ) 1 = ( 𝐕 ⋅ 𝐧 ) 2 , ( ( 2 . 1 3 ) 𝐕 ⋅ 𝐭 ) 1 = ( 𝐕 ⋅ 𝐭 ) 2 , ( 2 . 1 4 ) ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 2 = ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 1 + 𝛾 𝜅 , ( 2 . 1 5 ) ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 2 = ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 1 , ( 2 . 1 6 ) where 𝛾 is the interfacial tension between the two fluids, 𝐭 is the unit tangent to the interface (see Figure 1 ), and 𝜅 = 𝜕 2 ℎ ( 𝑥 ) / 𝜕 𝑥 2 is the mean curvature of the interface. We restrict our attention to two-dimensional systems which are invariant in the 𝑧 direction and denote the velocity components in the 𝑥 and 𝑦 directions by 𝑢 and 𝜐 , respectively; that is, 𝐕 𝑖 = ( 𝑢 𝑖 , 𝜐 𝑖 ) . Upon substituting the expression for the Maxwell stress tensor ( 2.8 ) in ( 2.15 ) and ( 2.16 ), respectively, we get ( 𝐧 ⋅ 𝐓 ⋅ 𝐧 ) 𝑖 = − 𝑝 𝑖 + 2 𝜇 𝑖 𝜕 𝜐 𝑖 𝜕 𝑦 + 𝜖 𝑖 𝜖 0  𝐸 𝑁 2 𝑖 − 1 2 𝐸 2 𝑖  + 𝜇 𝑖 𝑘 1 𝜑 𝑖 , ( 𝐭 ⋅ 𝐓 ⋅ 𝐧 ) 𝑖 = 𝜇 𝑖  𝜕 𝑢 𝑖 + 𝜕 𝑦 𝜕 𝜐 𝑖  𝜕 𝑥 + 𝜖 𝑖 𝜖 0 𝐸 𝑇 𝑖 𝐸 𝑁 𝑖 , ( 2 . 1 7 ) where 𝐸 𝑇 𝑖 and 𝐸 𝑁 𝑖 are the tangential and normal components of the electric field 𝐄 𝑖 in the fluid 𝑖 , respectively. Hence, we obtain the following conditions to be applied at the interface 𝑦 = ℎ ( 𝑥 ) :  𝑝 1 − 2 𝜇 1 𝜕 𝜐 1  −  𝑝 𝜕 𝑦 2 − 2 𝜇 2 𝜕 𝜐 2  + 1 𝜕 𝑦 𝑘 1  𝜇 2 𝜑 2 − 𝜇 1 𝜑 1  𝜖 = 𝛾 𝜅 + 1 𝜖 0 2   ∇ 𝜓 1  ⋅ 𝐧 2 −  ∇ 𝜓 1  ⋅ 𝐭 2  − 𝜖 2 𝜖 0 2   ∇ 𝜓 2  ⋅ 𝐧 2 −  ∇ 𝜓 2  ⋅ 𝐭 2  , ( 2 . 1 8 ) for the normal stress continuity, and 𝜇 2  𝜕 𝑢 2 + 𝜕 𝑦 𝜕 𝜐 2  𝜕 𝑥 − 𝜇 1  𝜕 𝑢 1 + 𝜕 𝑦 𝜕 𝜐 1   𝜕 𝑥 = − ∇ 𝜓 1  ⋅ 𝐭 𝑞 ( 𝑥 , 𝑡 ) , ( 2 . 1 9 ) for the tangential stress continuity. In the last equation, we have used the normal electric field continuity condition, ( 2.3 ), to simplify the right-hand side. The kinematic condition at the interface prescribes the evolution of the interface position ℎ ( 𝑥 , 𝑡 ) , 𝜐 1 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜐 2 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 ℎ , ( 2 . 2 0 ) where ∇ 𝑠 is the gradient operator along the interface 𝑦 = ℎ ( 𝑥 ) . Finally, the interfacial charge is governed by a conservation equation 𝜀 𝜕 𝑞 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 𝜎 𝑞 − 𝑞 𝐧 ⋅ ( 𝐧 ⋅ ∇ ) 𝐕 = 𝜀   2 𝐄 2  −  𝜎 ⋅ 𝐧 1 𝐄 1 ⋅ 𝐧   , ( 2 . 2 1 ) where the terms on the left-hand side represent, respectively, the accumulation, convection, and variation of the charge due to dilation of the interface, while the right side represents the migration of charge to or from the interface due to ion conduction in the bulk [ 3 ]. 3. Nondimensional Forms It is useful at this point to nondimensionalize the governing equations and boundary conditions by setting 𝜓  𝑖 = 𝜓 𝑖 𝜓 𝑏 , 𝑝  𝑖 = 𝑝 𝑖 𝐻 2 𝜖 0 𝜓 2 𝑏 , 𝑇  𝑖 = 𝑇 𝑖 𝐻 2 𝜖 0 𝜓 2 𝑏 , 𝑥  = 𝑥 𝐻 , 𝑦  = 𝑦 𝐻 , 𝐄  𝑖 = 𝐄 𝑖 𝐻 𝜓 𝑏 , 𝐕  𝑖 = 𝐕 𝑖 𝜇 2 𝐻 𝜖 0 𝜓 2 𝑏 , 𝑡  = 𝑡 𝜖 0 𝜓 2 𝑏 𝜇 2 𝐻 2 , 𝜑  𝑖 = 𝜇 2 𝜑 𝑖 𝜖 0 𝜓 2 𝑏 , 𝑞  = 𝑞 𝐻 𝜖 0 𝜓 𝑏 , 𝜅  = 𝜅 𝐻 , 𝑘  1 = 𝑘 1 𝐻 2 , ℎ  = ℎ 𝐻 . ( 3 . 1 ) Upon using these scales to nondimensionalize the above governing equations and boundary conditions, we end up with the following nondimensional set of equations. Without loss of clarity, and for the sake of brevity, we represent nondimensional variables with the same notation in the ensuing discussion by dropping dashes. The nondimensional governing equations for the electric potentials 𝜓 𝑖 are ∇ 2 𝜓 𝑖 = 0 , 𝑖 = 1 , 2 , ( 3 . 2 ) with the following boundary conditions at the interface 𝑦 = ℎ ( 𝑥 ) ,  𝜖 2 ∇ 𝜓 2  −  𝜖 ⋅ 𝐧 1 ∇ 𝜓 1  ⋅ 𝐧 = 𝑞 , 𝜓 1 = 𝜓 2 , ( 3 . 3 ) and the following boundary conditions at the top and bottom boundaries 𝜓 1 ( 𝑦 = 𝛽 ) = 0 , 𝜓 2 ( 𝑦 = − 1 ) = 1 . ( 3 . 4 ) Similarly, the nondimensional equations governing the fluid motion are ∇ ⋅ 𝐕 1 = 0 , ∇ ⋅ 𝐕 2 = 0 , ( 3 . 5 ) ∇ 𝑝 1 − 𝜇 𝑟 𝜀 ∇ 2 𝐕 1 + 𝜇 𝑟 𝑘 1 𝐕 1 = 0 , ∇ 𝑝 2 − 1 𝜀 ∇ 2 𝐕 2 + 1 𝑘 1 𝐕 2 = 0 . ( 3 . 6 ) with 𝜇 𝑟 = 𝜇 1 / 𝜇 2 being the ratio of viscosities of the two fluids. The nondimensional normal and tangential stress continuity conditions at the interface become  𝑝 1 − 2 𝜇 𝑟 𝜕 𝜐 1  −  𝑝 𝜕 𝑦 2 − 2 𝜕 𝜐 2  + 1 𝜕 𝑦 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝜖 𝛾 𝜅 + 1 2   ∇ 𝜓 1  ⋅ 𝐧 2 −  ∇ 𝜓 1  ⋅ 𝐭 2  − 𝜖 2 2   ∇ 𝜓 2  ⋅ 𝐧 2 −  ∇ 𝜓 2  ⋅ 𝐭 2  ,  ( 3 . 7 ) 𝜕 𝑢 2 + 𝜕 𝑦 𝜕 𝜐 2  𝜕 𝑥 − 𝜇 𝑟  𝜕 𝑢 1 + 𝜕 𝑦 𝜕 𝜐 1   𝜕 𝑥 = − ∇ 𝜓 1  ⋅ 𝐭 𝑞 ( 𝑥 , 𝑡 ) , ( 3 . 8 ) where 𝛾 = 𝛾 𝐻 / 𝜖 0 𝜓 2 𝑏 is the nondimensional interfacial tension. The boundary conditions for the velocities at the top and bottom plates become 𝑢 1 𝜐 ( 𝑦 = 𝛽 ) = 0 , 1 𝑢 ( 𝑦 = 𝛽 ) = 0 , 2 ( 𝑦 = − 1 ) = 0 , 𝜐 2 ( 𝑦 = − 1 ) = 0 . ( 3 . 9 ) The nondimensional kinematic condition at the interface is 𝜐 1 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜐 2 ( 𝑦 = ℎ ( 𝑥 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 ℎ , ( 3 . 1 0 ) and the nondimensional charge conservation equation at the interface is 𝜀 𝜕 𝑞 𝜕 𝑡 + 𝐕 ⋅ ∇ 𝑠 𝑆 𝑞 − 𝑞 𝐧 ⋅ ( 𝐧 ⋅ ∇ ) 𝐕 = 𝜀   1 ∇ 𝜓 1  −  𝑆 ⋅ 𝐧 2 ∇ 𝜓 2 ⋅ 𝐧   , ( 3 . 1 1 ) where 𝑆 𝑖 = 𝜎 𝑖 𝜇 2 𝐻 2 / 𝜖 2 0 𝜓 2 𝑏 , 𝑖 = 1 , 2 are the nondimensional conductivities in the two fluids. This completes the specification of the governing equations and boundary conditions, which are highly coupled. Due to the negligible effect of gravity at length scales of interest here, the above system of equations undergoes a long-wave instability. We now carry out a long-wave asymptotic analysis to make the above system of equations tractable, and thereby derive coupled nonlinear evolution equations for the interface position ℎ ( 𝑥 , 𝑡 ) and charge 𝑞 ( 𝑥 , 𝑡 ) . While the main focus of this paper is to analyze the stability of the linearized equations, it is nonetheless useful to first derive the nonlinear evolution equations, since these equations can be used in future studies to understand (by numerical simulations) the nonlinear evolution processes that occur after the linear instability. 4. Long-Wave Asymptotic Analysis In the long-wave limit, the wavelength 𝐿 of the fastest growing modes is much larger than the transverse length scale 𝐻 in the system, and it is useful to define a small parameter 𝛿 = 𝐻 / 𝐿 ≪ 1 . The lateral length scale 𝐿 is determined self-consistently in the following analysis to be  𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 , and this is further estimated below to be much larger than 𝐻 . Similarly, a slow time scale is necessary to describe the dynamics of the interface motion at such large length scales, and this is introduced a little later in ( 4.15 ). In the limit 𝐻 ≪ 𝐿 , the derivatives in the 𝑥 direction should be scaled with 𝐿 . To this end, we define the slowly varying scale 𝜒 in the following manner: 𝜕 𝜕 𝜕 𝑥 = 𝛿 𝜕 𝜒 , ( 4 . 1 ) and 𝜕 / 𝜕 𝜒 ∼ 𝑂 ( 1 ) . When we apply the above scalings, ( 4.1 ), to the Laplace equation, ( 3.2 ), for the electric potential 𝜓 𝑖 , it simplifies in the limit 𝛿 ≪ 1 to 𝜕 2 𝜓 𝑖 𝜕 𝑦 2 = 0 , 𝑖 = 1 , 2 . ( 4 . 2 ) The continuity condition for the normal component of the electric field at the interface 𝑦 = ℎ ( 𝜒 ) , ( 3.3 ), is similarly simplified in the long-wave limit as 𝜖 2 𝜕 𝜓 2 𝜕 𝑦 − 𝜖 1 𝜕 𝜓 1 𝜕 𝑦 = 𝑞 , ( 4 . 3 ) while the other interface condition (second equation in ( 3.3 )) and the boundary conditions, ( 3.4 ), remain unchanged in the long-wave limit. We now turn to the simplification of the momentum equations ( 3.6 ) for the fluid motion in the long-wave limit. It is useful to define the variable 𝜇 𝑟 , 𝑖 such that 𝜇 𝑟 , 𝑖 = 𝜇 𝑟 for 𝑖 = 1 and 𝜇 𝑟 , 𝑖 = 1 for 𝑖 = 2 . The 𝑥 -momentum equation can be simplified in the long-wave limit as 𝛿 𝜕 𝑝 𝑖 − 𝜇 𝜕 𝜒 𝑟 , 𝑖 𝜀 𝜕 2 𝑢 𝑖 𝜕 𝑦 2 + 𝜇 𝑟 , 𝑖 𝑘 1 𝑢 𝑖 = 0 . ( 4 . 4 ) This suggests that 𝑢 𝑖 ∼ 𝑂 ( 𝛿 ) 𝑝 𝑖 . In order to make this explicit and to make ordering of various quantities simpler, we represent the pressure 𝑝 𝑖 and the 𝑥 -component velocity 𝑢 𝑖 in the following manner: 𝑝 𝑖 = 𝑝 𝑖 ( 0 ) , 𝑢 𝑖 = 𝛿 𝑢 𝑖 ( 0 ) . ( 4 . 5 ) The above variables are the leading order quantities in an asymptotic expansion in 𝛿 , and we will be concerned only with the leading order variables in this paper. The continuity equation in both fluids ( 3.5 ) becomes, upon using 𝜕 / 𝜕 𝑥 ∼ 𝛿 𝜕 / 𝜕 𝜒 and using the above expansion for 𝑢 𝑖 , 𝛿 2 𝜕 𝑢 𝑖 ( 0 ) + 𝜕 𝜒 𝜕 𝜐 𝑖 𝜕 𝑦 = 0 , ( 4 . 6 ) which suggests the following expansion for 𝜐 𝑖 : 𝜐 𝑖 = 𝛿 2 𝜐 𝑖 ( 0 ) . ( 4 . 7 ) Upon using this expansion, the nondimensional 𝑦 component of the momentum equation yields to leading order in 𝛿 , 𝜕 𝑝 𝑖 ( 0 ) / 𝜕 𝑦 = 0 , implying that the pressure is constant in both films across the 𝑦 direction, and so 𝑝 𝑖 = 𝑝 𝑖 ( 𝜒 , 𝑡 ) . The simplified 𝑥 -momentum equation, therefore, is given by 𝑑 𝑝 𝑖 ( 0 ) − 𝜇 𝑑 𝜒 𝑟 , 𝑖 𝜀 𝜕 2 𝑢 𝑖 ( 0 ) 𝜕 𝑦 2 + 𝜇 𝑟 , 𝑖 𝑘 1 𝑢 𝑖 ( 0 ) = 0 . ( 4 . 8 ) The normal stress condition at the interface, ( 3.7 ), simplifies to give the following equation in the long-wave limit: 𝑝 1 ( 0 ) − 𝑝 2 ( 0 ) + 1 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝛾 𝛿 2 𝜕 2 ℎ 𝜕 𝜒 2 + 𝜖 1 2  𝜕 𝜓 1  𝜕 𝑦 2 − 𝜖 2 2  𝜕 𝜓 2  𝜕 𝑦 2 . ( 4 . 9 ) In order for the interfacial tension to be of the same order as the other terms in the above equation, we require 𝛾 𝛿 2 ∼ 𝑂 ( 1 ) , where 𝛿 = 𝐻 / 𝐿 . We set 𝛾 ( 𝐻 / 𝐿 ) 2 = 1 , and from this relation, we determine the lateral length scale 𝐿 to be √ 𝐿 = 𝛾 𝐻 2 =  𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 . Upon using the relation 𝛾 𝛿 2 = 1 , ( 4.9 ) becomes 𝑝 1 ( 0 ) − 𝑝 2 ( 0 ) + 1 𝑘 1  𝜑 2 − 𝜇 𝑟 𝜑 1  = 𝜕 2 ℎ 𝜕 𝜒 2 + 𝜖 1 2  𝜕 𝜓 1  𝜕 𝑦 2 − 𝜖 2 2  𝜕 𝜓 2  𝜕 𝑦 2 . ( 4 . 1 0 ) The tangential stress continuity, ( 3.8 ), simplifies to the following condition in the long-wave limit: 𝜕 𝑢 2 ( 0 ) 𝜕 𝑦 − 𝜇 𝑟 𝜕 𝑢 1 ( 0 ) 𝜕 𝑦 = − 𝑞 𝜕 𝜓 1 𝜕 𝜒 . ( 4 . 1 1 ) The nondimensional kinematic condition at the interface, ( 3.10 ), after using the asymptotic expansion for 𝜐 𝑖 ( 4.7 ), yields 𝛿 2 𝜐 1 ( 0 ) ( 𝑦 = ℎ ( 𝜒 ) ) = 𝜀 𝜕 ℎ 𝜕 𝑡 + 𝛿 2 𝑢 1 ( 0 ) 𝜕 ℎ 𝜕 𝜒 . ( 4 . 1 2 ) In order for the time derivative term in the above equation to be of the same order as the other two terms, it is necessary to stipulate a slow time scale in the long-wave limit such that 𝜕 𝜕 𝑡 = 𝛿 2 𝜕 𝜕 𝜏 , ( 4 . 1 3 ) where 𝜕 / 𝜕 𝜏 is 𝑂 ( 1 ) . The kinematic condition thus gives 𝜐 1 ( 0 ) ( 𝑦 = ℎ ( 𝜒 ) ) = 𝜀 𝜕 ℎ 𝜕 𝜏 + 𝑢 1 ( 0 ) 𝜕 ℎ 𝜕 𝜒 . ( 4 . 1 4 ) The dimensional slow time scale is obtained as follows: 𝜕 = 𝜇 𝜕 𝑡 2 𝐻 2 𝜖 0 𝜓 2 𝑏 𝜕 𝜕 𝑡 d i m = 𝛿 2 𝜕 , 𝜕 𝜏 ( 4 . 1 5 ) where 𝑡 d i m is the dimensional time. Thus, in order to nondimensionalize the dimensional time 𝑡 d i m in long-wave limit, the appropriate time scale is 𝜇 2 𝐻 2 / ( 𝜖 0 𝜓 2 𝑏 𝛿 2 ) . After using 𝛾 𝛿 2 = 1 to eliminate 𝛿 2 , the time scale becomes 𝜇 2 𝐻 3 𝛾 / ( 𝜖 0 𝜓 2 𝑏 ) 2 . Finally, the nondimensional interfacial charge balance, ( 3.11 ) is simplified in the long-wave limit as 𝜀 𝛿 2 𝜕 𝑞 𝜕 𝜏 + 𝛿 2 𝑢 1 ( 0 ) 𝜕 𝑞 𝜕 𝜒 − 𝛿 2 𝑞 𝜕 𝜐 1 ( 0 )  𝑆 𝜕 𝑦 = 𝜀 1 𝜕 𝜓 1 𝜕 𝑦 − 𝑆 2 𝜕 𝜓 2  . 𝜕 𝑦 ( 4 . 1 6 ) The above equation suggests that the nondimensional conductivities 𝑆 1 and 𝑆 2 both should scale as 𝛿 2 in order to balance the left side of the equation. So, we let 𝑆 𝑖 = 𝛿 2 𝑆 𝑖 ( 0 ) , 𝑖 = 1 , 2 , where 𝑆 𝑖 ( 0 ) = 𝜎 𝑖 𝜇 2 𝛾 𝐻 3 / ( 𝜖 3 0 𝜓 4 𝑏 ) . Upon using these rescaled conductivities, the charge conservation equation becomes, after using the continuity equation, ( 4.6 ), 𝜀 𝜕 𝑞 + 𝜕  𝑢 𝜕 𝜏 1 ( 0 ) 𝑞   𝑆 𝜕 𝜒 = 𝜀 1 ( 0 ) 𝜕 𝜓 1 𝜕 𝑦 − 𝑆 2 ( 0 ) 𝜕 𝜓 2  . 𝜕 𝑦 ( 4 . 1 7 ) This completes the derivation of the simplified governing equations in the long-wave limit. 5. Nonlinear Evolution Equations We now outline the derivation of the nonlinear evolution equations for the interfacial position ℎ ( 𝜒 , 𝜏 ) and surface charge density 𝑞 ( 𝜒 , 𝜏 ) . The simplified Laplacian for the potential 𝜓 𝑖 , ( 4.2 ), is easily solved along with the boundary conditions, ( 4.3 ), to give the following expressions for the potentials 𝜓 𝑖 ( 𝑖 = 1 , 2 ): 𝜓 1  𝜖 ( 𝜒 , 𝑦 ) = ( 𝛽 − 𝑦 ) 2 +  ( 1 + ℎ ( 𝜒 ) ) 𝑞 ( 𝜒 ) 𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  , 𝜓 ℎ ( 𝜒 ) 2 ( 𝜒 , 𝑦 ) = 𝛽 𝜖 2 − 𝜖 1  𝜖 𝑦 + 𝛽 ( 1 + 𝑦 ) 𝑞 ( 𝜒 ) + ℎ ( 𝜒 ) 1 − 𝜖 2 −  ( 1 + 𝑦 ) 𝑞 ( 𝜒 ) 𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  , ℎ ( 𝜒 ) ( 5 . 1 ) where the interfacial charge density 𝑞 ( 𝜒 , 𝜏 ) is determined below by the interface charge conservation equation ( 4.17 ). The simplified 𝑥 -momentum equation ( 4.8 ) can be integrated with respect to 𝑦 , since 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 is independent of 𝑦 and the two constants of integration that arise are determined by the boundary conditions 𝑢 1 ( 0 ) ( 𝛽 ) = 0 , 𝑢 1 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝑢 ( 0 ) i n t 𝑢 ( 𝜒 ) , 2 ( 0 ) ( − 1 ) = 0 , 𝑢 2 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝑢 ( 0 ) i n t ( 𝜒 ) . ( 5 . 2 ) Here, 𝑢 ( 0 ) i n t ( 𝜒 ) is the 𝑥 component of the velocity at the interface 𝑦 = ℎ ( 𝜒 ) , and this quantity will eventually be determined by using the tangential stress continuity condition ( 4.11 ). For the purposes of keeping the algebra tractable, it is found convenient to keep 𝑢 ( 0 ) i n t ( 𝜒 ) undetermined at present. The solutions for the 𝑥 -component velocities 𝑢 𝑖 ( 0 ) ( 𝑖 = 1 , 2 ) thus obtained are 𝑢 1 ( 0 ) 𝑘 ( 𝜒 , 𝑦 ) = 1 𝜇 𝑟 𝑑 𝑝 1 ( 0 ) ⎡ ⎢ ⎢ ⎣  √ 𝑑 𝜒 s i n h 𝜀 / 𝑘 1   √ ( 𝑦 − 𝛽 ) − s i n h 𝜀 / 𝑘 1  ( 𝑦 − ℎ ( 𝜒 ) )  √ s i n h 𝜀 / 𝑘 1 [ ]  ⎤ ⎥ ⎥ ⎦ ℎ ( 𝜒 ) − 𝛽 − 1 + 𝑢 ( 0 ) i n t  √ ( 𝜒 ) s i n h 𝜀 / 𝑘 1  ( 𝑦 − 𝛽 )  √ s i n h 𝜀 / 𝑘 1 [ ]  , 𝑢 ℎ ( 𝜒 ) − 𝛽 ( 5 . 3 ) 2 ( 0 ) ( 𝜒 , 𝑦 ) = 𝑘 1 𝑑 𝑝 2 ( 0 ) ⎡ ⎢ ⎢ ⎣  √ 𝑑 𝜒 s i n h 𝜀 / 𝑘 1   √ ( 𝑦 + 1 ) − s i n h 𝜀 / 𝑘 1  ( 𝑦 − ℎ ( 𝜒 ) )  √ s i n h 𝜀 / 𝑘 1  ⎤ ⎥ ⎥ ⎦ ( ℎ ( 𝜒 ) + 1 ) − 1 + 𝑢 ( 0 ) i n t  √ ( 𝜒 ) s i n h 𝜀 / 𝑘 1  ( 𝑦 + 1 )  √ s i n h 𝜀 / 𝑘 1  , ( ℎ ( 𝜒 ) + 1 ) ( 5 . 4 ) where the pressure gradients 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 ( 𝑖 = 1 , 2 ) are determined below. The continuity equation ( 4.6 ), after substituting the asymptotic expansion for 𝜐 𝑖 , ( 4.7 ), simplifies to 𝜕 𝑢 𝑖 ( 0 ) + 𝜕 𝜒 𝜕 𝜐 𝑖 ( 0 ) 𝜕 𝑦 = 0 . ( 5 . 5 ) The above equation is integrated with respect to 𝑦 from 𝑦 = 𝛽 to 𝑦 = ℎ ( 𝜒 ) for fluid 1 and 𝑦 = − 1 to 𝑦 = ℎ ( 𝜒 ) for fluid 2, to yield the following expressions for the normal velocities at the interface 𝜐 𝑖 ( 0 ) [ 𝑦 = ℎ ( 𝜒 ) ] for 𝑖 = 1 , 2 , where the integration constant is set to zero in order to satisfy the no-penetration condition at 𝑦 = 𝛽 and 𝑦 = − 1 : 𝜐 1 ( 0 ) [ ]  𝑦 = ℎ ( 𝜒 ) = − 𝛽 ℎ ( 𝜒 ) 𝜕 𝑢 1 ( 0 ) 𝜐 𝜕 𝜒 𝑑 𝑦 2 ( 0 ) [ ]  𝑦 = ℎ ( 𝜒 ) = − ℎ ( 𝜒 ) − 1 𝜕 𝑢 2 ( 0 ) 𝜕 𝜒 𝑑 𝑦 . ( 5 . 6 ) Note that the normal velocities of the two fluids are equal at the interface (normal velocity continuity condition), and so, 𝜐 1 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝜐 2 ( 0 ) [ ] 𝑦 = ℎ ( 𝜒 ) = 𝜐 ( 0 ) i n t . ( 5 . 7 ) Therefore, we equate the two integrals in the above equation and apply Leibnitz rule 𝜕  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 − 𝜕 ℎ 𝑢 𝜕 𝜒 1 ( 0 ) ℎ 𝜕 ( 𝜒 ) =  𝜕 𝜒 ℎ ( 𝜒 ) − 1 𝑢 2 ( 0 ) 𝑑 𝑦 − 𝜕 ℎ 𝑢 𝜕 𝜒 2 ( 0 ) ℎ ( 𝜒 ) . ( 5 . 8 ) Noting that the 𝑥 -component velocities are equal at the interface, 𝑢 1 ( 0 ) [ ℎ ( 𝜒 ) ] = 𝑢 2 ( 0 ) [ ℎ ( 𝜒 ) ] , the above equation is simplified to 𝜕   𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 )  𝑑 𝑦 − ℎ ( 𝜒 ) − 1 𝑢 2 ( 0 )  𝑑 𝑦 = 0 . ( 5 . 9 ) We next integrate the above equation with respect to 𝜒 , and set the integration constant (which is at most a function of time) to zero. The constant of integration is zero, because the pressure gradients 𝑑 𝑝 𝑖 ( 0 ) / 𝑑 𝜒 in the two fluids should be zero in the absence of electric field, and when the interfaces are flat. We then substitute the expressions for 𝑢 𝑖 ( 0 ) ( 𝑦 ) , ( 5.3 ) and ( 5.4 ), in the above equation, carry out the integrations with respect to 𝑦 , and substitute the simplified normal stress continuity condition ( 4.10 ) to eliminate 𝑑 𝑝 2 ( 0 ) / 𝑑 𝜒 in terms of 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 . Prior to determining 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 , it is useful to determine the 𝑥 -component of the fluid velocity at the interface 𝑢 ( 0 ) i n t from the simplified tangential stress continuity condition, ( 4.11 ). Once 𝑢 ( 0 ) i n t is determined, the pressure gradient 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 is determined from the integrated version of ( 5.9 ), and thus the velocity profile (( 5.3 ) and ( 5.4 )), is known completely. Therefore, we get 𝑑 𝑝 1 ( 0 ) =  𝑘 𝑑 𝜒 1  ℎ ( 1 + ℎ ( 𝑥 ) )      ( 𝜒 ) − 𝑀 ( 𝜒 ) − 𝐺 𝐹 𝑘 1 𝜀 − 𝑘 1  ℎ    (    𝜒 ) − 𝑀 ( 𝜒 ) 𝐺 + 4 𝑘 1 𝜀    t a n h 𝜀 4 𝑘 1 [ ]    𝑘 1 + ℎ ( 𝜒 ) 1  1 1 + ℎ ( 𝜒 ) + 𝜇 𝑟 [ ]    𝛽 − ℎ ( 𝜒 ) − 𝐺 t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) − t a n h 𝜀 4 𝑘 1 [ ] −  ℎ ( 𝜒 ) − 𝛽   4 𝑘 1 𝜀    t a n h 𝜀 4 𝑘 1 [ ]  − 1 1 + ℎ ( 𝜒 ) 𝜇 𝑟   t a n h 𝜀 4 𝑘 1 [ ] ℎ ( 𝜒 ) − 𝛽     − 1 , ( 5 . 1 0 ) 𝑑 𝑝 2 ( 0 ) = 𝑑 𝜒 𝑑 𝑝 1 ( 0 ) ( 𝜒 ) − 𝜕 𝑑 𝜒 3 ℎ ( 𝜒 ) 𝜕 𝜒 3 + 𝑢 ( 0 ) i n t ( 𝜒 ) 𝑘 1  𝜇 𝑟  𝑢 − 1 + 𝑀 ( 𝜒 ) , ( 5 . 1 1 ) ( 0 ) i n t   ( 𝜒 ) = 𝑘 1 𝜀 𝐹 − 𝑘 1 𝑑 𝑝 1 ( 0 )    𝑑 𝜒 t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) − t a n h 𝜀 4 𝑘 1 [ ] ℎ ( 𝜒 ) − 𝛽   + 𝑘 1   t a n h 𝜀 4 𝑘 1 [ ] 𝜕 1 + ℎ ( 𝜒 )   3 ℎ ( 𝜒 ) 𝜕 𝜒 3 ×    − 𝑀   c o t h 𝜀 𝑘 1 [ ]  1 + ℎ ( 𝜒 ) − 𝜇 𝑟   c o t h 𝜀 𝑘 1 [ ]  +  𝜇 ℎ ( 𝜒 ) − 𝛽 𝑟    − 1 t a n h 𝜀 4 𝑘 1 [ ] 1 + ℎ ( 𝜒 )   − 1 , ( 5 . 1 2 ) where 𝐺 , 𝑀 , and 𝐹 are functions of 𝜒 , and they are defined as    𝐺 ( 𝜒 ) = c o t h 𝜀 𝑘 1 [ ]  1 + ℎ ( 𝜒 ) − 𝜇 𝑟   c o t h 𝜀 𝑘 1 [ ]  +  𝜇 ℎ ( 𝜒 ) − 𝛽 𝑟    − 1 t a n h 𝜀 4 𝑘 1 [ ] 1 + ℎ ( 𝜒 )   − 1 ×   𝜇 𝑟  [ ] +  − 1 1 + ℎ ( 𝜒 ) 𝑘 1 𝜀   1 − 2 𝜇 𝑟    t a n h 𝜀 4 𝑘 1 [ ]    1 + ℎ ( 𝜒 ) + t a n h 𝜀 4 𝑘 1 [ ]  𝜖 ℎ ( 𝜒 ) − 𝛽    𝑀 ( 𝜒 ) = 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) − 3  𝜖 1 𝜖 2 ℎ   ( 𝜒 ) 𝑞 ( 𝜒 ) ( 1 + 𝛽 ) + 𝜖 2 − 𝜖 1 ×  𝜖   2 + 𝜖 1 [ ]  +  𝜖 + 𝑞 ( 𝜒 ) 1 + 2 ℎ ( 𝜒 ) − 𝛽 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) − 2 ×  𝑞 ( 𝜒 ) 𝑞   𝜖 ( 𝜒 ) 1 [ ] 1 + ℎ ( 𝜒 ) 2 + 𝜖 2 [ ] 𝛽 − ℎ ( 𝜒 ) 2  + 𝜖 1 𝜖 2 𝑞  [ ]  , 𝐹 ( 𝜒 ) 1 + 2 ℎ ( 𝜒 ) − 𝛽 ( 𝜒 ) = − 𝑞 ( 𝜒 ) ℎ  ( 𝜒 ) 𝜖 2  ( 𝛽 − ℎ ( 𝜒 ) ) 𝑞 ( 𝜒 ) ( 𝛽 + 1 ) − 𝜖 1 + 𝜖 2   𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   ℎ ( 𝜒 ) 2 − 𝑞 ( 𝜒 ) 𝑞  ( 𝜒 ) ( 𝛽 − ℎ ( 𝜒 ) ) ( 1 + ℎ ( 𝜒 ) )  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2   . ℎ ( 𝜒 ) ( 5 . 1 3 ) Finally, the kinematic condition at the interface, ( 4.14 ), is used to derive the evolution equation for ℎ ( 𝜒 , 𝜏 ) , as follows. We first substitute the expressions for the normal fluid velocities at the interface, ( 5.6 ), after using Leibnitz rule on the integral 𝜀 𝜕 ℎ 𝜕 𝜏 + 𝑢 1 ( 0 ) [ ] ℎ ( 𝜒 ) 𝜕 ℎ 𝜕 𝜕 𝜒 = −  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 + 𝑢 1 ( 0 ) [ ] ℎ ( 𝜒 ) 𝜕 ℎ , 𝜕 𝜒 ( 5 . 1 4 ) which is simplified to give 𝜀 𝜕 ℎ + 𝜕 𝜕 𝜏  𝜕 𝜒 𝛽 ℎ ( 𝜒 ) 𝑢 1 ( 0 ) 𝑑 𝑦 = 0 , ( 5 . 1 5 ) where 𝑢 1 ( 0 ) is substituted from ( 5.3 ), after using the expressions for 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 and 𝑢 ( 0 ) i n t determined using the procedure outlined above. The evolution equation for the interfacial charge density 𝑞 ( 𝜒 , 𝜏 ) is obtained from ( 4.17 ) after substituting the expressions for 𝑢 ( 0 ) i n t and the gradients of the potential (from ( 5.1 )) in that equation. The nonlinear evolution equations for ℎ ( 𝜒 , 𝜏 ) and 𝑞 ( 𝜒 , 𝜏 ) , respectively, take the forms 𝜀 𝜕 ℎ + 1 𝜕 𝜏 𝜇 𝑟 𝑘   1 𝑑 𝑝 1 ( 0 ) + 𝜇 𝑑 𝜒 𝑟 2 𝑢 ( 0 ) i n t  ( 𝜒 ) 𝜕 ℎ 𝜕 𝜒 s e c h 2  1 2  𝜀 𝑘 1  ( 𝛽 − ℎ ) + 𝑘 1  𝑑 ( 𝛽 − ℎ ) 2 𝑝 1 ( 0 ) 𝑑 𝜒 2 − 𝑑 𝑝 1 ( 0 ) 𝑑 𝜒 𝜕 ℎ  −  𝜕 𝜒 𝑘 1 𝜀  2 𝑘 1 𝑑 2 𝑝 1 ( 0 ) 𝑑 𝜒 2 + 𝜇 𝑟 𝑢  ( 0 ) i n t   1 ( 𝜒 ) × t a n h 2  𝜀 𝑘 1   𝜀 ( 𝛽 − ℎ ) = 0 , 𝜕 𝑞 𝜕 𝜏 + 𝑢 ( 0 ) i n t ( 𝜒 ) 𝜕 𝑞 𝜕 𝜒 − 𝑞 𝜕 ℎ  √ 𝜕 𝜒 𝜀 𝑘 1 𝜇 𝑟 𝑑 𝑝 1 ( 0 )  1 𝑑 𝜒 t a n h 2  𝜀 𝑘 1  +  ( ℎ − 𝛽 ) 𝜀 𝑘 1 𝑢 ( 0 ) i n t   ( 𝜒 ) c o t h 𝜀 𝑘 1   ( ℎ − 𝛽 ) + 𝑞 𝑢  ( 0 ) i n t + ( 𝜒 ) 𝜀 𝑆 1 ( 0 )  ( 1 + ℎ ) 𝑞 + 𝜖 2  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  ℎ + 𝜀 𝑆 2 ( 0 )  ( 𝛽 − ℎ ) 𝑞 − 𝜖 1  𝜖 1 + 𝛽 𝜖 2 +  𝜖 1 − 𝜖 2  ℎ = 0 . ( 5 . 1 6 ) These coupled nonlinear equations can be solved numerically with appropriate initial conditions to determine the evolution of the interface in the presence of electric fields and porous medium. In the present work, however, we restrict ourselves to studying the linear stability properties of these equations, an issue we turn to next. 6. Stability Analysis and Discussion Before linearizing the coupled nonlinear equations, it is necessary to first determine the base state about which we perturb. The steady base state we consider is that of stationary fluids with a flat interface ℎ ( 𝜒 , 𝜏 ) = 0 , and with a constant interfacial charge density 𝑞 0 which is independent of 𝜒 and 𝜏 . This base state interfacial charge 𝑞 0 is determined from ( 4.17 ) with the left side set to zero, since 𝜕 / 𝜕 𝜏 = 0 and 𝑢 𝑖 ( 0 ) = 0 in the base state 𝑆 1 ( 0 )  𝜕 𝜓 1  𝜕 𝑦 𝑦 = 0 = 𝑆 2 ( 0 )  𝜕 𝜓 2  𝜕 𝑦 𝑦 = 0 . ( 6 . 1 ) The derivatives of the potentials in the above equation are calculated from ( 5.1 ), with ℎ ( 𝜒 ) set to zero. This yields the following expression for the steady interfacial charge density: 𝑞 0 = 𝜖 1 𝑆 2 ( 0 ) − 𝜖 2 𝑆 1 ( 0 ) 𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 ) . ( 6 . 2 ) The variables ℎ and 𝑞 are now perturbed about their base state values ℎ ( 𝜒 , 𝜏 ) = ℎ 1 [ ] , 𝑞 e x p 𝑖 𝑘 𝜒 + 𝜔 𝜏 ( 𝜒 , 𝜏 ) = 𝑞 0 + 𝑞 1 [ ] , e x p 𝑖 𝑘 𝜒 + 𝜔 𝜏 ( 6 . 3 ) where ℎ 1 and 𝑞 1 are the amplitudes of the perturbations which are independent of 𝜒 and 𝜏 , 𝑘 is the nondimensional wavenumber based on the lateral length scale  𝐿 = 𝛾 𝐻 3 / 𝜖 0 𝜓 2 𝑏 and 𝜔 is the nondimensional growth rate based on the time scale 𝜇 2 𝐻 3 𝛾 / ( 𝜖 0 𝜓 2 𝑏 ) 2 . We substitute the expressions for 𝑑 𝑝 1 ( 0 ) / 𝑑 𝜒 and 𝑢 ( 0 ) i n t from ( 5.10 ), ( 5.12 ), and ( 6.3 ) in the nonlinear evolution equations ( 5.16 ), and apply Taylor expansion on the hyperbolic functions about ℎ ( 𝜒 ) = 0 . We then linearize the resulting coupled nonlinear evolution equations with respect to ℎ 1 and 𝑞 1 , to obtain a set linear homogeneous equations for ℎ 1 and 𝑞 1 of the form  𝐻 1 ℎ 1 +  𝑄 1 𝑞 1  𝐻 = 0 , 2 ℎ 1 +  𝑄 2 𝑞 1 = 0 , ( 6 . 4 ) where  𝐻 1  𝑅 = − 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  𝛽 4 𝜖 4 2 + 𝜖 4 1  𝑅 1 4 + 𝑘 4 𝑅 1 3 𝑅 1 − 𝑘 2 𝑅 1 3 𝑅 1 𝜖 2  − 𝜖 2 1 𝜖 2 2  6  𝑅 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  𝛽 2 + 𝑘 2 𝑅 1 3 𝑅 1 𝜖 2  + 𝜖 3 1 𝜖 2  𝑅 1 5  𝑅 − 4 𝛽 1 4 + 𝑘 4 𝑅 1 3 𝑅 1  + 𝑘 2 𝑅 1 3 𝑅 1 𝛽 𝜖 2  − 𝜖 1 𝜖 3 2  𝑅 1 5 + 4 𝛽 3  𝑅 1 4 + 𝑘 4 𝑅 1 3 R 1  + 𝑘 2 𝑅 1 3 𝑅 1 𝛽 𝜖 2  + 𝑞 2 0 ( 1 + 𝛽 ) 𝜖 2  − 𝑅 6 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  × 𝑅   1 𝐶 ( 𝛽 − 1 ) − 1 𝐵 1 𝛽  𝑘 1 𝜀  𝜖 1 − 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2   − 𝑞 0 𝜖 2  𝑅 8 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  ×  𝛽  2 𝑅 1 − 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 1 +  2 𝑅 1 − 𝐶 1 𝐵 1  𝛽 ( 𝛽 − 1 ) 𝑘 1 𝜀  𝜖 1 𝜖 2 + 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 2   + 𝐵 1 𝑅 5 𝜔 𝜀 3 / 2 𝜇 𝑟 ,  𝑄 1 = ( 𝛽 − 1 ) 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2 𝑅   1 5 + 𝑘 2 𝑅 1 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2   − 𝑞 0  𝑅 7 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3  𝜖 1 + 𝛽 𝜖 2  2 × 𝑅   1 𝐶 + 𝛽 1 𝐵 1  𝑘 1 𝜀  𝜖 1 + 𝛽 2  𝑅 1 + 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 ,  𝐻   2 =  𝜖 1 + 𝛽 𝜖 2   𝜀  𝜀 𝑅 2  𝑅 1 8 − 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 4 1 𝜇 r + 𝜖 3 1  𝑅 1 6 + 𝑅 3 𝑅 4 + 𝐵 1 𝑅 2 𝑅 4 𝑅 3 4  + 𝜖 2 1 𝜖 2  𝑅 1 7 + 3 𝛽 𝑅 3 𝑅 4 + 3 𝐵 1 𝑅 2 𝑅 4 𝑅 3 5  + 𝛽 𝜖 2  𝑅 3 0 𝑅 2 2  k 1 𝜀 q 0 𝜇 r + 𝛽 𝜖 2 2  𝑅 1 9 + 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 R 2 ×  𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑅 1 8 + 𝐵 1 𝑅 4 𝑆 1 ( 0 )  𝜖 2  𝜇 𝑟   + 𝜖 1  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2  𝑅 2 0 + 3 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  3 𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  4 𝛽 𝑅 1 8 + 3 𝐵 1 𝑅 4 𝑆 1 ( 0 ) − 𝛽 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 2  𝜇 𝑟    − 𝑅 2 𝑞 0 𝜇 𝑟  𝑅 9  𝑅 2 1 + 𝑅 2 3 + 𝑅 2 5 + 𝑅 2 8 − 𝑅 2 9  + 𝑅 3 2 𝑅 3 6   ,  𝑄 2 = 𝑅 2  𝜖 1 + 𝛽 𝜖 2  𝜇 𝑟  𝐵 1 𝑅 2 𝜀 2  𝜖 1 + 𝛽 𝜖 2  3  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  𝜖 + 𝜔 1 + 𝛽 𝜖 2   − 𝑞 0  𝑅 9  𝑅 2 2 + 𝑅 2 4 + 𝑅 2 6 + 𝑅 2 7  − 𝑅 3 1 𝑅 2 √ 𝜀 𝑘 1  𝜖 1 + 𝛽 𝜖 2  + 𝑅 3 3 𝑅 3 6 ,   ( 6 . 5 ) in which 𝐵 1 , 𝐵 2 , 𝐶 1 , 𝐶 2 , and 𝑅 1 − 𝑅 4 9 are given in the appendix. The set of linear homogeneous equations ( 6.4 ) for ℎ 1 and 𝑞 1 is written in the matrix form 𝑀 ⋅ 𝐶 𝑇 = 0 , where the vector 𝐶 = [ ℎ 1 , 𝑞 1 ] and the determinant of this matrix 𝑀 is set to zero for nontrivial solutions in order to obtain the characteristic equation for 𝜔 . This characteristic equation, which gives the growth rate 𝜔 as a function of 𝑘 , 𝑘 1 , 𝛽 , 𝜇 𝑟 , 𝑆 1 ( 0 ) , 𝑆 2 ( 0 ) , 𝜀 , 𝜖 1 , and 𝜖 2 , is a quadratic equation for 𝜔 . The roots of the characteristic equation for 𝜔 can be written as 𝜔 1 𝑅 = − 4 8 𝑅 4 9 −  𝑅 2 4 8  𝐻 + 2 2  𝑄 1 𝑅 4 9 − 2 𝑅 4 1 𝑅 4 7 𝑅 4 9 𝑅 4 9 , 𝜔 ( 6 . 6 ) 2 𝑅 = − 4 8 𝑅 4 9 +  𝑅 2 4 8  𝐻 + 2 2  𝑄 1 𝑅 4 9 − 2 𝑅 4 1 𝑅 4 7 𝑅 4 9 𝑅 4 9 . ( 6 . 7 ) The roots 𝜔 1 and 𝜔 2 are always real, with one of them 𝜔 1 is always negative, and the other one 𝜔 2 can be positive or negative depending on the choice of the system parameters. We can take 𝜔 = 𝜔 2 as a function of 𝑘 only in ( 6.7 ) if the values of the other parameters are known. Now, to see the effects of various parameters on the stability of the considered system, we calculate the growth rate 𝜔 given by ( 6.7 ) as a function of the wavenumber 𝑘 for different values of all physical parameters included in the analysis. These calculations are presented in Figures 2 – 7 for the general case of leaky-leaky dielectric fluids, Figures 8 and 9 for the limiting case of perfect-leaky dielectric fluids in which 𝑆 1 ( 0 ) = 0 , and Figures 10 – 13 for another limiting case of perfect-perfect dielectric fluids in which 𝑆 1 ( 0 ) = 𝑆 2 ( 0 ) = 0 , where we have given the growth rate 𝜔 against the wavenumber 𝑘 for the porosity of porous medium 𝜀 , medium permeability 𝑘 1 , nondimensional conductivities 𝑆 1 ( 0 ) , 𝑆 2 ( 0 ) , dielectric constants 𝜖 1 , 𝜖 2 , the ratio of the thicknesses of top and bottom fluids 𝛽 , and the ratio of viscosities of the two fluids 𝜇 𝑟 = 𝜇 1 / 𝜇 2 , respectively. Figure 2: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 2 0 and (b) 0 ≤ 𝑘 ≤ 5 . Figure 3: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑘 1 when 𝜀 = 0 . 2 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑘 1 = 5 , 1 0 , 1 2 , respectively. Figure 4: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑆 1 ( 0 ) when 𝜀 = 0 . 4 , 𝑘 1 = 1 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑆 1 ( 0 ) = 1 0 2 , 𝑆 1 ( 0 ) = 1 0 3 , 1 0 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 1 0 0 and (b) 0 ≤ 𝑘 ≤ 1 0 . Figure 5: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑆 2 ( 0 ) when 𝜀 = 0 . 4 , 𝑘 1 = 1 0 , 𝑆 1 ( 0 ) = 1 0 3 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑆 2 ( 0 ) = 1 0 2 , 1 0 4 , 1 0 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 1 0 0 and (b) 0 ≤ 𝑘 ≤ 3 . Figure 6: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 1 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 1 = 3 , 5 , 7 , respectively. Figure 7: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜇 𝑟 when 𝜀 = 0 . 4 , 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 1 0 3 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , and 𝜇 𝑟 = 1 , 2 , 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 3 0 , and (b) 0 ≤ 𝑘 ≤ 5 . Figure 8: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 2 0 , and (b) 0 ≤ 𝑘 ≤ 5 . Figure 9: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜇 𝑟 when 𝜀 = 0 . 4 , 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 1 0 4 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , and 𝜇 𝑟 = 1 , 2 , 4 , respectively, for the wavenumber ranges (a) 0 ≤ 𝑘 ≤ 3 0 and (b) 0 ≤ 𝑘 ≤ 1 0 . Figure 10: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜀 when 𝑘 1 = 3 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜀 = 0 . 1 , 0 . 3 , 0 . 5 , respectively. Figure 11: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝑘 1 when 𝜀 = 0 . 2 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝑘 1 = 5 , 1 0 , 1 2 , respectively. Figure 12: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 1 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 2 = 4 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 1 = 3 , 5 , 7 , respectively. Figure 13: Variation of growth rate 𝜔 with wavenumber 𝑘 for various values of 𝜖 2 when 𝜀 = 0 . 4 , 𝑘 1 = 5 , 𝑆 1 ( 0 ) = 0 , 𝑆 2 ( 0 ) = 0 , 𝜖 1 = 3 , 𝛽 = 0 . 5 , 𝜇 𝑟 = 1 , and 𝜖 2 = 4 , 6 , 1 0 , respectively. 6.1. Leaky Dielectric-Leaky Dielectric Interface This is the general case in which 𝑆 1 ( 0 ) ≠ 0 and 𝑆 2 ( 0 ) ≠ 0 , which is discussed in Figures 2 – 7 , when the conductivities of the upper and lower fluids are present in the analysis. Figures 2(a) and 2(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the porosity of porous medium 𝜀 . It is clear from Figure 2(b) that for small values of the porosity (e.g., 𝜀 = 0 . 1 ), the growth rate 𝜔 increases by increasing the wavenumber 𝑘 till a maximum value of 𝜀 after which 𝜔 decreases by increasing 𝑘 ; that is, there exists a maximum mode of instability only for small values of the porosity 𝜀 , while for any other value of 𝜀 , we found that 𝜔 decreases by increasing 𝑘 , that is. the system is always stable in this case. It is clear also from Figure 2(a) , by increasing the porosity values and at any wavenumber value, that the porosity of porous medium has a stabilizing effect for the wavenumber range 0 ≤ 𝑘 ≤ 1 5 , and it has a destabilizing effect for higher wavenumber values 𝑘 > 1 5 ; that is, the porosity of porous medium has a dual role on the stability of the considered system (stabilizing and then destabilizing) depends on the wavenumber range (less than or hifher than a critical wavenumber value 𝑘 = 1 5 , resp.). Figure 3 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the medium permeability 𝑘 1 . We conclude from this figure that for any wavenumber value 𝑘 ≤ 2 0 , the growth rate 𝜔 increases by increasing the medium permeability 𝑘 1 values, while for wavenumber values 𝑘 > 2 0 , all the curves correspond to different values of medium permeability 𝑘 1 coincide. This means that the medium permeability 𝑘 1 has a destabilizing effect for the wavenumber range 0 < 𝑘 ≤ 2 0 and it has no effect on the stability of the considered system for higher wavenumber values 𝑘 > 2 0 . Figures 4(a) and 4(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the conductivity of upper fluid 𝑆 1 ( 0 ) . It is clear from this figure that the conductivity 𝑆 1 ( 0 ) has a stabilizing effect for small wavenumber values and a destabilizing effect for high wavenumber values, as the growth rates 𝜔 decrease and increase by increasing the conductivity 𝑆 1 ( 0 ) values, respectively. It is also seen that between these two wavenumber ranges, the conductivity of upper fluid 𝑆 1 ( 0 ) has a dual role on the stability of the considered system; that is, it has a destabilizing effect for 𝑆 1 ( 0 ) values greater than 1 0 2 , while it has a stabilizing effect for 𝑆 1 ( 0 ) values greater than 1 0 4 . Therefore, we conclude that the conductivity of upper fluid 𝑆 1 ( 0 ) has different effects on the stability of the system depending on the chosen wavenumber range. Figures 5(a) and 5(b) shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the conductivity of lower fluid 𝑆 2 ( 0 ) . It indicates that the conductivity 𝑆 1 ( 0 ) has a destabilizing effect for the wavenumber range 0 < 𝑘 < 1 0 , while it has a stabilizing effect for higher wavenumber values 𝑘 ≥ 1 0 , since the growth rate 𝜔 increases in the first wavenumber range and decreases in the second wavenumber range by increasing the conductivity of lower fluid 𝑆 2 ( 0 ) . We conclude also from Figure 5(b) that there exists a mode of maximum instability for small values of the conductivity 𝑆 2 ( 0 ) which disappears for high values of the conductivity of lower fluid 𝑆 2 ( 0 ) ≥ 1 0 4 . Figure 6 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the dielectric constant of upper fluid 𝜖 1 . It is seen from this figure that there exists a critical wavenumber value 𝑘 = 4 . 6 , before which the growth rates decrease by increasing the dielectric constant 𝜖 1 , and after which they increase by increasing 𝜖 1 values. Thus, the dielectric constant of upper fluid 𝜖 1 has a stabilizing as well as a destabilizing effects, for wavenumber ranges before and after this critical wavenumber value, respectively. Similarly, the effects of both the dielectric constant of lower fluid 𝜖 2 , and the ratio of thicknesses of upper and lower fluids 𝛽 on the stability of the considered system are found to have opposite effects to the effect of the dielectric constant 𝜖 1 given by Figure 6 , but the corresponding figures are not given here. In other words, we conclude that the dielectric constant of lower fluid 𝜖 2 has a destabilizing as well as a stabilizing effects for wavenumber ranges before and after the same critical wavenumber value 𝑘 = 4 . 6 , respectively, while the ratio of thicknesses of upper and lower fluids 𝛽 has also a destabilizing as well as a stabilizing effects for wavenumber ranges before and after a less critical wavenumber value 𝑘 = 2 . 5 , respectively. Figures 7(a) and 7(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the ratio of viscosity of upper and lower fluids 𝜇 𝑟 . It indicates that the viscosity ratio 𝜇 𝑟 has a stabilizing effect for small wavenumber values, and it has also a destabilizing effect for higher wavenumber values, since the growth rate 𝜔 decreases and increases by increasing the increase of 𝜇 𝑟 , respectively. It is clear also from Figure 7(b) that for 𝜇 𝑟 > 1 , that is, when the viscosity of upper fluid is larger than the viscosity of lower fluid, there exists a mode of maximum instability which disappear for 𝜇 𝑟 ≤ 1 . 6.2. Perfect Dielectric-Leaky Dielectric Interface This is the limiting case in which 𝑆 1 ( 0 ) = 0 and 𝑆 2 ( 0 ) ≠ 0 , which is discussed in Figures 8 and 9 , when the conductivity of the upper fluid is not included in the analysis. Figures 8(a) and 8(b) show the variation of the growth rate 𝜔 against the wavenumber 𝑘 for various values of the porosity of porous medium 𝜀 . In comparison with Figures 2(a) and 2(b) , we conclude that the porosity of porous medium 𝜀 behaves in the same manner as in the previous case of leaky-leaky dielectric fluids, except that the values of growth rates 𝜔 are higher than their values in the previous case and the corresponding curves intersect the 𝑘 -axis at larger wavenumber values than in the previous case. We conclude also that for small wavenumber values, there exists a mode of maximum instability for porosity values 𝜀 ≤ 0 . 3 . Similarly, the effects of medium permeability 𝑘 1 , the conductivity of lower fluid 𝑆 2 ( 0 ) , the dielectric constants 𝜖 1 , 𝜖 2 , and the ratio of thicknesses of upper and lower fluids 𝛽 on the stability of the considered system are found to behave in the same manner as their effects in the previous case of leaky-leaky dielectric fluids, but figures are excluded to save space, except that in this case: for the effect of 𝑘 1 , the growth rate values are higher than their values in the previous case; for the effect of 𝑆 2 ( 0 ) , the obtained curves intersect the 𝑘 -axis at larger wavenumber values than in the previous case; for the effect of 𝜖 1 , there exists a mode of maximum instability for all values of 𝜖 1 ; for the effect of 𝜖 2 , the obtained curves are exactly similar to those obtained in the previous case; finally, for the effect of 𝛽 , there is a mode of maximum instability and the obtained curves intersect the 𝑘 -axis at bigger wavenumber values than those obtained in the previous case. Figures 9(a) and 9(b) shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for various values of the viscosity ratio of upper and lower fluids 𝜇 r . In comparison with Figures 7(a) and 7(b) , we conclude that the viscosity ratio 𝜇 𝑟 behaves in the same manner as in the previous case of leaky-leaky dielectric fluids with the only difference that in this case, for all values of 𝜇 𝑟 ⪋ 1 , there exists a mode of maximum instability and not only for 𝜇 𝑟 > 1 shown in the previous case. 6.3. Perfect Dielectric-Perfect Dielectric Interface This is the limiting case in which 𝑆 1 ( 0 ) = 𝑆 2 ( 0 ) = 0 , which is discussed in Figures 10 – 13 , when the conductivities of the upper and lower fluids are absent in the analysis. Figure 10 shows the variation of the growth rate 𝜔 versus the wavenumber 𝑘 for different values of the porosity of porous medium 𝜀 . In view of the above discussion, we conclude from this figure that the porosity 𝜀 has a slightly stabilizing effect for small wavenumber values 𝑘 ≤ 2 . 5 and that it has no effect on the stability of the considered system afterwards for higher wavenumber values. Figure 11 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for various values of the medium permeability 𝑘 1 . It is clear from this figure that the permeability 𝑘 1 has a destabilizing effect, since the growth rate 𝜔 increases by increasing the medium permeability values at any fixed wavenumber value. Similarly, the effect of dielectric constant of the upper fluid 𝜖 1 on the stability of the system is illustrated in Figure 12 , and it shows that the dielectric constant 𝜖 1 has a stabilizing effect since the growth rate 𝜔 decreases by increasing the dielectric constant 𝜖 1 at any wavenumber value. Figure 13 shows the variation of growth rate 𝜔 with the wavenumber 𝑘 for different values of the dielectric constant of lower fluid 𝜖 2 , and from which we conclude thatthe dielectric constant 𝜖 1 , and after which they increase by increasing 𝜖 2 has a destabilizing effect in the wavenumber range 0 < 𝑘 ≤ 4 . 6 , and it has no effect on the stability of the considered system afterwards for higher wavenumber values 𝑘 > 4 . 6 . Similarly, the effects of both the ratio of thicknesses of upper and lower fluids 𝛽 , and the ratio of viscosities of upper and lower fluids 𝜇 𝑟 , respectively, on the stability of the system are found to has the same and the opposite effects as the effect of the dielectric constant 𝜖 2 shown in Figure 13 , but the corresponding figures are not given here to avoid any kind of repitation; that is, the parameters 𝛽 and 𝜇 𝑟 have destabilizing and stabilizing effects for small wavenumber values, respectively. 7. Concluding Remarks In conclusion, we have provided a general formulation for analyzing the effect of an externally applied electric field on the stability and dynamics of the interface between two leaky dielectric fluids of arbitrary viscosities and conductivities in porous medium. A systematic long-wave asymptotic analysis was used to derive coupled nonlinear evolution equations for the position of the interface and free charge density at the interface. Attention was restricted to linearized stability of the coupled nonlinear equations and the effect of a variety of system parameters on the stability of the considered system. Two limiting cases are also studied, that is, the case of perfect-leaky dielectric fluids and the case of two perfect dielectric fluids, and recovered the previous studies in absence of porous medium. The obtained results in these limiting cases and the general case of two leaky dielectric fluids can be summarized as follows: (I) For perfect-perfect dielectric fluids, we conclude for small wavenumbers that (i) the porosity of porous medium 𝜀 , the dielectric constant of upper fluid 𝜖 1 , and the ratio of viscosity of upper and lower fluids 𝜇 𝑟 have stabilizing effects, (ii) the medium permeability 𝑘 1 , the dielectric constant of lower fluid 𝜖 2 , and the ratio of thickness of upper and lower fluids 𝛽 have destabilizing effects, (iii) at any value of these physical parameters there are no modes of maximum instability, that is. the system is always stable, (iv) these physical parameters have no effect on the stability of the system for high wavenumber values. (II) For perfect-leaky dielectric fluids, we found that (i) the conductivity of lower fluid S 2 ( 0 ) has a destabilizing effect for small wavenumbers and a stabilizing effect for high wavenumbers, (ii) for small wavenumber values, the physical parameters 𝜀 , 𝜖 1 , and 𝜇 𝑟 have stabilizing effects, while the parameters, 𝑘 1 , 𝜖 2 , and 𝛽 have destabilizing effects, as in the case of perfect-perfect dielectrics, (iii) for high wavenumber values, new regions of stability or instability appear; in other words, the physical parameters 𝜀 , 𝜖 1 , and 𝜇 r have destabilizing effects for high wavenumbers, while the parameters 𝜖 2 , and 𝛽 have stabilizing effects, and k 1 has no effect on the stability of the system for high wavenumber values, (iv) there exists a mode of maximum instability for some of these physical parameters which do not appear in the previous case of perfect-perfect dielectrics. (III) For leaky-leaky dielectric fluids, we found that (i) the conductivity of upper fluid S 1 ( 0 ) has a destabilizing effect for small wavenumbers 0 < 𝑘 < 𝑘 1 and also a stabilizing effect for high wavenumbers 𝑘 > 𝑘 2 , while it has a dual role on the stability of the considered system between them in the wavenumber range 𝑘 1 < 𝑘 < 𝑘 2 . (ii) the effects of all other physical parameters on the stability of the considered system behave in the same manner as their effects in the case of perfect-leaky dielectrics, except that in the case of perfect-leaky dielectrics the stability or instability regions occur more faster than the corresponding case of leaky-leaky dielectrics, and the maximum instability holds for more values of the physical parameters included in the analysis. It should be mentioned that the problem investigated in this article can be generalized to study the linear electrohydrodynamic instabilities at the interface between two immiscible fluids, either perfect or leaky dielectrics, subjected to alternating electric fields and moving through a porous medium in the limit of the electrode spacing being large compared to the wavelength of the perturbation using the Floquet theory analysis [ 37 ], and this case is now in a current research. Appendix 𝐵 1   = c o t h 𝜀 𝑘 1  + 𝜇 𝑟  𝛽  c o t h 𝜀 𝑘 1  +  𝜇 𝑟   1 − 1 t a n h 2  𝜀 𝑘 1  , 𝐵 2 =  𝜀 𝑘 1  − c s c h 2   𝜀 𝑘 1  + 𝜇 𝑟 c s c h 2  𝛽  𝜀 𝑘 1  + 1 2  𝜇 𝑟  − 1 s e c h 2  1 2  𝜀 𝑘 1 , 𝐶   1 = 𝜇 𝑟  − 1 − 𝑘 1 𝜀   2 𝜇 𝑟   1 − 1 t a n h 2  𝜀 𝑘 1   1 + t a n h 2 𝛽  𝜀 𝑘 1 , 𝐶   2 = 𝜇 𝑟 1 − 1 + 2   1 + 1 − 2 𝜇 𝑟  s e c h 2  1 2  𝜀 𝑘 1  − t a n h 2  1 2 𝛽  𝜀 𝑘 1 , 𝑅   1 = 𝑘 1 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀   1 t a n h 2  𝜀 𝑘 1   , 𝑅 − 1 2 = 1 𝜇 𝑟   𝛽 − 2 𝑘 1 𝜀  𝛽 t a n h 2  𝜀 𝑘 1  − 𝜇 𝑟  𝑅 1 𝑘 1 + 𝐶 1 𝐵 1  𝛽 t a n h 2  𝜀 𝑘 1  , 𝑅   3 = 𝐵 1 𝑅 2 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  − 𝜖 1 𝑆 2 ( 0 ) + 𝜖 2 𝑆 1 ( 0 )  , 𝑅 4  = − − 𝐵 2 𝐶 1 + 𝐵 1 𝐶 2 𝐵 2 1    1 t a n h 2  𝜀 𝑘 1   1 + t a n h 2 𝛽  𝜀 𝑘 1 , 𝑅   5 = 𝑅 2 2  𝜖 1 + 𝛽 𝜖 2  4 , 𝑅 6 = 𝜀 𝜖 1  1 ( 1 − 𝛽 ) s i n h 2  𝜀 𝑘 1   𝜖 + 𝛽 1 + 𝛽 𝜖 2   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1  , 𝑅 7 =  𝜖 1 + 𝛽 𝜖 2   𝜀  𝜖 1 + 𝛽 2 𝜖 2   1 s i n h 2  𝜀 𝑘 1   𝜖 + 𝛽 1 + 𝛽 𝜖 2   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1 , 𝑅   8 = 2 𝜀 𝜖 1  𝜖 2 + 𝛽 𝜖 1   1 s i n h 2  𝜀 𝑘 1   + 𝛽 𝛽 𝜖 2 + 𝜖 1 𝜖   2 − 𝜖 1   𝜀 𝑘 1  1 c o s h 2  𝜀 𝑘 1  , 𝑅 9  1 = t a n h 2  𝜀 𝑘 1   𝛽 + t a n h 2  𝜀 𝑘 1  , 𝑅 1 0 = 𝑘 2 𝑘 1 𝜀 𝑅 2  𝜖 1 + 𝛽 𝜖 2   1 s e c h 2  𝜀 𝑘 1  , 𝑅 1 1 = 𝐵 1  𝛽 √ √ 𝜀 − 2 𝑘 1  𝛽 t a n h 2  𝜀 𝑘 1 , 𝑅   1 2 = 𝑅 2 𝜇 𝑟 √ 𝑘 1  𝛽 t a n h 2  𝜀 𝑘 1  , 𝑅 1 3 = 𝑅 9 𝑅 1 2 + 𝑅 1 1 𝑅 2 , 𝑅 1 4 = 𝑅 2 𝑅 1 2 𝑘 5 1  1 t a n h 2  𝜀 𝑘 1  , 𝑅 1 5 = 𝜀 𝑅 1 0 𝑅 1 2  1 s i n h 2  𝜀 𝑘 1  , 𝑅 1 6 = 𝑅 2  4 𝑅 3 + 𝑞 0 𝜇 𝑟 𝑅 2  𝜀 𝐵 2  𝑆 1 ( 0 ) − 𝑆 2 ( 0 )  + 4 𝛽 𝑘 1 𝑘 4 𝜖 2  1 t a n h 2  𝜀 𝑘 1 , 𝑅    1 7 = 𝑅 2  4 𝑅 3 ( 2 𝛽 − 1 ) + 3 𝛽 q 0 𝜇 𝑟 𝑅 2  𝜀 𝐵 1  𝑆 1 ( 0 ) − 𝑆 2 ( 0 )  + 2 𝛽 𝑘 1 𝑘 4 𝜖 2  1 t a n h 2  𝜀 𝑘 1 , 𝑅    1 8 = 𝛽 𝑞 0 𝑘 1 𝑘 4 𝑅 2 𝜀  1 t a n h 2  𝜀 𝑘 1  , 𝑅 1 9 = 𝑅 2  − 4 𝑅 3 + 𝛽 𝜀 𝑞 0 𝜇 𝑟 𝐵 1 𝑅 2  𝑆 1 ( 0 ) − 𝑆 2 ( 0 ) , 𝑅   2 0 = 𝑅 2  4 ( 𝛽 − 2 ) 𝑅 2 + 3 𝛽 𝜀 𝑞 0 𝜇 𝑟 𝐵 1 𝑅 3  𝑆 1 ( 0 ) − 𝑆 2 ( 0 ) , 𝑅   2 1 = − 𝑘 1 𝑘 2 (  𝜖 1 + 𝛽 ) 1 + 𝛽 𝜖 2  𝑞 2 0 𝜖 2  ( 𝛽 − 1 ) 𝜀 𝜖 1 +  𝜖 1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 2 = 𝑞 0 𝑘 1 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2  𝜀  𝜖 1 + 𝛽 2 𝜖 2  −  𝜖 1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 3 = 𝑞 0 𝜖 2  𝜖 1 + 𝛽 𝜖 2  𝑘 1 𝑘 2  2 𝜀 𝜖 1  𝛽 𝜖 1 + 𝜖 2  +  𝜖 1 − 𝜖 2 𝜖   1 + 𝛽 𝜖 2  𝛽 𝐶 1 𝐵 1  𝜀 𝑘 1  , 𝑅 2 4 = 𝜀 ( 1 − 𝛽 ) 𝑘 1 𝑘 2 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2  2 , 𝑅 2 5 = 𝜀 𝑘 1 𝑘 2  𝜖 1 + 𝛽 𝜖 2 𝜖   3 1  𝑘 2 − 𝜖 2  + 3 𝛽 𝑘 2 𝜖 2 1 𝜖 2 + 𝛽 3 𝑘 2 𝜖 3 2 + 𝜖 1 𝜖 2 2  3 𝛽 2 𝑘 2 + 𝜖 2 , 𝑅   2 6 = − 𝜀 𝑞 0 𝑘 1 𝑘 2  𝐶 1 𝐵 1  + 2 𝑘 1 𝜀  𝜖 2 1  𝜖 1 + 𝛽 ( 2 + 𝛽 ) 𝜖 2   1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 7 = 𝜀 𝜖 2 𝑘 1 𝑘 2  𝐶 1 𝐵 1  + 2 𝑘 1 𝜀  ×  ( 𝛽 − 1 ) 𝜖 1  𝜖 1 + 𝛽 𝜖 2  2 − 𝛽 2 𝑞 0 𝜖 2  𝛽 2 𝜖 2 + ( 1 + 2 𝛽 ) 𝜖 1    1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 8 = 𝜀 𝑘 1 𝑘 2 𝜖 1 𝜖 2 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀  ×  𝛽  𝛽 2  𝑞 − 1 2 0 𝜖 2 +  𝜖 1 + 𝛽 𝜖 2 𝜖   2 1 − 𝜖 2 2  − 2 𝑞 0 𝜖 2   1 + 𝛽 2  𝜖 1 + 𝛽 𝜖 2   − 4 𝛽 𝑞 0 𝜖 2 1  𝑘 1 𝜀   1 t a n h 2  𝜀 𝑘 1  , 𝑅 2 9 = 𝜀 𝑘 1 𝑘 2 𝐶   1 𝐵 1  + 2 𝑘 1 𝜀  ×  𝑘 2  𝜖 4 1 + 4 𝛽 3 𝜖 1 𝜖 3 2 + 𝛽 4 𝜖 4 2  + 𝜖 2 1 𝜖 2   1 − 𝛽 2  𝑞 2 0 + 2 𝛽 𝑘 2  2 𝜖 1 + 3 𝛽 𝜖 2  𝐶    + 2 𝛽 1 𝐵 1  𝑞 0 𝜖 3 1 𝜖 2   1 t a n h 2  𝜀 𝑘 1  , 𝑅 3 0 = 𝛽 𝑞 0 𝑘 2 𝜖 2  𝜖 1 + 𝛽 𝜖 2   ( 1 + 𝛽 ) 𝑞 0 − 𝜖 1 + 𝜖 2  , 𝑅 3 1 = 𝛽 𝑞 0 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2 , 𝑅 3 2 = 𝑘 2 𝜖 1 𝜖 2  𝜖 1 + 𝛽 𝜖 2 (   𝛽 − 1 ) 𝑞 0 − 𝜖 1 − 𝜖 2 (   𝛽 + 1 ) 𝑞 0 − 𝜖 1 + 𝜖 2  , 𝑅 3 3 = 𝑘 2  𝜖 1 + 𝛽 𝜖 2  2  ( 𝛽 − 1 ) 𝜖 1 𝜖 2 − 𝑞 0  𝜖 1 + 𝛽 2 𝜖 2 , 𝑅   3 4 = 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑆 1 ( 0 ) − 3 𝛽 𝑆 2 ( 0 )  𝜖 2  , 𝑅 3 5 = 𝛽 𝜀 𝜇 𝑟  𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑆 1 ( 0 ) − 𝛽 𝑆 2 ( 0 )  𝜖 2  , 𝑅 3 6 = 𝜀 𝑘 1 𝑅 2  1 t a n h 2  𝜀 𝑘 1  , 𝑅 3 7 = 𝑅 1 4 + 𝑘 4 𝑅 1 𝑅 1 3 , 𝑅 3 8 =  𝑅 1 𝐶 ( 𝛽 − 1 ) − 1 𝐵 1 𝛽  𝑘 1 𝜀  𝜖 1 − 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 , 𝑅 3 9  = 𝛽 2 𝑅 1 − 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜖 2 1 +  2 𝑅 1 − 𝐶 1 𝐵 1  𝛽 ( 𝛽 − 1 ) 𝑘 1 𝜀  𝜖 1 𝜖 2 + 𝐶 1 𝐵 1 𝛽 2  𝑘 1 𝜀 𝜖 2 2 , 𝑅 4 0 =  𝑅 1 𝐶 + 𝛽 1 𝐵 1  𝑘 1 𝜀  𝜖 1 + 𝛽 2  𝑅 1 + 𝐶 1 𝐵 1  𝑘 1 𝜀  𝜀 2 , 𝑅 4 1 = − 𝑅 3 7 𝛽 4 𝜖 4 2 + 𝜖 4 1  𝑅 3 7 − 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  − 𝜖 2 1 𝜖 2 2  6 𝑅 3 7 𝛽 2 + 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  + 𝜖 3 1 𝜖 2  𝑅 1 5 − 4 𝛽 𝑅 3 7 + 𝛽 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  − 𝜖 1 𝜖 3 2  𝑅 1 5 + 4 𝛽 3 𝑅 3 7 + 𝛽 𝑘 2 𝑅 1 𝑅 1 3 𝜖 2  + ( 1 + 𝛽 ) 𝑞 2 0 𝜖 2  − 𝑅 6 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3 𝑅 3 8  𝜖 1 + 𝛽 𝜖 2   − 𝑞 0 𝜖 2  𝑅 8 𝑅 1 0 𝑅 1 2 + 𝑘 2 𝑅 1 3 𝑅 3 9  𝜖 1 + 𝛽 𝜖 2 , 𝑅   4 2 = 𝛽 𝜖 2  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2 ×  𝑅 1 9 + 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  𝑅 1 8 + 𝐵 1 𝑅 4 𝑆 1 ( 0 )  𝜖 2  𝜇 𝑟   , 𝑅 4 3 = 𝜖 1  𝑅 3 0 𝑅 2 2  𝑘 1 𝜀 𝑞 0 𝜇 𝑟 + 𝛽 𝜖 2 2 ×  𝑅 2 0 + 3 𝛽 𝑅 3 𝑅 4 + 𝜀 𝛽 𝑅 2  3 𝐵 1 𝑅 4 𝑞 0  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  +  4 𝛽 𝑅 1 8 + 3 𝐵 1 𝑅 4 𝑆 1 ( 0 ) − 𝛽 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 2  𝜇 𝑟   , 𝑅 4 4 = 𝜀 𝑅 2  𝑅 1 8 − 𝐵 1 𝑅 4 𝑆 2 ( 0 )  𝜖 4 1 𝜇 𝑟 + 𝜖 3 1  𝑅 1 6 + 𝑅 3 𝑅 4 + 𝐵 1 𝑅 2 𝑅 4 𝑅 3 4  + 𝜖 2 1 𝜖 2  𝑅 1 7 + 3 𝛽 𝑅 3 𝑅 4 + 3 𝐵 1 𝑅 2 𝑅 4 𝑅 3 5  , 𝑅 4 5 = − 𝑅 2  𝜖 1 + 𝛽 𝜖 2  𝜇 𝑟 𝑞 0  𝑅 9  𝑅 2 2 + 𝑅 2 4 + 𝑅 2 6 + 𝑅 2 7  − 𝑅 3 1 𝑅 2 √ 𝜀 𝑘 1  𝜖 1 + 𝛽 𝜖 2  + 𝑅 3 3 𝑅 3 6  , 𝑅 4 6 = 𝐵 1 𝑅 2 2 𝜀 2 𝜇 𝑟  𝜖 1 + 𝛽 𝜖 2  4 , 𝑅 4 7 = 𝑅 4 6  𝑆 1 ( 0 ) + 𝛽 𝑆 2 ( 0 )  + 𝑅 4 5 , 𝑅 4 8 = 𝑅 4 1 𝑅 4 6  𝜖 1 + 𝛽 𝜖 2  + 𝐵 1 𝑅 5 𝑅 4 7 𝜇 𝑟 𝜀 3 / 2 , 𝑅 4 9 = 2 𝐵 1 𝑅 5 𝑅 4 6 𝜇 𝑟 𝜀 3 / 2  𝜖 1 + 𝛽 𝜖 2  . ( A . 1 ) <h4>References</h4> G. I. Taylor and A. D. McEwan, “ The stability of a horizontal fluid interface in a vertical electric field ,” The Journal of Fluid Mechanics , vol. 22, pp. 1–15, 1965. J. R. Melcher and C. V. Smith, “ Electrohydrodynamic charge relaxation and interfacial perpendicular-field instability ,” Physics of Fluids , vol. 12, no. 4, pp. 778–790, 1969. D. A. Saville, “ Electrohydrodynamics: the Taylor-Melcher leaky dielectric model ,” Annual Review of Fluid Mechanics , vol. 29, pp. 27–64, 1997. D. J. Griffiths, Introduction to Electrohydrodynamics , Pearson Education, Delhi, India, 3rd edition, 2006. A. A. Mohamed, E. F. Elshehawey, and M. F. El-Sayed, “ Electrohydrodynamic stability of two superposed viscous fluids ,” Journal of Colloid And Interface Science , vol. 169, no. 1, pp. 65–78, 1995. N. T. M. El-Dabe, “ Electrohydrodynamic stability of two superposed elasticoviscous liquids in plane Couette flow ,” Journal of Mathematical Physics , vol. 28, no. 11, pp. 2791–2800, 1987. D. D. Joseph and Y. Y. Renardy, Fundamentals of Two-Fluid Dynamics—Part I: Mathematical Theory and Application , vol. 3 of Interdisciplinary Applied Mathematics , Springer, New York, NY, USA, 1993. K. Abdella and H. Rasmussen, “ Electrohydrodynamic instability of two superposed fluids in normal electric fields ,” Journal of Computational and Applied Mathematics , vol. 78, no. 1, pp. 33–61, 1997. J. R. Melcher, Continuum Electromechanics , MIT Press, Cambridge, UK, 1981. E. Schaffer, T. Thurn-Albrecht, T. P. Russell, and U. Steiner, “ Electrohydrodynamic instabilities in polymer films ,” Europhysics Letters , vol. 53, no. 4, pp. 518–523, 2001. Z. Lin, T. Kerle, S. M. Baker, et al., “ Electric field induced instabilities at liquid/liquid interfaces ,” Journal of Chemical Physics , vol. 114, no. 5, pp. 2377–2381, 2001. M. D. Morariu, N. E. Voicu, E. Schaffer, Z. Lin, T. P. Russell, and U. Steiner, “ Hierarchical structure formation and pattern replication induced by an electric field ,” Nature Materials , vol. 2, no. 1, pp. 48–52, 2003. O. Ozen, N. Aubry, D. T. Papageorgiou, and P. G. Petropoulos, “ Electrohydrodynamic linear stability of a two immiscible fluids in channel flow ,” Electrochimica Acta , vol. 51, pp. 5316–5323, 2006. L. Wu and S. Y. Chou, “ Electrohydrodynamic instability of a thin film of viscoelastic polymer underneath a lithographically manufactured mask ,” Journal of Non-Newtonian Fluid Mechanics , vol. 125, no. 2-3, pp. 91–99, 2005. R. G. Larson, Constitutive Equations for Polymer Melts and Solutions , Butterworth, Stoneham, Mass, USA, 1988. R. Verma, A. Sharma, K. Kargupta, and J. Bhaumik, “ Electric field induced instability and pattern formation in thin liquid films ,” Langmuir , vol. 21, no. 8, pp. 3710–3721, 2005. D. T. Papageorgiou and P. G. Petropoulos, “ Generation of interfacial instabilities in charged electrified viscous liquid films ,” Journal of Engineering Mathematics , vol. 50, no. 2-3, pp. 223–240, 2004. V. Shankar and A. Sharma, “ Instability of the interface between thin fluid films subjected to electric fields ,” Journal of Colloid and Interface Science , vol. 274, pp. 294–308, 2004. R. V. Craster and O. K. Matar, “ Electrically induced pattern formation in thin leaky dielectric films ,” Physics of Fluids , vol. 17, no. 3, Article ID 032104, 2005. F. Li, O. Ozen, N. Aubry, D. T. Papageorgiou, and P. G. Petropoulos, “ Linear stability of a two-fluid interface for electrohydrodynamic mixing in a channel ,” Journal of Fluid Mechanics , vol. 583, pp. 347–377, 2007. G. Tomar, V. Shankar, A. Sharma, and G. Biswas, “ Electrohydrodynamic instability of a confined viscoelastic liquid film ,” Journal of Non-Newtonian Fluid Mechanics , vol. 143, no. 2-3, pp. 120–130, 2007. G. Supeene, C. R. Koch, and S. Bhattacharjee, “ Deformation of a droplet in an electric field: nonlinear transient response in perfect and leaky dielectric media ,” Journal of Colloid and Interface Science , vol. 318, no. 2, pp. 463–476, 2008. S. Herminghaus, “Dynamical instability of thin liquid films between conducting media,” Physical Review Letters , vol. 83, pp. 2539–2542, 1999. L. F. Pease III and W. B. Russel, “ Linear stability analysis of thin leaky dielectric films subjected to electric fields ,” Journal of Non-Newtonian Fluid Mechanics , vol. 102, no. 2, pp. 233–250, 2002. D. B. Ingham and I. Pop, Transport Phenomena in Porous Media , Pergamon Press, Oxford UK, 1998. K. Vafai, Handbook of Porous Media , Marcel Dekker, New York, NY, USA, 2000. J. A. Del Rio and S. Whitaker, “ Electrohydrodynamics in porous media ,” Transport in Porous Media , vol. 44, no. 2, pp. 385–405, 2001. I. Pop and D. B. Ingham, Convevtive Heat Transfer: Mathematical and Computational Modeling of Viscous Fluids and Porous Media , Pergamon Press, Oxford, UK, 2001. D. A. Nield and A. Bejan, Convection in Porous Media , Springer, New York, NY, USA, 3rd edition, 2006. M. F. El-Sayed, “ Electrohydrodynamic instability of two superposed viscous streaming fluids through porous media ,” Canadian Journal of Physics , vol. 75, no. 7, pp. 499–508, 1997. M. F. El-Sayed, “ Effect of normal electric fields on Kelvin-Helmholtz instability for porous media with Darcian and Forchheimer flows ,” Physica A , vol. 255, no. 1-2, pp. 1–14, 1998. M. F. El-Sayed, “ Electrohydrodynamic instability of dielectric fluid layer between two semi-infinite identical conducting fluids in porous medium ,” Physica A , vol. 367, pp. 25–41, 2006. M. F. El-Sayed, “ Instability of two streaming conducting and dielectric bounded fluids in porous medium under time-varying electric field ,” Archive of Applied Mechanics , vol. 79, no. 1, pp. 19–39, 2009. M. F. El-Sayed, “ Thermohydrodynamic instability of viscous rotating dielectric fluid layer in porous medium with vertical ac electric field ,” Special Topics & Reviews in Porous Media , vol. 1, pp. 15–29, 2010. M. F. El-Sayed, “ Electrohydrodynamic instability conditions for two superposed dielectric bounded fluids streaming with fine dust in a porous medium ,” Journal of Porous Media , vol. 13, pp. 457–473, 2010. M. F. El-Sayed, G. M. Moatimid, and T. M. N. Metwaly, “Nonlinear electrohydrodynamic stability of two superposed streaming finite dielectric fluids in porous medium with interfacial surface charges,” Transport in Porous Media , vol. 86, pp. 559–578, 2011. P. Gambhire and R. M. Thaokar, “Electrohydrodynamic instabilities at interfaces subjected to alternating electric field,” Physics of Fluids , vol. 22, Article ID 064103, 2010. //

Journal

ISRN Applied MathematicsHindawi Publishing Corporation

Published: Aug 24, 2011

There are no references for this article.