PHYSICAL REVIEW E 89, 032409 (2014)
Electrokinetic model for electric-field-induced interfacial instabilities Priya Gambhire and Rochish Thaokar* Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai 400 076, India (Received 18 October 2013; revised manuscript received 19 December 2013; published 26 March 2014) Technology based on electric-field-induced instabilities on thin polymer film surfaces has emerged as a promising candidate for soft lithography. Typically, the instability is modeled using the perfect dielectric (PD) or the leaky dielectric (LD) model. These assume the electric diffuse layer to be infinitesimally large or small, respectively. In the present work we conduct stability analysis assuming a PD-electrolyte solution interface. The concentration of ions and, hence, the diffuse layer thickness is in general assumed to be of the same order as the electrolyte film thickness. The PD-LD models are then realized as limiting cases of the ratio of the double layer thickness to the film thickness. DOI: 10.1103/PhysRevE.89.032409
PACS number(s): 68.03.−g, 47.57.jd, 47.65.−d
I. INTRODUCTION
An interface between two fluids destabilizes under electric fields in the presence of a contrast in the electrical properties of the two fluids such as dielectric constant and electrical conductivity. This phenomenon has been termed electrohydrodynamic instability (EHDI) [1] and finds applications in techniques such as electrospraying (drops under electric field) [2,3], electrospinning (cylindrical jets under electric field) [4], electropatterning (or soft lithography, planar fluid interfaces under electric fields) [5], and so on. Pioneering work in this field was conducted in the 1960s with contributions from G. I. Taylor and J. R. Melcher [6–9], who studied various aspects of EHD such as drops and planar fluid deformation under electric field. While dealing with instabilities in dielectric materials, the leaky dielectric (LD) model was proposed [1,5,10–13]. It overcame the fallacies in the instability predictions obtained assuming either perfectly dielectric (PD) or perfectly conducting (PC) fluids. According to the leaky dielectric theory, the fluids are assumed to have an infinitesimal amount of free charge confined to the interface while the bulk fluid is free of charge. Moreover, the thermal motion of ions is neglected. This model can explain deformation patterns which are beyond the scope of the perfect dielectric or the perfect conductor theory, increasing its utility in the field of electrohydrodynamics. In 2000, following the work by Schaffer et al. [5] wherein a submicron thin polymer film is observed to deform into ordered patterns under electric field, the leaky dielectric model is modified to suit these thin-film systems (negligible gravity). The leaky dielectric model combined with a linear stability analysis yields a wavelength of the instability which is in good agreement with the literature. However, it has been recently shown that in certain parametric regimes, the assumptions made in the leaky dielectric model may be invalid [14]. The wavelength of the instability depends on parameters such as conductivity, applied electric field, and frequency which alter the inherent time scales in the system. Nonlinear analysis is used to seamlessly obtain wavelengths at all the ratios of time scales under both dc and ac fields.
*
[email protected]
1539-3755/2014/89(3)/032409(9)
Although general in approach, the above studies neglect the thermal motion of the free charge present in the system. In a real system, free ions, albeit infinitesimal, will still exhibit diffusion leading to a diffuse layer of counterions around charged surfaces or interfaces. This layer is characterized by its thickness, called the Debye length, κ −1 . The Poisson-Nernst-Planck equation governs the dynamics of the charges within electrolytic materials. This equation leads to the Boltzmann distribution at equilibrium when charge dynamics is unimportant. For the case of the potential being very less than kB T , the Debye-Huckel approximation is appropriate and has been used quite extensively in the field of electrokinetics to study membrane or fluid interfaces [15]. The coupling of the curvature of a bilayer membrane to its surface charges where the charges can flip to the other side and where they can diffuse on the same side of an interface results in a curvature-surface charge density coupling that leads to instabilities [16]. The Poisson-Boltzmann equation with the Debye-Huckel approximation is used to obtain an expression for the electrostatic energy of the system and different modes of instabilities are analyzed. A linear stability analysis of the electrokinetic model [17] of a floating bilayer with surface charges on either side exhibits an instability. The presence of the bilayer membrane adds a dielectric layer between two conducting layers of fluid. A diffusion equation for the surface charges enables a coupling between the surface undulations of the membrane and the charges on the surface which beyond a threshold lead to an instability. In another set of studies, the effect of the Debye layers on the electrostatic properties of an artificial cell membrane have been investigated [18]. In the limit of vanishing thickness and dielectric constant and in the absence of inherent charges on the membrane, the system reduces to that of a fluid-fluid interface. The system differs from the present one in that ion channels are assumed and a discontinuity is imposed on the electric field at the interface. The authors study the effect of a perturbation on this system in order to obtain the corrections to the elastic modulii. They also characterize the induced-charge electro-osmotic flow around the perturbed membrane. In a follow-up study [19], the authors proposed an improved zero-thickness electrokinetic model under dc fields. Using the Poisson-Boltzmann equation with the Debye-Huckel approximation in conjunction with Robin-type of boundary
032409-1
©2014 American Physical Society
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
conditions, they obtain a dispersion relation for the perturbed membrane along with expressions for the threshold voltages for the instability in the case of a conductive and a nonconductive membrane. They present this model as a basis for further development to systems of various sizes and charge densities. The zero-thickness membrane model discussed above differs from the present work in a couple of aspects. The membrane is conducting while its distance from the electrodes is quite large compared to the Debye layer thickness. Thereby a quasineutral environment is assumed. The bending modulus is suppressed if the electrode spacing is decreased to the order of microns. The objective of the present work is to study the instability at a fluid-liquid interface as a function of the thickness of the liquid film. As the thickness changes, the assumption of electroneutrality in the bulk of the fluid (that is, the fluid away from the Debye layers) becomes invalid. Therefore, the main objective of this work is to study the scenarios of overlapped Debye layers. A detailed study of system with overlapped Debye layers is reported by Qu and Li [20]. They study a fluid confined between two flat plates such that the double layers formed on the plates overlap. They argue that the Boltzmann approximation cannot be used for a system with overlapped Debye layers as this violates the assumption of electroneutrality in the bulk liquid. Using suitable boundary conditions they develop governing equations for this case. They show that at low plate distances the predictions from the classical theory differ substantially from those predicted by their model. Although they study the system of overlapped Debye layers, the system consists of a single fluid under a field. A study of the limits where the electrokinetic model might match the predictions from the electrohydrodynamic model is absent. Electrokinetic models vis-`a-vis electrohydrodynamic models have been analyzed in the study of drop deformation. Hua et al. [21] use the electrokinetic model in the context of conducting drops suspended in a dielectric medium. They argue that these systems hitherto studied via a conductor model can be approached using an electrokinetic model. They motivate the use of this model in systems with high curvature such as small droplets. They show in this system the conductor model results as a limiting case (for thin double layers) to the electrokinetic model. Zholkovskij et al. [22] studied the deformation of drops under electric field taking into account the presence of the double layers. In this generalized model, termed the electrokinetic (EK) model, the thickness of the Debye layers occurs as a parameter which affects the deformation of the drops. They show that at very low Debye layer thicknesses, the drop deformation matches that predicted by the electrohydrodynamic theory assuming perfectly conducting materials while at very high Debye thicknesses it matches that of the perfect dielectric theory. In the present work, the system under study is that of two stratified planar fluids. Under the action of an applied electric field, the free charges in the system migrate and accumulate at the interface and on the electrode surfaces. Under the combined action of the electric field and thermal motion, these charges are assumed to form diffuse layers at these surfaces. The field and the instability in such a system is
studied using the “electrokinetic model.” The objective of the present work is to study the effect of confinement of the fluid films to thicknesses as those practically used in experiments such as by Schaffer et al. [5]. The confinement of the liquid films result in the overlapping of the Debye layers formed within the fluid. This in turn violates the electroneutrality condition. The concentration and the potential within the bulk solution are not considered zero but are derived consistently (refer Appendix). To our knowledge, a thorough study of how the EK model predicts various regimes resulting from the steadily increasing confinement for a simple fluid-fluid interface has not been carried out. Apart from carrying out the above, a systematic evidence of how the EK model matches the predictions made by the leaky dielectric model in the limiting cases gives universality to the model studied. II. THE MODEL A. System description
The system consists of two fluid layers held between planar electrodes (Fig. 1). The longitudinal direction is denoted by x while the voltage is applied along the transverse direction, i.e., the y direction. Variations in the third direction (z) are not considered. The upper fluid is referred to as 1 and the lower fluid as 2. The viscosity, density, conductivity, and dielectric constant of the upper fluid are denoted by μ1 , ρ1 , σ1 , and 1 , respectively, and those of the lower fluid by μ2 , ρ2 , σ2 , and 2 . The position of the interface is denoted by y ∗ = 0 while the electrodes are held at y ∗ = βh0 and y ∗ = −h0 . The superscript ∗ denotes dimensional quantities. B. Governing equations
The continuity and the momentum balance equations are used to model the hydrodynamics in the system. In the dimensional form, in the presence of an applied electric field they can be written as ∇ ∗ · V ∗i = 0,
(1)
ρ(∂t∗ V ∗i + V ∗i · ∇ ∗ V ∗i ) = −∇ ∗ pi∗ + μi ∇ ∗2 V ∗i − ρf∗i ∇ ∗ φi∗ , (2)
FIG. 1. A schematic of the interface between two dielectric fluids subjected to a normal electric field. The dashed line is a schematic representation of the potential distribution in the two fluids.
032409-2
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
where the subscript i = 1,2 refers to fluid 1 and fluid 2, respectively; pi is the pressure; Vi = ueˆ x + v eˆ y is the velocity vector; and ρ and μ are the density and the viscosity of the fluids, respectively. The final term on the right-hand side of the momentum balance equation gives the contribution of the applied electric field to fluid flow. ρfi and φi are the free charge density and the electric potential in the fluid, respectively. When electroneutrality is assumed to be valid in the bulk solution, this term becomes negligible and the electrohydrodynamic coupling comes only through interfacial conditions. To derive the equation governing the charges in the system, the species conservation equation can be written as ∂n∗± = −∇ ∗ · (∓em± n∗± ∇ ∗ φ ∗ − m± kB T ∇ ∗ n∗± ), ∂t ∗
(3)
where e is the elementary charge; n∗+ ,n∗− are the number concentrations of the free charges in the fluid; and m± is the mobility. Subtracting the positive and negative number concentrations, and using ρf∗ = e(n∗+ − n∗− ),σ = e(m∗+ n∗+ + m∗− n∗− ),
(4)
potential, at a reference position in the fluid y = b, and kB and T are the Boltzmann constant and temperature, respectively. The details of the derivation of Eq. (8) are given in the Supplementary Material [23]. It also shows that for viscous time scales, the dynamics of the charges can be ignored, and, hence, the Poisson Boltzmann equation for potential is decoupled. For the DH approximation to be valid, the term eφ ∗ /(kB T ) needs typically φ = 25 mV for it to be O(1). The voltages applied in EHD-based soft lithography are typically much higher, and ideally the Poisson-Boltzmann (PB) equations should be used. However, in the past, the DH model has been used as a substitute to the more complicated, nonlinear PB model to understand the physics analytically, which is also the endeavor in the present study. For further reading on the DH approximation and solving the nonlinear PB equation the authors cite the work by Lekkerkerker [24], Goldstein et al. [25], Vlahovska et al. [26]. C. Boundary conditions
The following set of boundary conditions are used to solve the governing equations. At the interface, the velocities are continuous in the normal and tangential directions,
results in the following equation: ∂ ∗ ρf∗ ∂t ∗
= −∇ ∗ · (σ ∇ ∗ φ ∗ − D∇ ∗ ρf∗ ).
(5)
process time scale τp =
μ2 h20 . 0 φ02 2
(6) 0 σ
and the
For a typical air-silicone oil
system, σ = 10−12 C2 /Nm s, 0 = 8.85 × 10−12 C2 /Nm2 , h0 = 100 × 10−6 m, μ2 = 30 Pa s, and φ0 = 30 V. Using these values τc ∼ 1 s and τp ∼ 37.6 s. For an air-water system, with the typical values of σ = 5 × 10−6 C2 /Nm2 s, 0 = 8.85 × 10−12 C2 /Nm2 , h0 = 100 × 10−6 m, μ2 = 10−3 Pa s, and φ0 = 30 V, the values of τc and τp are 1.7 × 10−6 s and 1.26 × 10−3 , respectively. Therefore, the left-hand sides of Eqs. (3) and (6) are equivalent to zero and the electrostatic equation therefore reduces to the Poisson equation [23], 0 i ∇ ∗2 φi∗ = −ρf∗i ,
(V1 ∗ · tˆ∗ ) = (V2 ∗ · tˆ∗ ),
(10)
1+∂x h
(8)
where κi∗ is the length. For fluid 1, κ1∗ = 0, while for √ Debye 2e2 no ∗ fluid 2, κ2 = 0 2 kB T , n0 , and φb are the uniform number concentration of the ions (in the absence of field) and the
1+∂x h
subscripts x, y, or t denotes differentiation with respect to that variable, h∗ (x,t) denotes the position of the interface, and V1 ∗ and V2 ∗ are the velocity vectors for the two fluids. For an unperturbed interface nˆ = eˆ y and tˆ = eˆ x . The stresses are balanced in the normal and tangential directions. In the normal direction the stresses satisfy [[nˆ ∗ · τi∗ · nˆ ∗ ]] = γ C ∗ ,
(11)
where the operator [[Xi ]] denotes the jump X1 –X2 across the interface h(x,t). C ∗ is the curvature given by C ∗ = ∇ ∗ · nˆ ∗ = −∂x∗2 h∗ /(1 + ∂x∗ h∗2 )3/2 . The term on the right-hand side in Eq. (11) represents the stabilizing force due to surface tension. The balance of stresses in the tangential direction is given by [[ tˆ∗ · τi ∗ · nˆ ∗ ]] = 0,
(7)
where 0 is the permittivity of vacuum and i is the dielectric constant of the two fluid layers i = 1,2, respectively. For the case of an air-liquid interface, ρfi , the free charge density in the top fluid is 0 while that in the lower liquid is given by ρf2 = e(n∗+ − n∗− ). This equation under the Debye-Huckel (DH) approximation reduces to ∇ ∗2 φi∗ = κi∗2 (φi∗ − φb∗ ),
(9)
where nˆ ∗ and tˆ∗ are the unit normal and the unit tangential vectors, respectively, as indicated in Fig. 1, given by −∂ ∗ h∗ eˆ +eˆ eˆ +∂ ∗ h∗ eˆ nˆ ∗ = √x x∗ ∗2y and tˆ∗ = √x x ∗ ∗2y . The operator ∂ with the
Upon scaling, the equation reduces to τc ∂ρf 0 D = −∇ · ∇φ − ∇ρ f , τp ∂t σ h20 where the charge relaxation time scale τc =
(V1 ∗ · nˆ ∗ ) = (V2 ∗ · nˆ ∗ ),
(12)
where τi ∗ , the total stress, i.e., the sum of the hydrodynamic and the electrical stresses, for the fluid i, is given by τi ∗ = −pi ∗ I + μi (∇ ∗ Vi ∗ + ∇ ∗ Vi ∗T ) + Mi ∗ .
(13)
The superscript T indicates transpose, i = 1,2 denotes the two fluid layers and Mi is the Maxwell stress tensor [27], Mi ∗ = i 0 E ∗i E ∗i − 12 (E ∗i · E ∗i )I . (14) A balance of the tangential component of the field gives the continuity of the potentials across the interface
032409-3
φ1∗ = φ2∗
(15)
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
and the balance of the normal component of electric field gives the electric displacement continuity equation [[0 i (−∇ ∗ φi∗ · nˆ ∗ )]] = 0.
(16)
The boundary condition assumes no adsorbed charges at the interface. It should be noted that although the bulk charge is nonzero the surface charge is zero. In fact, the net bulk charge integrated over the double layer (in the electrokinetic model) is really the surface charge integrated over the area (in the perfect dielectric-leaky dielectric model). This clearly differs from the free charge BC in PD-LD systems, where [[0 i (−∇ ∗ φi∗ · nˆ ∗ )]] = q aided by the equation ∂q = [[σi∗ ∇ ∗ φi∗ · nˆ ∗ ]]. The ∂t dynamics of the position of the interface is given by the kinematic condition ∂t h∗ + V1 ∗ · ∇ ∗ s h∗ = v2∗ = v1∗ .
steady state from the perturbed state. Any random function can be simply expressed as a sum of Fourier modes. As in the linear limit only O(1) terms need to be considered, LSA gives a neat, simplified mathematical formalism to study the stability of any system. In the following sections the methodology of LSA applied to the present system is described. A. Addition of perturbations
The system is subjected to an infinitesimal arbitrary perturbation. The system variables are therefore expressed as a sum of a mean value (denoted by an overbar on the variable) and a perturbation term (denoted by the superscript tilde). Thus, ¯ ˜ m(x,y,t) = m(x,y,t) + ηm,
(17)
D. Scalings
The parameters in the above equations are scaled using their characteristic values in the system. That is, length, in both the x and y directions, is scaled using the thickness of the thin film h0 . Potential is scaled by φ0 , the applied potential. Velocity is scaled using the term 0 φ02 /(μ2 h0 ) and pressure is scaled using 0 φ02 / h20 while the scaling for time is obtained to be μ2 h20 /(0 φ02 ).
where η is a small variable which gives the deviation from the base state. The perturbations are incorporated in the governing equations. The leading-order terms give the base state of the system while order η variables are considered in the linear stability analysis. 1. Base state
In the present case, in the base state the interface is flat with no perturbations and, hence, the mean velocities are zero,
E. Governing equations
The equations, after incorporating the dimensionless variables listed in Sec. II D are as follows: ∇ · V i = 0,
(18)
umi = 0, vmi = 0.
pm1 = 0,
(19)
−∇pi + μr ∇ 2 V i + i ∇ 2 φi ∇φi = 0.
(20)
pm2 =
∇ φi = κ (φi − φb ), 2
1 2 2
∂φ ∂y
(24)
2 + Pc .
Note that the hydrostatic equations for PD and LD cases is ∂y pi = 0 while for the EK model it is ∂y p = ρf E = κ 2 φ∇φ. The mean potential or the base-state potential is calculated by solving the following equation (see the Appendix for derivation) for fluids 1 and 2: ∂y2 φmi = κi2 (φmi − φb ).
(21)
(25)
The constant Pc is calculated from a stress balance such that ∂φm2 2 ∂φm1 2 1 pm2 = − 1 . (26) 2 2 ∂y ∂y
The nondimensional Poisson equation is 2
(23)
The base-state pressure is obtained from the (y-direction) momentum balance equation in both fluids,
Rei (∂t V i + V i · ∇V i ) = −∇pi + μr ∇ 2 V i + i ∇ 2 φi ∇φi , where μr = μ1 /μ2 (note: the nondimensional variables are denoted without the superscript ∗). The Reynold’s number (Re) is given by Rei = ρi 0 φ02 /μ22 for dc fields. For a typical system of air (fluid 1) and silicone oil (fluid 2) the properties are ρ1 ∼ 1 kg/m3 , ρ2 ∼ 970 kg/m3 , μ2 ∼ 0.97 Pa s, and h0 ∼ 100 × 10−9 m and, using φ0 ∼ 30 V, the Reynolds number for the dc case is of the order of 10−8 –10−5 . Hence, in the ensuing analysis, we consider the case where Rei are approximately zero. The momentum equations with the above parameter values are simplified to
(22)
(27)
where κ is the nondimensional inverse Debye layer thickness κ = κ ∗ h0 .
The following (scaled) boundary conditions, which include the potential at the electrodes and the continuity of potential and electric displacement at the interface, are used to solve the above equation:
III. LINEAR STABILITY ANALYSIS
φm1 |y=β = 0,
(28)
The linear stability analysis (LSA) is an important tool to study the response of a system to random disturbances. The basic principle of this theory is that any system always undergoes infinitesimal perturbations in its steady or equilibrium state. A system can be called linearly stable if it returns to its
φm2 |y=−1 = 1,
(29)
φm1 = φm2 |y=0 ,
(30)
[[−i ∂y φmi ]] = 0.
(31)
032409-4
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
The reference potential φb is determined using the condition 0 0 ρ2 dy = 2 κ22 (φm2 − φb )dy = 0. −1
−1
The base state values of the potential in the two fluids assuming κ1 = 0 and κ2 = κ are given by (1 + eκ )(y − β)2 κ , (32) 2(−1 + eκ )1 + (1 + eκ )β2 κ
β2 κcosh(κ/2) + 1 sinh(κ/2) − sinh 12 + y κ . = β2 κcosh(κ/2) + 21 sinh(κ/2) (33) φm1 = −
φm2
In Fig. 2, the potential distribution in the fluids in the base state is plotted along the normal direction y. The top fluid is bounded at y = 1 and y = 0 to −1 denotes the lower fluid. The top fluid is assumed to be air (a perfect dielectric) while the potential in the lower fluid is studied as a function of the concentration of ions through the parameter, inverse Debye length κ. The origin of the Debye layer in this case is due to the charge separation under electric field. Ions are attracted to respective electrodes while the counterions are repelled and accumulate at the interface. In the presence of thermal forces, these charges do not exist in a thin layer at the interface but form a diffuse layer with the concentration of ions decaying away from the interface or electrode. This leads to a Debye layer formation at both the interface and the electrode in the lower fluid layer giving rise to a potential distribution that is dependent upon the concentration of ions in the layer. As discussed earlier, increase in the concentration of ions leads to the formation of thinner, mroe compact Debye layers while lower concentration of ions leads to thicker, more diffuse Debye layers. In Fig. 2, the potential distribution for κ = 10−3 (solid line), κ = 5 (dashed line), κ = 10 (dotted line), and κ = 103 (dash-dot line) progressively representing thinner
bilayers, in the bottom fluid, is plotted along y. In the top fluid (y = 0 − 1), the potential distribution is always linear 32, a characteristic of the potential drop across a dielectric medium which is governed by the Laplacian of the potential being zero, while in the lower fluid (y = 0 to −1), the potential drop is given by Eq. (33) which shows the characteristic exponential drop within the Debye layers at the interface (y = 0) and the electrode surface (y = −1). The value κ = 1 denotes Debye layers that are as thick as the film thickness. In such a situation, the Debye layers below the interface and the layer above the electrode surface overlap. At low concentrations of ions, the overlap of potentials from the two Debye layers is significant to the extent that the potential drop becomes linear and the system behaves like a dielectric material until the concentration is so low that the potential drop in the system is equivalent to that in a PD material [refer Fig. 3]. As κ increases beyond 1, Debye layers become thinner, the overlap decreases, and the field at midplane of the fluid layer becomes zero. With increase in the concentration of ions (conductivity) or κ the Debye layers asymptote to zero thicknesses and, correspondingly, show perfect conductor-like behavior. The intricate physics of compact Debye layers which are referred to as Stern layers is beyond the scope of this paper and more details can be found in Refs. [19,28]. To get a perspective of the practicality of the κ values used in the study consider if the film thickness of 100 μm is a fixed parameter also typically encountered in experiments, the range of κ studied in the present work corresponds to the Debye layer thicknesses of 0.1 μm to 0.1 m. Although the value of 0.1 m appears unrealistic from an electrokinetics point-of-view the high thickness of the Debye layer means that the fluid has a very low concentration of charge spread over a large distance and the characteristics tend to that of a perfect dielectric material. A Debye layer thickness of 0.1 μm, on the other hand, is quite realistic. To add another perspective to the term κ. For a fixed fluid-fluid system, say, air-silicone oil, having thinner fluid film changes the wavelength of the EHD instability [5]. The silicone
1 0
0.5 y
−0.2
0
y
−0.4 −0.6
−0.5
−0.8
−1 0
0.2
0.4
φ
m
0.6
0.8
1 −1 0.75
FIG. 2. Potential distribution (φm ) in the two fluid layers [y = 1 to 0 indicates the top fluid (1) and y = 0 to −1 indicates the lower fluid (2)] in the base state at different values of κ. The solid line indicates κ = 10−3 , the dashed line indicates κ = 5, the dotted line indicates κ = 10 and the dash-dot line indicates κ = 1000. The potential distributions obtained from the PD-PD and the PD-PC theories coincide exactly with the κ = 10−3 and κ = 1000 curves, respectively, and, hence, cannot be distinguished. The values of parameters used are 1 = 1, 2 = 3, σ1 = 0 (implies κ1 = 0), μr = 0, and β = 1.
0.8
0.85
φ
m2
0.9
0.95
1
FIG. 3. Potential distribution (φm2 ) in the lower fluid layer, in the base state, at different values of κ. The dashed line indicates κ = 10−4 , the dotted line indicates κ = 1, the dash-dot line curve indicates κ = 5, the curve with ◦ legends indicates κ = 10, the solid line curve indicates κ = 102 , and the curve with the + legends indicates κ = 103 . The potential distributions obtained from the PD-PD and the PD-PC theories coincide exactly with the κ = 10−4 and κ = 1000 curves, respectively, and, hence, have not be plotted. The values of parameters used are 1 = 1, 2 = 3, μr = 0, and β = 1.
032409-5
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
IV. STABILITY ANALYSIS
The equations governing the perturbation variables [obtained by collecting the order η terms from Eqs. (18), (20), and (21)] after substituting the normal mode form for the variables are (34) ∂x u˜ i + ∂y v˜i = 0, 2
−∂x p˜i + Ci ∂x u˜ i + ∂y2 u˜i + i ∂y2 φmi ∂x φ˜ i = 0, (35) 2
2
2 2 −∂y p˜i + Ci ∂x v˜i + ∂y v˜i + i ∂x φ˜ i + ∂y φ˜ i ∂y φmi + ∂y2 φmi ∂y φ˜ i = 0, ∂x2 φ˜i + ∂y2 φ˜ i = κi2 φ˜ i .
(36) (37)
(a) 1
q
0.8 0.6 0.4 0.2 0
−3
10
(b)
−2
−1
10
κ
10
0
1
10
2
10
3
10
10
0.5
0.4
τ
yy
m
0.45
0.35 0.3 −3
10
(c) 0.5
−2
−1
10
1
κ
10
2
10
10
3
10
(a) (b) (c)
0.25 pm and τe
oil with an approximate conductivity of ∼10−12 S/m has a Debye layer of ∼10−6 m. The value κ = 103 corresponds to a 4-mm-thick film which will show a conductor-like behavior because the κ −1 h0 , while κ = 10−2 corresponds to an approximately O(10)-nm film. For the same fluids (silicone oil), this system would now behave as a perfect dielectric. Thus for a given fluid with a κ ∗ , the relative thickness of the film determines the extent of bulk charges. The gradual variation in the mean potential distribution φm2 with changing inverse Debye length κ in the lower fluid is shown in the Fig. 3. The infinitesimal exponential potential drops at the surfaces for κ = 102 and increases with decreasing κ until they merge (dash-dot line κ = 5 in Fig. 3) and become linear at κ 1. Initially behaving like two PD fluids, a PD-LD interface actually behaves like a PD-PC (under dc fields and at long times) interface as the potential drop in the lower fluid is balanced by the charge separation across it and the effective field becomes zero. The high κ limit of the EK theory exactly matches the results from the PD-PC theory. The free charge at the interface between the two fluids, given 0 by −1/2 2 κ 2 (φm2 − φb ), can be used to confirm the limiting behavior of the EK model. In Fig. 4(a), the free charge is plotted as a function of κ denoted by the solid line. At very low κ, the free charge is almost zero. But as κ increases, the free charge in the Debye layer increases until, at high-enough κ, it matches the accumulated interfacial charge calculated using a perfect conductor model [denoted by the dashed line in Fig. 4(a)]. In Fig. 4(b), the Maxwell’s stress in the top fluid in the base state is plotted as a function of κ. In the top fluid, the Maxwell’s stress shows a similar trend as the free charge at the interface. It asymptotes to the value given by the PD-PD theory at low κ and that of a PD-PC theory at high κ. But in the lower fluid, the presence of the Debye layer leads to a different scenario [see Fig. 4(c)]. The base-state pressure given by Eq. (26) balances the Maxwell’s stress such that the total stress [Eq. (13)] follows a similar trend as that of the Maxwell’s stress in the top fluid. That is, the total stress from the EK theory [curve (b) in Fig. 4(c)] asymptotes to the limiting cases of PD-PD and PD-PC theories at low and high κ, respectively.
(d) (e) (f)
0
(g)
−0.25 −0.5 −4 10
(h) (i) −3
10
−2
10
−1
10
0
κ 10
1
10
2
10
3
10
FIG. 4. (a) A comparison of charge (q) accumulated at the interface obtained using the PD-PC model (indicated by the dashed line) and the free charge present in the Debye layer near the interface, in the lower fluid (ρ) obtained using the electrokinetic model (indicated by the solid line) varying as a function of the inverse Debye length κ. (b) Variation of Maxwell’s stress in the top fluid (τyy1 ), at the interface (y = 0), in the base state, as a function of the inverse Debye length κ. The dotted line indicates results from the PD-PD theory, the dashed line indicates results from the PD-PC theory and the solid lined curve indicates results from the electrokinetic model. (c) Variation of mean pressure [pm2 , curves (d), (e), and (f)], the electric stress [τe curves (g), (h), and (i)], and the total stress [τyy2 , curves (a), (b), and (c)] in the lower fluid, in the base state, as a function of the inverse Debye length κ. The dotted lines indicate results from the PD-PD theory, the dashed lines indicate results from the PD-PC theory, and the solid lined curve indicates results from the electrokinetic model. The values of parameters used are 1 = 1, 2 = 3, μr = 0, and β = 1.
The perturbation variables are now expressed in the normal ˜ is the perturbation variable mode form. Thus, if m ˜ ˆ m(x,y,t) = m(y,t)exp(ikx)dk, (38) ˆ where m(y,t) is the amplitude of the perturbation and k is the wave number of the perturbation. Substituting this form of the perturbation term in the Eqs. (34)–(37) and simplifying
032409-6
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
Eq. (37) can be written as 2 2
∂y − k − κi2 φˆi = 0.
PHYSICAL REVIEW E 89, 032409 (2014)
φˆ 1 = (e−ky − ek(y−2β) )c5 , (39)
Equations (39) and (40) are the governing equations of the system. The solutions to these are sought with appropriate boundary conditions in terms of perturbation variables. Expressing the boundary conditions listed in Sec. (II C) in terms of perturbations in the normal mode form, the no slip boundary condition at the interface becomes vˆ1 = vˆ2 ,
(41)
uˆ 1 = uˆ 2 .
(42)
The normal stress balance and the tangential stress balance upon substituting the stress definitions from Eqs. (13) and (14) become [[−pˆi + 2Ci ∂y vˆi + i ∂y φmi ∂y φˆ i ]] = −Bk 2 hˆ 2 ˆ y φmi = 0, (43) Ci (ik vˆi + ∂y uˆ i ) + i ik φˆ i ∂y φmi + i ik h∂ where B is a nondimensional number given by B = continuity of potential is given by
γ h0 . 0 φ02
ˆ y φmi ]] = 0, [[φˆ i + h∂
The (44)
where the second term indicates contribution because of the gradient in the mean potential. The normal component of the electric field, balanced across the interface, is written as [[−i ∂y φˆ i ]] = 0.
(45)
The dynamics of the position of the interface and the interfacial charge is given by the kinematic condition and the charge conservation equation, respectively. The kinematic condition becomes ∂t hˆ = vˆ2 .
(46)
A. Dispersion equation
For the case of time-invariant fields, the time dependence of the perturbation amplitude is assumed to be of the form est , where s is the growth rate of the perturbation. The sign of s decides the stability of the system. A negative s indicates a stable system, whereas the system is rendered unstable if s is positive. The second- and fourth-order ordinary differential equations (ODEs) in terms of φˆ i [Eq. (39)] and vˆi [Eq. (40)] are solved using the following boundary conditions: y → β or −1 vˆi = 0, ∂y vˆi = 0, and φˆ i = 0, respectively. The expressions obtained for vˆi and φˆ i are − (y + 2kyβ − 2kβ )c2 )], + [(−1 + e2k(1+y) )y + 2k(1 + y)c4 ],
(50)
V. RESULTS AND DISCUSSION
The objective of the present work is to overcome the limitation of electroneutrality of the bulk fluid imposed in the EHD models by accounting for the thermal motion of ions in the bulk. The electrokinetic model (EK) where the electric potential is modeled using the Debye-Huckel theory is used to overcome this limitation. Thereby, we assume that the dynamics of the ions is quite fast and the ion distribution is Boltzmann and the voltage or potential is low. The results from the EK theory are compared to those from the EHD model. The dynamics of the system when subjected to infinitesimal random perturbations can be gauged by determining the growth rate or the growth exponent s. The dispersion relation obtained from the linear stability analysis discussed in Sec. IV A is solved to obtain an expression for s the growth rate of the perturbation. This expression is plotted for different values of κ in Fig. 5. It can be seen from the figure that the growth rate of the perturbation depends on its wave number and the inverse Debye layer thickness κ and exhibits a maxima (smax ) for a particular value of wave number kmax which is the signature of the instability. An increase in κ results in 0.1
κ=0.01 κ=1
(47)
(48)
)c6 ,
ˆ where C is the column matrix of the constants c1 –c6 and h. This system of linear equations can have a nontrivial solution only when the coefficient matrix is singular. Therefore, equating the determinant of matrix A to zero we get a characteristic equation in terms of s and k and other system parameters like B, β, μr , 1 , 2 , σ1 , and σ2 . This equation is then solved using MATHEMATICA 7.0 for s, the growth rate of the perturbation.
0.05
κ=10 κ=100
0 −0.05 −0.1 0
vˆ2 = e−k(2+y) [(−1 + e2k(1+y) − 2k(1 + y)]c3
(−1 + e
(49)
AC = 0,
vˆ1 = e−ky [e2ky (c1 + yc2 ) + e2kβ ((−1 − 2ky + 2kβ)c1 2
k 2 +κ 2
√ 2(1+y) k 2 +κ 2
uˆ i and pˆ i are obtained in terms of vˆi from Eqs. (34), (35), and (36). These expressions for vˆi , uˆ i , pˆ i , and φˆ i , in terms of the constants c1 –c6 are substituted in the eight boundary conditions to get a system of homogeneous linear equations ˆ These equations in terms of the constants and the unknown h. can be represented by the matrix equation
s (growth rate)
Similarly, upon eliminating pˆ from Eqs. (35) and (36), we get a fourth-order differential equation in v which is given by 2
2 ∂y − k 2 vˆi = 0. (40)
φˆ 2 = −e−y
√
0.2
0.4 0.6 k (wavenumber)
0.8
1
FIG. 5. Variation of the growth rate s as a function of the wave number k at different values of κ. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
032409-7
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014) TABLE I. Debye layer thicknesses for a given conductivity of the fluid.
0.7
k
max
0.6
S. no
0.5 0.4 0.3 −4 10
−3
10
−2
10
−1
10
κ
1
10
2
10
3
10
4
10
S. no
FIG. 6. Variation of kmax as a function of the inverse Debye length κ. The dotted line indicates results obtained from the PD-PD theory, the dashed line indicated those from the PD-PC, while the solid line indicates results from the electrokinetic theory. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
an increase in the growth rate due to a dominance of free charge. In Figs. 6 and 7, kmax and smax have been plotted as a function of κ. A similar trend as seen with the base state variables is seen with the instability parameters. The values of kmax and smax plateau to limiting values of PD-PD and PD-PC theories at low and high values of κ, respectively. For the case of an air-water system, where the conductivity of pure deionized water is taken to be 5.5 × 10−6 Sm−1 , the Debye length is obtained using the relationship σ = 0 Dκ 2 , where the diffusion coefficient D = 2.8 × 10−9 m2 s−1 is λD = κ −1∗ = ∼9 × 10−7 m. Typical electrode spacing in EHD-based soft lithography processes vary from 100 to 1000 nm [5,29,30] with a few groups reporting use of 80 and 100 μm [31]. Under these conditions, the thickness of the Debye layer is almost of the order of the electrode spacing. For other fluids like silicone oil and polymer melts with very low conductivities ranging from 10−9 to 10−12 S/m, the Debye length will further increase, although a low diffusion coefficient characteristic of these fluids ensures that it ranges from the submicron range to being of the order of electrode spacing. Therefore, an EK model is very necessary to describe the dynamics of the diffuse layers in these fluid-fluid systems.
0.06 0.05
s
max
0.04 0.03 0.02 0.01 0 −4 10
−3
10
−2
10
−1
10
κ
1
10
2
10
3
10
1 2 3 4 5
1 2 3 4 5
Conductivity (H2 O) S/m
Debye length (H2 O) m
1.0225×10−2 1.0225×10−3 1.0225×10−4 1.0225×10−5 1.0225×10−6
8.69293 × 10−9 2.748 × 10−8 8.69293 × 10−8 2.74894 × 10−7 8.69293 × 10−7
Conductivity (PDMS) S/m
Debye length (PDMS) m
−8
1.0225×10 1.0225×10−9 1.0225×10−10 1.0225×10−11 1.0225×10−12
4.364 × 10−8 1.374 × 10−7 4.364 × 10−7 1.374 × 10−6 4.364 × 10−6
VI. CONCLUSIONS
As one goes to lower thicknesses, the interfacial instability in thin fluid films under an electric field can be best described by an electrokinetic model which takes into account the dynamics of diffuse layers formed at the surfaces in the system. It is shown that the wavelength of the instability depends significantly on the Debye length, which is comparable to the fluid thickness for the low-conductivity fluids used in these studies. The PD-PD and PD-PC models are only limiting cases when the Debye length is very large and very small, respectively, when compared to film thickness. Thus it can be summarized that, in a perfect PD-PD system, free charges are absent and, hence, the electric body force term in the momentum balance equation is negligible. Even at the interface, free charges are absent, hence, the electric displacement continuity equation is used as a boundary condition. For a PD-PC system, with a large conductivity of the PC fluid, under a dc field the charges present in the conducting fluid move instantaneously to the surface of the fluid and reorganize into small Debye layers whose thickness is of the order of nanometers. The field within the fluid is nullified and the surface becomes an isopotential. The bulk of the fluid therefore can be considered as electroneutral. In case of leaky dielectrics, for conductivities of the order of 10−11 S/m, the Debye layer thickness is of the order of 1 to 5 μm, especially for viscous oils. When the thickness of the fluid film is much larger than a micron, the system is essentially a charged Debye layer separated by a core layer which is electrically neutral. An estimate of Debye layer thicknesses is given in Table I for typical aqueous and nonaqueous systems. Thus one can determine when such a system would behave as a perfect or a leaky dielectric (conductor) if the confinement is known.
4
10
FIG. 7. Variation of smax as a function of the inverse Debye length κ. The dotted line indicates results obtained from the PD-PD theory, the dashed line indicates those from the PD-PC, while the solid line indicates results from the electrokinetic theory. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
ACKNOWLEDGMENTS
The authors acknowledge the financial assistance provided by the Department of Science and Technology, India. The authors also acknowledge useful comments provided by V. Kumaran (IISc, Bangalore) and the reviewer.
032409-8
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
APPENDIX: RELATIONSHIP BETWEEN THE INVERSE DEBYE LENGTH κ AND CONDUCTIVITY σ2
Using the DH approximation, −e[φ(y) − φb ] σ = e m+ n0 1 − kB T e[φ(y) − φb ] , + m− n0 1 + kB T
To be able to compare the electrokinetic model to the electrohydrodynamic model a relationship between the concentration of ions in the liquid and the conductivity of the liquid has to be determined. The electric conductivity (σ ) is related to the cationic and anionic molar concentrations c± by the following relation [22]: ˆ + c+ + m ˆ − c− ), σ = F (m 2
Assuming m+ = m− = m, −e[φ(y) − φb ] +1 = emn0 1 − kB T e[φ(y) − φb ] + , kB T
(A1)
where F is the Faraday’s constant given by F = eNa , where ˆ± e is the elementary charge and Na is Avagadro’s number. m are the mobilities of ions (units: mol N−1 ms−1 ) and c± are the molar concentrations, respectively, ˆ + c+ + m ˆ − c− ), σ = e2 Na2 (m ˆ + n+ + m ˆ − n− ), = e Na (m 2
= e(m+ n+ + m− n− ),
(A2)
where m± are the ionic mobilities and n± are the number concentration of charges, respectively. Substituting the expression for charge concentration from the equation below,
= 2emn0 .
2e2 n0 0 kB T and the Einstein’s relation D = mkB T /e, where D is the diffusivity (m2 s−1 ),
Using the expression for inverse Debye length κ 2 =
σ = D0 κ 2 .
∓e[φ(y)−φb ] kB T
n± = n0 e , −e[φ(y)−φb ] e[φ(y)−φb ] ⇒ σ = e μ+ n0 e kB T + μ− n0 e kB T .
(A4)
(A5)
(A3)
Therefore, the conductivity σ depends on the diffusivity, the dielectric constant, and the square of the Debye length.
[1] D. A. Saville, Annu. Rev. Fluid Mech. 29, 27 (1997). [2] S. Torza, R. G. Cox, and S. G. Mason, Philos. T. R. Soc. A 269, 295 (1971). [3] O. Vizika and D. A. Saville, J. Fluid Mech. 239, 1 (1992). [4] G. I. Taylor, Proc. R. Soc. London A 313, 453 (1969). [5] E. Schaffer, T. Thurn-Albrecht, T. P. Russell, and U. Steiner, Nature 403, 874 (2000). [6] G. I. Taylor and A. D. McEwan, J. Fluid Mech. 22, 1 (1965). [7] E. B. Devitt and J. R. Melcher, Phys. Fluids 8, 1193 (1965). [8] D. H. Michael, Proc. Camb. Philos. Soc. 64, 1203 (1968). [9] J. R. Melcher and C. V. Smith, Phys. Fluids 12, 778 (1969). [10] E. Schaffer, T. T. Albrecht, T. P. Russell, and U. Steiner, Europhys. Lett. 53, 518 (2001). [11] Z. Lin, T. Kerle, S. M. Baker, D. A. Hoagland, E. Schaffer, U. Steiner, and T. P. Russell, J. Chem. Phys. 114, 2377 (2001). [12] P. Gambhire and R. M. Thaokar, Phys. Fluids 22, 064103 (2010). [13] P. Gambhire and R. M. Thaokar, Eur. Phys. J. E 34, 84 (2011). [14] P. Gambhire and R. M. Thaokar, Phys. Rev. E 86, 036301 (2012). [15] W. B. Russel, D. A. Saville, and W. R. Scholwalter, Colloidal Dispersions (Cambridge University Press, Cambridge, 1989). [16] V. Kumaran, Phys. Rev. Lett. 85, 4996 (2000). [17] V. Kumaran, Phys. Rev. E 64, 011911 (2001). [18] D. Lacoste, G. I. Menon, M. Z. Bazant, and J. F. Joanny, Eur. Phys. J. E 28, 243 (2009).
[19] F. Ziebert, M. Z. Bazant, and D. Lacoste, Phys. Rev. E 81, 031912 (2010). [20] W. Qu and D. Li, J. Colloid Interface Sci. 224, 397 (2000). [21] C. K. Hua, D. W. Lee, and I. S. Kang, Colloid. Surface. A 372, 86 (2010). [22] E. K. Zholkovskij, J. A. Masliyah, and J. Czarnecki, J. Fluid Mech. 472, 1 (2002). [23] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevE.89.032409 for the derivation of the Debye-Huckel equation and calculation of ion concentration. [24] H. N. W. Lekkerkerker, Physica A 167, 384 (1990). [25] R. E. Goldstein, A. I. Pesci, and V. Romero-Rochin, Phys. Rev. A 41, 5504 (1990). [26] P. M. Vlahovska, K. D. Danov, A. Mehreteab, and G. Broze, J. Colloid Interface Sci. 192, 194 (1997). [27] A. Ramos, Electrokinetics and Electrohydrodynamics in Microsystems (Springer, New York, 2011). [28] M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, New J. Phys. 11, 075016 (2009). [29] M. D. Morariu, Z. Lin, T. P. Russell, and U. Steiner, Nature 2, 48 (2002). [30] P. Goldberg-Oppenheimer and U. Steiner, Small 6, 1248 (2010). [31] N. Arun, A. Sharma, P. S. G. Pattader, I. Banerjee, H. M. Dixit, and K. S. Narayan, Phys. Rev. Lett. 102, 254502 (2009).
032409-9
Electrokinetic model for electric-field-induced interfacial instabilities Priya Gambhire and Rochish Thaokar* Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai 400 076, India (Received 18 October 2013; revised manuscript received 19 December 2013; published 26 March 2014) Technology based on electric-field-induced instabilities on thin polymer film surfaces has emerged as a promising candidate for soft lithography. Typically, the instability is modeled using the perfect dielectric (PD) or the leaky dielectric (LD) model. These assume the electric diffuse layer to be infinitesimally large or small, respectively. In the present work we conduct stability analysis assuming a PD-electrolyte solution interface. The concentration of ions and, hence, the diffuse layer thickness is in general assumed to be of the same order as the electrolyte film thickness. The PD-LD models are then realized as limiting cases of the ratio of the double layer thickness to the film thickness. DOI: 10.1103/PhysRevE.89.032409
PACS number(s): 68.03.−g, 47.57.jd, 47.65.−d
I. INTRODUCTION
An interface between two fluids destabilizes under electric fields in the presence of a contrast in the electrical properties of the two fluids such as dielectric constant and electrical conductivity. This phenomenon has been termed electrohydrodynamic instability (EHDI) [1] and finds applications in techniques such as electrospraying (drops under electric field) [2,3], electrospinning (cylindrical jets under electric field) [4], electropatterning (or soft lithography, planar fluid interfaces under electric fields) [5], and so on. Pioneering work in this field was conducted in the 1960s with contributions from G. I. Taylor and J. R. Melcher [6–9], who studied various aspects of EHD such as drops and planar fluid deformation under electric field. While dealing with instabilities in dielectric materials, the leaky dielectric (LD) model was proposed [1,5,10–13]. It overcame the fallacies in the instability predictions obtained assuming either perfectly dielectric (PD) or perfectly conducting (PC) fluids. According to the leaky dielectric theory, the fluids are assumed to have an infinitesimal amount of free charge confined to the interface while the bulk fluid is free of charge. Moreover, the thermal motion of ions is neglected. This model can explain deformation patterns which are beyond the scope of the perfect dielectric or the perfect conductor theory, increasing its utility in the field of electrohydrodynamics. In 2000, following the work by Schaffer et al. [5] wherein a submicron thin polymer film is observed to deform into ordered patterns under electric field, the leaky dielectric model is modified to suit these thin-film systems (negligible gravity). The leaky dielectric model combined with a linear stability analysis yields a wavelength of the instability which is in good agreement with the literature. However, it has been recently shown that in certain parametric regimes, the assumptions made in the leaky dielectric model may be invalid [14]. The wavelength of the instability depends on parameters such as conductivity, applied electric field, and frequency which alter the inherent time scales in the system. Nonlinear analysis is used to seamlessly obtain wavelengths at all the ratios of time scales under both dc and ac fields.
*
[email protected]
1539-3755/2014/89(3)/032409(9)
Although general in approach, the above studies neglect the thermal motion of the free charge present in the system. In a real system, free ions, albeit infinitesimal, will still exhibit diffusion leading to a diffuse layer of counterions around charged surfaces or interfaces. This layer is characterized by its thickness, called the Debye length, κ −1 . The Poisson-Nernst-Planck equation governs the dynamics of the charges within electrolytic materials. This equation leads to the Boltzmann distribution at equilibrium when charge dynamics is unimportant. For the case of the potential being very less than kB T , the Debye-Huckel approximation is appropriate and has been used quite extensively in the field of electrokinetics to study membrane or fluid interfaces [15]. The coupling of the curvature of a bilayer membrane to its surface charges where the charges can flip to the other side and where they can diffuse on the same side of an interface results in a curvature-surface charge density coupling that leads to instabilities [16]. The Poisson-Boltzmann equation with the Debye-Huckel approximation is used to obtain an expression for the electrostatic energy of the system and different modes of instabilities are analyzed. A linear stability analysis of the electrokinetic model [17] of a floating bilayer with surface charges on either side exhibits an instability. The presence of the bilayer membrane adds a dielectric layer between two conducting layers of fluid. A diffusion equation for the surface charges enables a coupling between the surface undulations of the membrane and the charges on the surface which beyond a threshold lead to an instability. In another set of studies, the effect of the Debye layers on the electrostatic properties of an artificial cell membrane have been investigated [18]. In the limit of vanishing thickness and dielectric constant and in the absence of inherent charges on the membrane, the system reduces to that of a fluid-fluid interface. The system differs from the present one in that ion channels are assumed and a discontinuity is imposed on the electric field at the interface. The authors study the effect of a perturbation on this system in order to obtain the corrections to the elastic modulii. They also characterize the induced-charge electro-osmotic flow around the perturbed membrane. In a follow-up study [19], the authors proposed an improved zero-thickness electrokinetic model under dc fields. Using the Poisson-Boltzmann equation with the Debye-Huckel approximation in conjunction with Robin-type of boundary
032409-1
©2014 American Physical Society
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
conditions, they obtain a dispersion relation for the perturbed membrane along with expressions for the threshold voltages for the instability in the case of a conductive and a nonconductive membrane. They present this model as a basis for further development to systems of various sizes and charge densities. The zero-thickness membrane model discussed above differs from the present work in a couple of aspects. The membrane is conducting while its distance from the electrodes is quite large compared to the Debye layer thickness. Thereby a quasineutral environment is assumed. The bending modulus is suppressed if the electrode spacing is decreased to the order of microns. The objective of the present work is to study the instability at a fluid-liquid interface as a function of the thickness of the liquid film. As the thickness changes, the assumption of electroneutrality in the bulk of the fluid (that is, the fluid away from the Debye layers) becomes invalid. Therefore, the main objective of this work is to study the scenarios of overlapped Debye layers. A detailed study of system with overlapped Debye layers is reported by Qu and Li [20]. They study a fluid confined between two flat plates such that the double layers formed on the plates overlap. They argue that the Boltzmann approximation cannot be used for a system with overlapped Debye layers as this violates the assumption of electroneutrality in the bulk liquid. Using suitable boundary conditions they develop governing equations for this case. They show that at low plate distances the predictions from the classical theory differ substantially from those predicted by their model. Although they study the system of overlapped Debye layers, the system consists of a single fluid under a field. A study of the limits where the electrokinetic model might match the predictions from the electrohydrodynamic model is absent. Electrokinetic models vis-`a-vis electrohydrodynamic models have been analyzed in the study of drop deformation. Hua et al. [21] use the electrokinetic model in the context of conducting drops suspended in a dielectric medium. They argue that these systems hitherto studied via a conductor model can be approached using an electrokinetic model. They motivate the use of this model in systems with high curvature such as small droplets. They show in this system the conductor model results as a limiting case (for thin double layers) to the electrokinetic model. Zholkovskij et al. [22] studied the deformation of drops under electric field taking into account the presence of the double layers. In this generalized model, termed the electrokinetic (EK) model, the thickness of the Debye layers occurs as a parameter which affects the deformation of the drops. They show that at very low Debye layer thicknesses, the drop deformation matches that predicted by the electrohydrodynamic theory assuming perfectly conducting materials while at very high Debye thicknesses it matches that of the perfect dielectric theory. In the present work, the system under study is that of two stratified planar fluids. Under the action of an applied electric field, the free charges in the system migrate and accumulate at the interface and on the electrode surfaces. Under the combined action of the electric field and thermal motion, these charges are assumed to form diffuse layers at these surfaces. The field and the instability in such a system is
studied using the “electrokinetic model.” The objective of the present work is to study the effect of confinement of the fluid films to thicknesses as those practically used in experiments such as by Schaffer et al. [5]. The confinement of the liquid films result in the overlapping of the Debye layers formed within the fluid. This in turn violates the electroneutrality condition. The concentration and the potential within the bulk solution are not considered zero but are derived consistently (refer Appendix). To our knowledge, a thorough study of how the EK model predicts various regimes resulting from the steadily increasing confinement for a simple fluid-fluid interface has not been carried out. Apart from carrying out the above, a systematic evidence of how the EK model matches the predictions made by the leaky dielectric model in the limiting cases gives universality to the model studied. II. THE MODEL A. System description
The system consists of two fluid layers held between planar electrodes (Fig. 1). The longitudinal direction is denoted by x while the voltage is applied along the transverse direction, i.e., the y direction. Variations in the third direction (z) are not considered. The upper fluid is referred to as 1 and the lower fluid as 2. The viscosity, density, conductivity, and dielectric constant of the upper fluid are denoted by μ1 , ρ1 , σ1 , and 1 , respectively, and those of the lower fluid by μ2 , ρ2 , σ2 , and 2 . The position of the interface is denoted by y ∗ = 0 while the electrodes are held at y ∗ = βh0 and y ∗ = −h0 . The superscript ∗ denotes dimensional quantities. B. Governing equations
The continuity and the momentum balance equations are used to model the hydrodynamics in the system. In the dimensional form, in the presence of an applied electric field they can be written as ∇ ∗ · V ∗i = 0,
(1)
ρ(∂t∗ V ∗i + V ∗i · ∇ ∗ V ∗i ) = −∇ ∗ pi∗ + μi ∇ ∗2 V ∗i − ρf∗i ∇ ∗ φi∗ , (2)
FIG. 1. A schematic of the interface between two dielectric fluids subjected to a normal electric field. The dashed line is a schematic representation of the potential distribution in the two fluids.
032409-2
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
where the subscript i = 1,2 refers to fluid 1 and fluid 2, respectively; pi is the pressure; Vi = ueˆ x + v eˆ y is the velocity vector; and ρ and μ are the density and the viscosity of the fluids, respectively. The final term on the right-hand side of the momentum balance equation gives the contribution of the applied electric field to fluid flow. ρfi and φi are the free charge density and the electric potential in the fluid, respectively. When electroneutrality is assumed to be valid in the bulk solution, this term becomes negligible and the electrohydrodynamic coupling comes only through interfacial conditions. To derive the equation governing the charges in the system, the species conservation equation can be written as ∂n∗± = −∇ ∗ · (∓em± n∗± ∇ ∗ φ ∗ − m± kB T ∇ ∗ n∗± ), ∂t ∗
(3)
where e is the elementary charge; n∗+ ,n∗− are the number concentrations of the free charges in the fluid; and m± is the mobility. Subtracting the positive and negative number concentrations, and using ρf∗ = e(n∗+ − n∗− ),σ = e(m∗+ n∗+ + m∗− n∗− ),
(4)
potential, at a reference position in the fluid y = b, and kB and T are the Boltzmann constant and temperature, respectively. The details of the derivation of Eq. (8) are given in the Supplementary Material [23]. It also shows that for viscous time scales, the dynamics of the charges can be ignored, and, hence, the Poisson Boltzmann equation for potential is decoupled. For the DH approximation to be valid, the term eφ ∗ /(kB T ) needs typically φ = 25 mV for it to be O(1). The voltages applied in EHD-based soft lithography are typically much higher, and ideally the Poisson-Boltzmann (PB) equations should be used. However, in the past, the DH model has been used as a substitute to the more complicated, nonlinear PB model to understand the physics analytically, which is also the endeavor in the present study. For further reading on the DH approximation and solving the nonlinear PB equation the authors cite the work by Lekkerkerker [24], Goldstein et al. [25], Vlahovska et al. [26]. C. Boundary conditions
The following set of boundary conditions are used to solve the governing equations. At the interface, the velocities are continuous in the normal and tangential directions,
results in the following equation: ∂ ∗ ρf∗ ∂t ∗
= −∇ ∗ · (σ ∇ ∗ φ ∗ − D∇ ∗ ρf∗ ).
(5)
process time scale τp =
μ2 h20 . 0 φ02 2
(6) 0 σ
and the
For a typical air-silicone oil
system, σ = 10−12 C2 /Nm s, 0 = 8.85 × 10−12 C2 /Nm2 , h0 = 100 × 10−6 m, μ2 = 30 Pa s, and φ0 = 30 V. Using these values τc ∼ 1 s and τp ∼ 37.6 s. For an air-water system, with the typical values of σ = 5 × 10−6 C2 /Nm2 s, 0 = 8.85 × 10−12 C2 /Nm2 , h0 = 100 × 10−6 m, μ2 = 10−3 Pa s, and φ0 = 30 V, the values of τc and τp are 1.7 × 10−6 s and 1.26 × 10−3 , respectively. Therefore, the left-hand sides of Eqs. (3) and (6) are equivalent to zero and the electrostatic equation therefore reduces to the Poisson equation [23], 0 i ∇ ∗2 φi∗ = −ρf∗i ,
(V1 ∗ · tˆ∗ ) = (V2 ∗ · tˆ∗ ),
(10)
1+∂x h
(8)
where κi∗ is the length. For fluid 1, κ1∗ = 0, while for √ Debye 2e2 no ∗ fluid 2, κ2 = 0 2 kB T , n0 , and φb are the uniform number concentration of the ions (in the absence of field) and the
1+∂x h
subscripts x, y, or t denotes differentiation with respect to that variable, h∗ (x,t) denotes the position of the interface, and V1 ∗ and V2 ∗ are the velocity vectors for the two fluids. For an unperturbed interface nˆ = eˆ y and tˆ = eˆ x . The stresses are balanced in the normal and tangential directions. In the normal direction the stresses satisfy [[nˆ ∗ · τi∗ · nˆ ∗ ]] = γ C ∗ ,
(11)
where the operator [[Xi ]] denotes the jump X1 –X2 across the interface h(x,t). C ∗ is the curvature given by C ∗ = ∇ ∗ · nˆ ∗ = −∂x∗2 h∗ /(1 + ∂x∗ h∗2 )3/2 . The term on the right-hand side in Eq. (11) represents the stabilizing force due to surface tension. The balance of stresses in the tangential direction is given by [[ tˆ∗ · τi ∗ · nˆ ∗ ]] = 0,
(7)
where 0 is the permittivity of vacuum and i is the dielectric constant of the two fluid layers i = 1,2, respectively. For the case of an air-liquid interface, ρfi , the free charge density in the top fluid is 0 while that in the lower liquid is given by ρf2 = e(n∗+ − n∗− ). This equation under the Debye-Huckel (DH) approximation reduces to ∇ ∗2 φi∗ = κi∗2 (φi∗ − φb∗ ),
(9)
where nˆ ∗ and tˆ∗ are the unit normal and the unit tangential vectors, respectively, as indicated in Fig. 1, given by −∂ ∗ h∗ eˆ +eˆ eˆ +∂ ∗ h∗ eˆ nˆ ∗ = √x x∗ ∗2y and tˆ∗ = √x x ∗ ∗2y . The operator ∂ with the
Upon scaling, the equation reduces to τc ∂ρf 0 D = −∇ · ∇φ − ∇ρ f , τp ∂t σ h20 where the charge relaxation time scale τc =
(V1 ∗ · nˆ ∗ ) = (V2 ∗ · nˆ ∗ ),
(12)
where τi ∗ , the total stress, i.e., the sum of the hydrodynamic and the electrical stresses, for the fluid i, is given by τi ∗ = −pi ∗ I + μi (∇ ∗ Vi ∗ + ∇ ∗ Vi ∗T ) + Mi ∗ .
(13)
The superscript T indicates transpose, i = 1,2 denotes the two fluid layers and Mi is the Maxwell stress tensor [27], Mi ∗ = i 0 E ∗i E ∗i − 12 (E ∗i · E ∗i )I . (14) A balance of the tangential component of the field gives the continuity of the potentials across the interface
032409-3
φ1∗ = φ2∗
(15)
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
and the balance of the normal component of electric field gives the electric displacement continuity equation [[0 i (−∇ ∗ φi∗ · nˆ ∗ )]] = 0.
(16)
The boundary condition assumes no adsorbed charges at the interface. It should be noted that although the bulk charge is nonzero the surface charge is zero. In fact, the net bulk charge integrated over the double layer (in the electrokinetic model) is really the surface charge integrated over the area (in the perfect dielectric-leaky dielectric model). This clearly differs from the free charge BC in PD-LD systems, where [[0 i (−∇ ∗ φi∗ · nˆ ∗ )]] = q aided by the equation ∂q = [[σi∗ ∇ ∗ φi∗ · nˆ ∗ ]]. The ∂t dynamics of the position of the interface is given by the kinematic condition ∂t h∗ + V1 ∗ · ∇ ∗ s h∗ = v2∗ = v1∗ .
steady state from the perturbed state. Any random function can be simply expressed as a sum of Fourier modes. As in the linear limit only O(1) terms need to be considered, LSA gives a neat, simplified mathematical formalism to study the stability of any system. In the following sections the methodology of LSA applied to the present system is described. A. Addition of perturbations
The system is subjected to an infinitesimal arbitrary perturbation. The system variables are therefore expressed as a sum of a mean value (denoted by an overbar on the variable) and a perturbation term (denoted by the superscript tilde). Thus, ¯ ˜ m(x,y,t) = m(x,y,t) + ηm,
(17)
D. Scalings
The parameters in the above equations are scaled using their characteristic values in the system. That is, length, in both the x and y directions, is scaled using the thickness of the thin film h0 . Potential is scaled by φ0 , the applied potential. Velocity is scaled using the term 0 φ02 /(μ2 h0 ) and pressure is scaled using 0 φ02 / h20 while the scaling for time is obtained to be μ2 h20 /(0 φ02 ).
where η is a small variable which gives the deviation from the base state. The perturbations are incorporated in the governing equations. The leading-order terms give the base state of the system while order η variables are considered in the linear stability analysis. 1. Base state
In the present case, in the base state the interface is flat with no perturbations and, hence, the mean velocities are zero,
E. Governing equations
The equations, after incorporating the dimensionless variables listed in Sec. II D are as follows: ∇ · V i = 0,
(18)
umi = 0, vmi = 0.
pm1 = 0,
(19)
−∇pi + μr ∇ 2 V i + i ∇ 2 φi ∇φi = 0.
(20)
pm2 =
∇ φi = κ (φi − φb ), 2
1 2 2
∂φ ∂y
(24)
2 + Pc .
Note that the hydrostatic equations for PD and LD cases is ∂y pi = 0 while for the EK model it is ∂y p = ρf E = κ 2 φ∇φ. The mean potential or the base-state potential is calculated by solving the following equation (see the Appendix for derivation) for fluids 1 and 2: ∂y2 φmi = κi2 (φmi − φb ).
(21)
(25)
The constant Pc is calculated from a stress balance such that ∂φm2 2 ∂φm1 2 1 pm2 = − 1 . (26) 2 2 ∂y ∂y
The nondimensional Poisson equation is 2
(23)
The base-state pressure is obtained from the (y-direction) momentum balance equation in both fluids,
Rei (∂t V i + V i · ∇V i ) = −∇pi + μr ∇ 2 V i + i ∇ 2 φi ∇φi , where μr = μ1 /μ2 (note: the nondimensional variables are denoted without the superscript ∗). The Reynold’s number (Re) is given by Rei = ρi 0 φ02 /μ22 for dc fields. For a typical system of air (fluid 1) and silicone oil (fluid 2) the properties are ρ1 ∼ 1 kg/m3 , ρ2 ∼ 970 kg/m3 , μ2 ∼ 0.97 Pa s, and h0 ∼ 100 × 10−9 m and, using φ0 ∼ 30 V, the Reynolds number for the dc case is of the order of 10−8 –10−5 . Hence, in the ensuing analysis, we consider the case where Rei are approximately zero. The momentum equations with the above parameter values are simplified to
(22)
(27)
where κ is the nondimensional inverse Debye layer thickness κ = κ ∗ h0 .
The following (scaled) boundary conditions, which include the potential at the electrodes and the continuity of potential and electric displacement at the interface, are used to solve the above equation:
III. LINEAR STABILITY ANALYSIS
φm1 |y=β = 0,
(28)
The linear stability analysis (LSA) is an important tool to study the response of a system to random disturbances. The basic principle of this theory is that any system always undergoes infinitesimal perturbations in its steady or equilibrium state. A system can be called linearly stable if it returns to its
φm2 |y=−1 = 1,
(29)
φm1 = φm2 |y=0 ,
(30)
[[−i ∂y φmi ]] = 0.
(31)
032409-4
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
The reference potential φb is determined using the condition 0 0 ρ2 dy = 2 κ22 (φm2 − φb )dy = 0. −1
−1
The base state values of the potential in the two fluids assuming κ1 = 0 and κ2 = κ are given by (1 + eκ )(y − β)2 κ , (32) 2(−1 + eκ )1 + (1 + eκ )β2 κ
β2 κcosh(κ/2) + 1 sinh(κ/2) − sinh 12 + y κ . = β2 κcosh(κ/2) + 21 sinh(κ/2) (33) φm1 = −
φm2
In Fig. 2, the potential distribution in the fluids in the base state is plotted along the normal direction y. The top fluid is bounded at y = 1 and y = 0 to −1 denotes the lower fluid. The top fluid is assumed to be air (a perfect dielectric) while the potential in the lower fluid is studied as a function of the concentration of ions through the parameter, inverse Debye length κ. The origin of the Debye layer in this case is due to the charge separation under electric field. Ions are attracted to respective electrodes while the counterions are repelled and accumulate at the interface. In the presence of thermal forces, these charges do not exist in a thin layer at the interface but form a diffuse layer with the concentration of ions decaying away from the interface or electrode. This leads to a Debye layer formation at both the interface and the electrode in the lower fluid layer giving rise to a potential distribution that is dependent upon the concentration of ions in the layer. As discussed earlier, increase in the concentration of ions leads to the formation of thinner, mroe compact Debye layers while lower concentration of ions leads to thicker, more diffuse Debye layers. In Fig. 2, the potential distribution for κ = 10−3 (solid line), κ = 5 (dashed line), κ = 10 (dotted line), and κ = 103 (dash-dot line) progressively representing thinner
bilayers, in the bottom fluid, is plotted along y. In the top fluid (y = 0 − 1), the potential distribution is always linear 32, a characteristic of the potential drop across a dielectric medium which is governed by the Laplacian of the potential being zero, while in the lower fluid (y = 0 to −1), the potential drop is given by Eq. (33) which shows the characteristic exponential drop within the Debye layers at the interface (y = 0) and the electrode surface (y = −1). The value κ = 1 denotes Debye layers that are as thick as the film thickness. In such a situation, the Debye layers below the interface and the layer above the electrode surface overlap. At low concentrations of ions, the overlap of potentials from the two Debye layers is significant to the extent that the potential drop becomes linear and the system behaves like a dielectric material until the concentration is so low that the potential drop in the system is equivalent to that in a PD material [refer Fig. 3]. As κ increases beyond 1, Debye layers become thinner, the overlap decreases, and the field at midplane of the fluid layer becomes zero. With increase in the concentration of ions (conductivity) or κ the Debye layers asymptote to zero thicknesses and, correspondingly, show perfect conductor-like behavior. The intricate physics of compact Debye layers which are referred to as Stern layers is beyond the scope of this paper and more details can be found in Refs. [19,28]. To get a perspective of the practicality of the κ values used in the study consider if the film thickness of 100 μm is a fixed parameter also typically encountered in experiments, the range of κ studied in the present work corresponds to the Debye layer thicknesses of 0.1 μm to 0.1 m. Although the value of 0.1 m appears unrealistic from an electrokinetics point-of-view the high thickness of the Debye layer means that the fluid has a very low concentration of charge spread over a large distance and the characteristics tend to that of a perfect dielectric material. A Debye layer thickness of 0.1 μm, on the other hand, is quite realistic. To add another perspective to the term κ. For a fixed fluid-fluid system, say, air-silicone oil, having thinner fluid film changes the wavelength of the EHD instability [5]. The silicone
1 0
0.5 y
−0.2
0
y
−0.4 −0.6
−0.5
−0.8
−1 0
0.2
0.4
φ
m
0.6
0.8
1 −1 0.75
FIG. 2. Potential distribution (φm ) in the two fluid layers [y = 1 to 0 indicates the top fluid (1) and y = 0 to −1 indicates the lower fluid (2)] in the base state at different values of κ. The solid line indicates κ = 10−3 , the dashed line indicates κ = 5, the dotted line indicates κ = 10 and the dash-dot line indicates κ = 1000. The potential distributions obtained from the PD-PD and the PD-PC theories coincide exactly with the κ = 10−3 and κ = 1000 curves, respectively, and, hence, cannot be distinguished. The values of parameters used are 1 = 1, 2 = 3, σ1 = 0 (implies κ1 = 0), μr = 0, and β = 1.
0.8
0.85
φ
m2
0.9
0.95
1
FIG. 3. Potential distribution (φm2 ) in the lower fluid layer, in the base state, at different values of κ. The dashed line indicates κ = 10−4 , the dotted line indicates κ = 1, the dash-dot line curve indicates κ = 5, the curve with ◦ legends indicates κ = 10, the solid line curve indicates κ = 102 , and the curve with the + legends indicates κ = 103 . The potential distributions obtained from the PD-PD and the PD-PC theories coincide exactly with the κ = 10−4 and κ = 1000 curves, respectively, and, hence, have not be plotted. The values of parameters used are 1 = 1, 2 = 3, μr = 0, and β = 1.
032409-5
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014)
IV. STABILITY ANALYSIS
The equations governing the perturbation variables [obtained by collecting the order η terms from Eqs. (18), (20), and (21)] after substituting the normal mode form for the variables are (34) ∂x u˜ i + ∂y v˜i = 0, 2
−∂x p˜i + Ci ∂x u˜ i + ∂y2 u˜i + i ∂y2 φmi ∂x φ˜ i = 0, (35) 2
2
2 2 −∂y p˜i + Ci ∂x v˜i + ∂y v˜i + i ∂x φ˜ i + ∂y φ˜ i ∂y φmi + ∂y2 φmi ∂y φ˜ i = 0, ∂x2 φ˜i + ∂y2 φ˜ i = κi2 φ˜ i .
(36) (37)
(a) 1
q
0.8 0.6 0.4 0.2 0
−3
10
(b)
−2
−1
10
κ
10
0
1
10
2
10
3
10
10
0.5
0.4
τ
yy
m
0.45
0.35 0.3 −3
10
(c) 0.5
−2
−1
10
1
κ
10
2
10
10
3
10
(a) (b) (c)
0.25 pm and τe
oil with an approximate conductivity of ∼10−12 S/m has a Debye layer of ∼10−6 m. The value κ = 103 corresponds to a 4-mm-thick film which will show a conductor-like behavior because the κ −1 h0 , while κ = 10−2 corresponds to an approximately O(10)-nm film. For the same fluids (silicone oil), this system would now behave as a perfect dielectric. Thus for a given fluid with a κ ∗ , the relative thickness of the film determines the extent of bulk charges. The gradual variation in the mean potential distribution φm2 with changing inverse Debye length κ in the lower fluid is shown in the Fig. 3. The infinitesimal exponential potential drops at the surfaces for κ = 102 and increases with decreasing κ until they merge (dash-dot line κ = 5 in Fig. 3) and become linear at κ 1. Initially behaving like two PD fluids, a PD-LD interface actually behaves like a PD-PC (under dc fields and at long times) interface as the potential drop in the lower fluid is balanced by the charge separation across it and the effective field becomes zero. The high κ limit of the EK theory exactly matches the results from the PD-PC theory. The free charge at the interface between the two fluids, given 0 by −1/2 2 κ 2 (φm2 − φb ), can be used to confirm the limiting behavior of the EK model. In Fig. 4(a), the free charge is plotted as a function of κ denoted by the solid line. At very low κ, the free charge is almost zero. But as κ increases, the free charge in the Debye layer increases until, at high-enough κ, it matches the accumulated interfacial charge calculated using a perfect conductor model [denoted by the dashed line in Fig. 4(a)]. In Fig. 4(b), the Maxwell’s stress in the top fluid in the base state is plotted as a function of κ. In the top fluid, the Maxwell’s stress shows a similar trend as the free charge at the interface. It asymptotes to the value given by the PD-PD theory at low κ and that of a PD-PC theory at high κ. But in the lower fluid, the presence of the Debye layer leads to a different scenario [see Fig. 4(c)]. The base-state pressure given by Eq. (26) balances the Maxwell’s stress such that the total stress [Eq. (13)] follows a similar trend as that of the Maxwell’s stress in the top fluid. That is, the total stress from the EK theory [curve (b) in Fig. 4(c)] asymptotes to the limiting cases of PD-PD and PD-PC theories at low and high κ, respectively.
(d) (e) (f)
0
(g)
−0.25 −0.5 −4 10
(h) (i) −3
10
−2
10
−1
10
0
κ 10
1
10
2
10
3
10
FIG. 4. (a) A comparison of charge (q) accumulated at the interface obtained using the PD-PC model (indicated by the dashed line) and the free charge present in the Debye layer near the interface, in the lower fluid (ρ) obtained using the electrokinetic model (indicated by the solid line) varying as a function of the inverse Debye length κ. (b) Variation of Maxwell’s stress in the top fluid (τyy1 ), at the interface (y = 0), in the base state, as a function of the inverse Debye length κ. The dotted line indicates results from the PD-PD theory, the dashed line indicates results from the PD-PC theory and the solid lined curve indicates results from the electrokinetic model. (c) Variation of mean pressure [pm2 , curves (d), (e), and (f)], the electric stress [τe curves (g), (h), and (i)], and the total stress [τyy2 , curves (a), (b), and (c)] in the lower fluid, in the base state, as a function of the inverse Debye length κ. The dotted lines indicate results from the PD-PD theory, the dashed lines indicate results from the PD-PC theory, and the solid lined curve indicates results from the electrokinetic model. The values of parameters used are 1 = 1, 2 = 3, μr = 0, and β = 1.
The perturbation variables are now expressed in the normal ˜ is the perturbation variable mode form. Thus, if m ˜ ˆ m(x,y,t) = m(y,t)exp(ikx)dk, (38) ˆ where m(y,t) is the amplitude of the perturbation and k is the wave number of the perturbation. Substituting this form of the perturbation term in the Eqs. (34)–(37) and simplifying
032409-6
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
Eq. (37) can be written as 2 2
∂y − k − κi2 φˆi = 0.
PHYSICAL REVIEW E 89, 032409 (2014)
φˆ 1 = (e−ky − ek(y−2β) )c5 , (39)
Equations (39) and (40) are the governing equations of the system. The solutions to these are sought with appropriate boundary conditions in terms of perturbation variables. Expressing the boundary conditions listed in Sec. (II C) in terms of perturbations in the normal mode form, the no slip boundary condition at the interface becomes vˆ1 = vˆ2 ,
(41)
uˆ 1 = uˆ 2 .
(42)
The normal stress balance and the tangential stress balance upon substituting the stress definitions from Eqs. (13) and (14) become [[−pˆi + 2Ci ∂y vˆi + i ∂y φmi ∂y φˆ i ]] = −Bk 2 hˆ 2 ˆ y φmi = 0, (43) Ci (ik vˆi + ∂y uˆ i ) + i ik φˆ i ∂y φmi + i ik h∂ where B is a nondimensional number given by B = continuity of potential is given by
γ h0 . 0 φ02
ˆ y φmi ]] = 0, [[φˆ i + h∂
The (44)
where the second term indicates contribution because of the gradient in the mean potential. The normal component of the electric field, balanced across the interface, is written as [[−i ∂y φˆ i ]] = 0.
(45)
The dynamics of the position of the interface and the interfacial charge is given by the kinematic condition and the charge conservation equation, respectively. The kinematic condition becomes ∂t hˆ = vˆ2 .
(46)
A. Dispersion equation
For the case of time-invariant fields, the time dependence of the perturbation amplitude is assumed to be of the form est , where s is the growth rate of the perturbation. The sign of s decides the stability of the system. A negative s indicates a stable system, whereas the system is rendered unstable if s is positive. The second- and fourth-order ordinary differential equations (ODEs) in terms of φˆ i [Eq. (39)] and vˆi [Eq. (40)] are solved using the following boundary conditions: y → β or −1 vˆi = 0, ∂y vˆi = 0, and φˆ i = 0, respectively. The expressions obtained for vˆi and φˆ i are − (y + 2kyβ − 2kβ )c2 )], + [(−1 + e2k(1+y) )y + 2k(1 + y)c4 ],
(50)
V. RESULTS AND DISCUSSION
The objective of the present work is to overcome the limitation of electroneutrality of the bulk fluid imposed in the EHD models by accounting for the thermal motion of ions in the bulk. The electrokinetic model (EK) where the electric potential is modeled using the Debye-Huckel theory is used to overcome this limitation. Thereby, we assume that the dynamics of the ions is quite fast and the ion distribution is Boltzmann and the voltage or potential is low. The results from the EK theory are compared to those from the EHD model. The dynamics of the system when subjected to infinitesimal random perturbations can be gauged by determining the growth rate or the growth exponent s. The dispersion relation obtained from the linear stability analysis discussed in Sec. IV A is solved to obtain an expression for s the growth rate of the perturbation. This expression is plotted for different values of κ in Fig. 5. It can be seen from the figure that the growth rate of the perturbation depends on its wave number and the inverse Debye layer thickness κ and exhibits a maxima (smax ) for a particular value of wave number kmax which is the signature of the instability. An increase in κ results in 0.1
κ=0.01 κ=1
(47)
(48)
)c6 ,
ˆ where C is the column matrix of the constants c1 –c6 and h. This system of linear equations can have a nontrivial solution only when the coefficient matrix is singular. Therefore, equating the determinant of matrix A to zero we get a characteristic equation in terms of s and k and other system parameters like B, β, μr , 1 , 2 , σ1 , and σ2 . This equation is then solved using MATHEMATICA 7.0 for s, the growth rate of the perturbation.
0.05
κ=10 κ=100
0 −0.05 −0.1 0
vˆ2 = e−k(2+y) [(−1 + e2k(1+y) − 2k(1 + y)]c3
(−1 + e
(49)
AC = 0,
vˆ1 = e−ky [e2ky (c1 + yc2 ) + e2kβ ((−1 − 2ky + 2kβ)c1 2
k 2 +κ 2
√ 2(1+y) k 2 +κ 2
uˆ i and pˆ i are obtained in terms of vˆi from Eqs. (34), (35), and (36). These expressions for vˆi , uˆ i , pˆ i , and φˆ i , in terms of the constants c1 –c6 are substituted in the eight boundary conditions to get a system of homogeneous linear equations ˆ These equations in terms of the constants and the unknown h. can be represented by the matrix equation
s (growth rate)
Similarly, upon eliminating pˆ from Eqs. (35) and (36), we get a fourth-order differential equation in v which is given by 2
2 ∂y − k 2 vˆi = 0. (40)
φˆ 2 = −e−y
√
0.2
0.4 0.6 k (wavenumber)
0.8
1
FIG. 5. Variation of the growth rate s as a function of the wave number k at different values of κ. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
032409-7
PRIYA GAMBHIRE AND ROCHISH THAOKAR
PHYSICAL REVIEW E 89, 032409 (2014) TABLE I. Debye layer thicknesses for a given conductivity of the fluid.
0.7
k
max
0.6
S. no
0.5 0.4 0.3 −4 10
−3
10
−2
10
−1
10
κ
1
10
2
10
3
10
4
10
S. no
FIG. 6. Variation of kmax as a function of the inverse Debye length κ. The dotted line indicates results obtained from the PD-PD theory, the dashed line indicated those from the PD-PC, while the solid line indicates results from the electrokinetic theory. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
an increase in the growth rate due to a dominance of free charge. In Figs. 6 and 7, kmax and smax have been plotted as a function of κ. A similar trend as seen with the base state variables is seen with the instability parameters. The values of kmax and smax plateau to limiting values of PD-PD and PD-PC theories at low and high values of κ, respectively. For the case of an air-water system, where the conductivity of pure deionized water is taken to be 5.5 × 10−6 Sm−1 , the Debye length is obtained using the relationship σ = 0 Dκ 2 , where the diffusion coefficient D = 2.8 × 10−9 m2 s−1 is λD = κ −1∗ = ∼9 × 10−7 m. Typical electrode spacing in EHD-based soft lithography processes vary from 100 to 1000 nm [5,29,30] with a few groups reporting use of 80 and 100 μm [31]. Under these conditions, the thickness of the Debye layer is almost of the order of the electrode spacing. For other fluids like silicone oil and polymer melts with very low conductivities ranging from 10−9 to 10−12 S/m, the Debye length will further increase, although a low diffusion coefficient characteristic of these fluids ensures that it ranges from the submicron range to being of the order of electrode spacing. Therefore, an EK model is very necessary to describe the dynamics of the diffuse layers in these fluid-fluid systems.
0.06 0.05
s
max
0.04 0.03 0.02 0.01 0 −4 10
−3
10
−2
10
−1
10
κ
1
10
2
10
3
10
1 2 3 4 5
1 2 3 4 5
Conductivity (H2 O) S/m
Debye length (H2 O) m
1.0225×10−2 1.0225×10−3 1.0225×10−4 1.0225×10−5 1.0225×10−6
8.69293 × 10−9 2.748 × 10−8 8.69293 × 10−8 2.74894 × 10−7 8.69293 × 10−7
Conductivity (PDMS) S/m
Debye length (PDMS) m
−8
1.0225×10 1.0225×10−9 1.0225×10−10 1.0225×10−11 1.0225×10−12
4.364 × 10−8 1.374 × 10−7 4.364 × 10−7 1.374 × 10−6 4.364 × 10−6
VI. CONCLUSIONS
As one goes to lower thicknesses, the interfacial instability in thin fluid films under an electric field can be best described by an electrokinetic model which takes into account the dynamics of diffuse layers formed at the surfaces in the system. It is shown that the wavelength of the instability depends significantly on the Debye length, which is comparable to the fluid thickness for the low-conductivity fluids used in these studies. The PD-PD and PD-PC models are only limiting cases when the Debye length is very large and very small, respectively, when compared to film thickness. Thus it can be summarized that, in a perfect PD-PD system, free charges are absent and, hence, the electric body force term in the momentum balance equation is negligible. Even at the interface, free charges are absent, hence, the electric displacement continuity equation is used as a boundary condition. For a PD-PC system, with a large conductivity of the PC fluid, under a dc field the charges present in the conducting fluid move instantaneously to the surface of the fluid and reorganize into small Debye layers whose thickness is of the order of nanometers. The field within the fluid is nullified and the surface becomes an isopotential. The bulk of the fluid therefore can be considered as electroneutral. In case of leaky dielectrics, for conductivities of the order of 10−11 S/m, the Debye layer thickness is of the order of 1 to 5 μm, especially for viscous oils. When the thickness of the fluid film is much larger than a micron, the system is essentially a charged Debye layer separated by a core layer which is electrically neutral. An estimate of Debye layer thicknesses is given in Table I for typical aqueous and nonaqueous systems. Thus one can determine when such a system would behave as a perfect or a leaky dielectric (conductor) if the confinement is known.
4
10
FIG. 7. Variation of smax as a function of the inverse Debye length κ. The dotted line indicates results obtained from the PD-PD theory, the dashed line indicates those from the PD-PC, while the solid line indicates results from the electrokinetic theory. The values of parameters used are 1 = 1, 2 = 3, μr = 0, β = 1, and B = 1.
ACKNOWLEDGMENTS
The authors acknowledge the financial assistance provided by the Department of Science and Technology, India. The authors also acknowledge useful comments provided by V. Kumaran (IISc, Bangalore) and the reviewer.
032409-8
ELECTROKINETIC MODEL FOR ELECTRIC-FIELD- . . .
PHYSICAL REVIEW E 89, 032409 (2014)
APPENDIX: RELATIONSHIP BETWEEN THE INVERSE DEBYE LENGTH κ AND CONDUCTIVITY σ2
Using the DH approximation, −e[φ(y) − φb ] σ = e m+ n0 1 − kB T e[φ(y) − φb ] , + m− n0 1 + kB T
To be able to compare the electrokinetic model to the electrohydrodynamic model a relationship between the concentration of ions in the liquid and the conductivity of the liquid has to be determined. The electric conductivity (σ ) is related to the cationic and anionic molar concentrations c± by the following relation [22]: ˆ + c+ + m ˆ − c− ), σ = F (m 2
Assuming m+ = m− = m, −e[φ(y) − φb ] +1 = emn0 1 − kB T e[φ(y) − φb ] + , kB T
(A1)
where F is the Faraday’s constant given by F = eNa , where ˆ± e is the elementary charge and Na is Avagadro’s number. m are the mobilities of ions (units: mol N−1 ms−1 ) and c± are the molar concentrations, respectively, ˆ + c+ + m ˆ − c− ), σ = e2 Na2 (m ˆ + n+ + m ˆ − n− ), = e Na (m 2
= e(m+ n+ + m− n− ),
(A2)
where m± are the ionic mobilities and n± are the number concentration of charges, respectively. Substituting the expression for charge concentration from the equation below,
= 2emn0 .
2e2 n0 0 kB T and the Einstein’s relation D = mkB T /e, where D is the diffusivity (m2 s−1 ),
Using the expression for inverse Debye length κ 2 =
σ = D0 κ 2 .
∓e[φ(y)−φb ] kB T
n± = n0 e , −e[φ(y)−φb ] e[φ(y)−φb ] ⇒ σ = e μ+ n0 e kB T + μ− n0 e kB T .
(A4)
(A5)
(A3)
Therefore, the conductivity σ depends on the diffusivity, the dielectric constant, and the square of the Debye length.
[1] D. A. Saville, Annu. Rev. Fluid Mech. 29, 27 (1997). [2] S. Torza, R. G. Cox, and S. G. Mason, Philos. T. R. Soc. A 269, 295 (1971). [3] O. Vizika and D. A. Saville, J. Fluid Mech. 239, 1 (1992). [4] G. I. Taylor, Proc. R. Soc. London A 313, 453 (1969). [5] E. Schaffer, T. Thurn-Albrecht, T. P. Russell, and U. Steiner, Nature 403, 874 (2000). [6] G. I. Taylor and A. D. McEwan, J. Fluid Mech. 22, 1 (1965). [7] E. B. Devitt and J. R. Melcher, Phys. Fluids 8, 1193 (1965). [8] D. H. Michael, Proc. Camb. Philos. Soc. 64, 1203 (1968). [9] J. R. Melcher and C. V. Smith, Phys. Fluids 12, 778 (1969). [10] E. Schaffer, T. T. Albrecht, T. P. Russell, and U. Steiner, Europhys. Lett. 53, 518 (2001). [11] Z. Lin, T. Kerle, S. M. Baker, D. A. Hoagland, E. Schaffer, U. Steiner, and T. P. Russell, J. Chem. Phys. 114, 2377 (2001). [12] P. Gambhire and R. M. Thaokar, Phys. Fluids 22, 064103 (2010). [13] P. Gambhire and R. M. Thaokar, Eur. Phys. J. E 34, 84 (2011). [14] P. Gambhire and R. M. Thaokar, Phys. Rev. E 86, 036301 (2012). [15] W. B. Russel, D. A. Saville, and W. R. Scholwalter, Colloidal Dispersions (Cambridge University Press, Cambridge, 1989). [16] V. Kumaran, Phys. Rev. Lett. 85, 4996 (2000). [17] V. Kumaran, Phys. Rev. E 64, 011911 (2001). [18] D. Lacoste, G. I. Menon, M. Z. Bazant, and J. F. Joanny, Eur. Phys. J. E 28, 243 (2009).
[19] F. Ziebert, M. Z. Bazant, and D. Lacoste, Phys. Rev. E 81, 031912 (2010). [20] W. Qu and D. Li, J. Colloid Interface Sci. 224, 397 (2000). [21] C. K. Hua, D. W. Lee, and I. S. Kang, Colloid. Surface. A 372, 86 (2010). [22] E. K. Zholkovskij, J. A. Masliyah, and J. Czarnecki, J. Fluid Mech. 472, 1 (2002). [23] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevE.89.032409 for the derivation of the Debye-Huckel equation and calculation of ion concentration. [24] H. N. W. Lekkerkerker, Physica A 167, 384 (1990). [25] R. E. Goldstein, A. I. Pesci, and V. Romero-Rochin, Phys. Rev. A 41, 5504 (1990). [26] P. M. Vlahovska, K. D. Danov, A. Mehreteab, and G. Broze, J. Colloid Interface Sci. 192, 194 (1997). [27] A. Ramos, Electrokinetics and Electrohydrodynamics in Microsystems (Springer, New York, 2011). [28] M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, New J. Phys. 11, 075016 (2009). [29] M. D. Morariu, Z. Lin, T. P. Russell, and U. Steiner, Nature 2, 48 (2002). [30] P. Goldberg-Oppenheimer and U. Steiner, Small 6, 1248 (2010). [31] N. Arun, A. Sharma, P. S. G. Pattader, I. Banerjee, H. M. Dixit, and K. S. Narayan, Phys. Rev. Lett. 102, 254502 (2009).
032409-9
