Langmuir 2009, 25, 147-152
147
On the Connection between Dielectric Breakdown Strength, Trapping of Charge, and Contact Angle Saturation in Electrowetting Antonis I. Drygiannakis, Athanasios G. Papathanasiou,* and Andreas G. Boudouvis School of Chemical Engineering, National Technical UniVersity of Athens, Athens GR-15780, Greece ReceiVed August 6, 2008. ReVised Manuscript ReceiVed October 2, 2008 Electrowetting on dielectric (EWOD) is simulated by solving the equations of capillary electrohydrostatics, by the Galerkin/finite element method. Aiming to provide reliable predictions of the voltage dependence of the apparent contact angle, close to or beyond the saturation limit, special attention is given in the treatment of the dielectric properties of the solid dielectric where the liquid sits. It is proposed that in regions where the electric field strength locally exceeds the material breakdown strength, the dielectric locally switches to a conductor. Without using any fitting parameter, the implementation of the proposed phenomenological idea realized a surprising matching of published experimental data concerning materials ranging from SiO2 to Parylene N and Teflon. Charge trapping is naturally connected to the field-induced transition, and its distribution as well as its dependence on the applied voltage is calculated.
1. Introduction Tunable charge distribution on a solid/liquid interface through the application of an external voltage realizes an attractive way of dynamically modifying the wetting properties of solids. Electrostatically assisted wetting of dielectric solids by conducting liquids (electrowetting on dielectric, EWOD) is recently applicable to many microfluidic devices implementing microdroplets manipulation.1 Constructed using microlithography techniques and operating without moving mechanical parts, electrowetting devices combine low cost, small size, and robustness in terms of energy consumption and reliability. Promising electrowetting applications, however, are limited not only from the well-known saturation phenomenon but also from the limited applicability of the Lippmann’s equation at low applied voltages.2,3 Lippmann’s equation gives the dependence of the apparent contact angle, θv, on the voltage, V, applied between a conducting liquid droplet and a dielectric-coated planar electrode where the liquid sits
cos θv ) cos θy + ε0εrV2 ⁄ (2δγ)
(1)
where θy is the Young’s contact angle and γ the liquid surface tension. δ is the thickness of the dielectric with constant εr, and ε0 is the permittivity of vacuum. It is a statement of the total interfacial energy minimization, namely, the sum of the liquid/ vapor (surrounding the drop), the solid/vapor, and the solid/ liquid interfacial energy augmented with the electrostatic energy stored in the solid/liquid interface. Neglecting the charge distribution at the liquid/vapor interface and assuming uniform electric field at the solid/liquid interface due to the ideal parallel plate capacitor approximation, Lippmann’s equation predictions are limited to low applied voltages. The desirable complete wetting, i.e., θv ) 0, although predicted by the Lippmann equation, is experimentally unreachable since wetting enhancement stops at a contact angle limit due to the saturation phenomenon. * Corresponding author: E-mail:
[email protected]. Web page: http://www.chemeng.ntua.gr/people/pathan. Phone: +30210 772 3290. Fax: +30 210 772 3155. (1) Shamai, R.; Andelman, D.; Berge, B.; Hayes, R. Soft Matter 2008, 4(1), 38–45. (2) Quilliet, C.; Berge, B. Curr. Opin. Colloid Interface Sci. 2001, 6(1), 34– 39. (3) Vallet, M.; Vallade, M.; Berge, B. Eur. Phys. J. B 1999, 11(4), 583–591.
The origins of the contact angle saturation are not yet fully understood, although many theories have been proposed to explain various related observations.4 The high field strength developed at the wedgelike geometry of the liquid at the TPL appears to be a key issue since it is related to charge trapping in the dielectric5 and charge leakage through the dielectric6-8 or through the surrounding air due to ionization.3 According to a completely different approach, saturation is a thermodynamic limit reached when an effective solid/liquid interfacial tension accounting for the electric field vanishes.9 Despite the plethora of the proposed theories, available theoretical predictions are still far from accurately simulating electrowetting close to the onset or beyond the saturation limit. In many cases theoretical predictions are still drawn from variations of Lippmann’s simplified equation augmented with fitting parameters introduced for matching experimental data. Here, the governing equations of the capillary electrohydrostatics are applied to the classic electrowetting system of the axisymmetric sessile droplet. Without making any simplification concerning the droplet shape or the field distribution especially close to the three phase contact line (TPL), a nonlinear and free boundary problem is formed which is discretized by the Galerkin/ finite element method. The nonlinear algebraic equation set that arises is solved by Newton iteration. Field distribution and droplet shape are computed simultaneously, and their dependence on physical parameters, like the applied voltage or the solid dielectric properties, is studied by parameter continuation. Special care is given in the computation of the field distribution and its effect on the solid dielectric properties in the vicinity of the TPL. By locally treating the dielectric as a conductor in regions where the field strength exceeds its breakdown strength, electrowetting is successfully simulated even beyond the saturation onset as the matching of published experimental data showed. The proposed transition is connected with the charge trapping (4) Mugele, F.; Baret, J. C. J. Phys.: Condens. Matter 2005, 17(28), R705. (5) Verheijen, H. J. J.; Prins, M. W. J. Langmuir 1999, 15(20), 6616–6620. (6) Papathanasiou, A. G.; Boudouvis, A. G. Appl. Phys. Lett. 2005, 86(16), 164102. (7) Papathanasiou, A. G.; Papaioannou, A. T.; Boudouvis, A. G. J. Appl. Phys. 2008, 103(3), 034901. (8) Seyrat, E.; Hayes, R. A. J. Appl. Phys. 2001, 90(3), 1383–1386. (9) Quinn, A.; Sedev, R.; Ralston, J. J. Phys. Chem. B 2005, 109(13), 6268– 6275.
10.1021/la802551j CCC: $40.75 2009 American Chemical Society Published on Web 12/03/2008
148 Langmuir, Vol. 25, No. 1, 2009
Drygiannakis et al.
air, and K is a reference pressure (constant along the free surface). h is the vertical distance of the free surface from the reference level, and C is its local mean curvature
[
]
2f 2 sin2 θ + fθ2 sin2 θ ffθ sin2 θ ∂ 1 C) 2 f sin θ √sin2 θ(f 2 + fθ2) ∂θ √sin2 θ(f 2 + fθ2) (4) fθ ≡ ∂f ⁄ ∂θ The magnitude, E, of the electric field strength in eq 2 is calculated at the free surface of the liquid. The spatial field distribution is governed by the equations of electrostatics, which are, as they apply here, the Gauss law for the electricity
∇ · (εr ∇ u) ) 0 Figure 1. Electrowetting in an axisymmetric sessile drop.
initially proposed by Verheijen and Prins.5 Trapped charge distribution is naturally calculated resolving many issues regarding its magnitude and localization, raised by Verheijen and Prins.
2. Electrohydrostatics of Sessile Drops A typical system for studying electrowetting is an axisymmetric liquid drop sitting on a horizontal surface of a dielectric solid, which in turn coats a metal electrode as shown in Figure 1. The drop is conducting and is surrounded by a dielectric fluid, in this case air. When a voltage is applied between the metal electrode and an electrode submerged in the liquid, electric charge accumulates in the solid/liquid interface causing a decrease in the corresponding interfacial energy, thus giving rise to an enhancement of wetting of the solid by the liquid. The geometry of the system suggests the implementation of a composite coordinate system. Cartesian coordinates (resulting from the projection of cylindrical coordinates in a plane parallel to the axis of symmetry) (x,z) are used in the region occupied by the solid dielectric (region 1 in Figure 1), whereas polar coordinates (resulting from the projection of spherical coordinates in a plane parallel to the axis of symmetry) (r,θ) are used above the dielectric (region 2 in Figure 1). The unknown drop free surface shape is parametrized in terms of θ, i.e., r ) f(θ). Equilibrium deformations are sought for; thus force balance at the free surface is required. The balance between gravity, surface tension, and electrostatic stresses is stated by the Young-Laplace equation of capillary hydrostatics augmented with the electric stress term.10,11
-Ngh +
NeE2 +C)K 2
(2)
The equation above as well as all the following equations are dimensionless. Ng and Ne are the gravitational and the electric Bond numbers respectively
g∆FR02 Ng ≡ γ
ε0V2 Ne ≡ γR0
(3)
where R0 is the characteristic length, which is the radius of a sphere of volume equal to the drop volume. g is the acceleration of gravity, ∆F is the density difference between the liquid and (10) Kang, K. H. Langmuir 2002, 18(26), 10318–10322. (11) Basaran, O. A.; Scriven, L. E. Phys. Fluids A 1989, 1(5), 799–809.
(5)
where u is an electric potential (E ≡ ∇u) and εr is the dielectric constant (εr ) 1 for air, εr ) εd for the solid dielectric). Here, we do not account for the dielectric properties of the hydrophobic top coating since its thickness is usually of the order of 10-2 µm, a scale much smaller than the dielectric thickness. The liquid is considered as perfect conductor and therefore the potential, u, is constant in the entire drop. Due to the incompressibility of the liquid, the volume remains constant at any drop deformation, giving rise to the following integral constraint for the volume conservation
∫0π ⁄ 2 f 3 sin θ dθ ) 2
(6)
The differential equations (2) and (5) together with the integral equation (6) are solved simultaneously for the electrostatic potential, u, inside the dielectric and the surrounding air, the free surface shape, f, and the reference pressure K, accounting for the following boundary conditions:
u2 ) 1 u1 ) 1
at
at
z)0
r ) f(θ)
(7a)
and
(7b)
r e f(π ⁄ 2)
u1 ) 0
at
z ) -δ
(7c)
u1 ) u2,
εd
∂u1 ∂u2 ) ∂z ∂z
(8)
at
z)0
and
f(π ⁄ 2) < r < Rmax
∂u2 ) 0 at r ) Rmax ∂r ∂u1 ) 0 at x ) Rmax ∂x ∂u2 ∂f ) 0, ) 0 at θ ) 0 ∂θ ∂θ ∂u1 ) 0 at x ) 0 ∂x n · ez ) cos θy at θ ) π⁄2
(9a) (9b) (10a) (10b) (11)
ez ) cos θer - sin θeθ where u1 and u2 are the potential inside the regions 1 and 2, respectively, ez, er, eθ are the unit vectors in the directions z, r, and θ, respectively, and n is the unit normal to the liquid free surface
Breakdown, Charge Trapping, Saturation in Electrowetting
n)
fer - fθeθ
√(f 2 + fθ2)
Langmuir, Vol. 25, No. 1, 2009 149
(12)
The applied voltage difference between the electrodes is imposed through the eqs 7a-7c. Equation 8 states the continuity of the electric potential and of the normal component of the dielectric displacement across the dielectric/air interface. Equations 9a and 9b are applied at the boundary of the domain, which is located far enough from the liquid drop. The axial symmetry of the drop is imposed through the eqs 10a and 10b. The wetting properties of the liquid are accounted in eq 11 through Young’s contact angle, θy. θy depends on the liquid/ dielectric, liquid/air, and dielectric/air interfacial tensions and is taken as constant and independent of the electric field. Young’s angle used here is the one of the hydrophobic top coats (usually Teflon AF) widely used in electrowetting systems. The resulting mathematical problem is nonlinear and free boundary since the unknown free surface shape affects the electric field distribution and vice versa.
3. Computational Analysis The equations governing the field distribution, coupled with the force balance equation and the volume constraint, are solved simultaneously by a compact numerical scheme based on the Galerkin/finite element method.12 To facilitate reading, the discretization of the governing equations is presented in the Appendix. The exploration of the solution space (parametric analysis), i.e., the study of the dependence of the equilibrium states on the problem parameters, like the applied voltage, V, is performed by parameter continuation. Taking into account that the liquid surface forms a wedge at the three-phase contact line (TPL), and therefore the charge density is high, the required mesh (see a sample of the mesh in Figure 2) refinement significantly increases the computational power demands. A typical size of the computational problem is of the order of 60000 unknowns. The accuracy of the results is tested against discretization refinenement. Moreover, alternative methods are employed for the efficient computation of the field at the TPL, such as singular elements, boundary elements, and an efficient combination thereof; preliminary results are reported in ref 13. As it shown in section 2, the effect of the electric field is not localized at the TPL as the simplified Lippmann’s equation assumes.14 Here, the distribution of the electric stress tensor along the entire droplet free surface is accounted for. The contact angle appearing in the results that follow is the apparent contact angle, θv, calculated by a third-order polynomial fitting near the
Figure 2. A sample computational mesh. Local magnification shows mesh details close to the TPL.
Figure 3. The dependence of θv on the applied voltage for 1 µm SiO2. (εr ) 3.8 (SiO2), δ ) 1 µm, θy ) 114°, γ ) 0.072 N/m, Ebd ) 10 MV/cm; top coat, 20 nm AF Teflon.)
TPL; it is reminded that the Young contact angle, θy, is independent of the electric field.14
4. Results and Discussion According to our recent findings, the field distribution in the vicinity of the three-phase contact line determines the onset of the contact angle saturation.6 It is reported that electrowetting is limited at a critical applied voltage where the field strength, averaged in small quadrilateral elements of finite size 10 nm, locally exceeds the breakdown strength of the dielectric material where the liquid sits. Although the onset of the saturation and the corresponding contact angle are quite well predicted for many electrowetting systems, the proposed methodology failed to capture the contact angle dependence beyond the saturation. Moreover the selection of the element size for the field strength averaging is important since the prediction of the saturation onset is sensitive on this parameter. Keeping the dielectric properties of the solid dielectric independent of the field strength distribution, the solution of the electrohydrostatic problem did not show any asymptotic stabilization of the contact angle at high voltages. Buehrle et al. also report a similar result.15 The issue of implementing a constitutive equation that accounts for a possible dependence of the dielectric constant on the field strength distribution was raised in order to treat the deficiency of the electrohydrostatic solution.6 Here we adopt a simple phenomenological technique to deal with the dielectric behavior of the solid dielectric at field strengths higher than its breakdown strength. We propose that in regions where locally the breakdown strength is exceeded, the dielectric flips to a perfect conductor. Treating the corresponding region as a conductor implies that the potential in the entire region is set to a constant valuesthe potential of the liquid. As the applied voltage increases up to a critical voltage where the breakdown strength of the dielectric is surpassed, regions belonging to the solid dielectric and adjacent to the TPL gradually flip from the dielectric to the conductor state. Thus a conductivity front is created evolving from the TPL toward the solid dielectric phase. This transition limits electric field strength below the breakdown strength of the solid material, induces charge concentration in (12) Strang, G.; Fix, G. J. An Analysis of the Finite Element Method; PrenticeHall: Englewood Cliffs, NJ, 1973. (13) Pashos, G.; Papathanasiou, A. G.; Boudouvis, A. G. Treatment of electric field singularities at solid-liquid-air junctions in electrowetting-on-dielectric computations. In 8th World Congress on Computational Mechanics (WCCM8) & 5th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2008), Lido Island in Venice, Italy, 2008. (14) Mugele, F.; Buehrle, J. J. Phys.: Condens. Matter 2007, (37), 375112. (15) Buehrle, J.; Herminghaus, S.; Mugele, F. Phys. ReV. Lett. 2003, 91(8), 086101.
150 Langmuir, Vol. 25, No. 1, 2009
Drygiannakis et al.
Figure 4. Potential (a) and field strength (b) distribution in the vicinity of the TPL at selected points (shown in circles) of figure 3.
the dielectric region that changed state, and consequently limits wetting enhancement. This concentrated charge is directly associated with what is in the literature called trapped charge5 and its origin is considered to be the high electric field developed at the area of the TPL in connection with local dielectric breakdown. The implementation of the idea through the finite element method is quite conveniently incorporated in the parameter continuation. At any step of the continuation, namely, the simultaneous computation of the equilibrium shape and field distribution for a certain applied voltage, the magnitude of the field strength is monitored in the center of mass of each of the elements belonging to the area occupied by solid dielectric (region 1 in Figure 1). If the field strength in any element exceeds the dielectric strength, Ebd, of the solid, the corresponding element is then treated as a conductor. This is implemented by explicitly replacing the potential at all of its nodes with the potential value of the adjacent conductor; it is obvious that since the highest field strength is localized at the TPL the elements that change state get the potential of the liquid which is constant and equal to u ) 1 (cf. eqs 7a and 7b). Then the equilibrium computation is run again to recompute/correct the shape and field distribution and check whether the field strength in any other element exceeds the breakdown strength of the solid. The continuation proceeds further by increasing the applied voltage and repeating the steps outlined. The proposed technique enabled the modeling of electrowetting even beyond the saturation limit. In Figure 3, the computed dependence of the apparent contact angle, θv, on the applied voltage, V, is shown for a sessile droplet sitting on a SiO2 dielectric layer (d ) 1 µm). It is clearly seen that the proposed computations not only predict the asymptotic lower limit of the contact angle but also successfully capture the deviation of the experimental results (which are shown in Figure 3) from the ideal Lippmann behavior. For this system, the computations showed that the maximum field strength, localized at the TPL, exceeds the SiO2 breakdown strength at V ) 64 V (see Figure 3). It is evident from the proposed idea that increased dielectric strength results in better electrowetting performance and vice versa (see also ref 7). Computational predictions showed limiting contact angles of 83° and 59° if the dielectric strength of SiO2 is chosen equal to 8 and 12 MV/cm, respectively (instead of 10 MV/ cm). Along selected points of the computational curve in Figure 3, the potential and the field strength distribution close to the TPL are shown in panels a and b of Figure 4, respectively. Figure 4a
Figure 5. Charge density distribution, |σ|, along the line AB shown in the inset. r measures from the axis of symmetry. (εr ) 3.8 (SiO2), δ ) 1 µm, θy ) 114°, γ ) 0.072 N/m, Ebd ) 10 MV/cm.)
shows a growing conductive region nucleating from the TPL (see the second snapshot from the left) end expanding along the solid/air interface as the applied voltage increases up to 80 V. Consequently the location of the highest field strength is gradually shifted away from the TPL (see Figure 4b). The field strength at the TPL is limited accordingly resulting to the limiting of electrowetting; the apparent contact angle (cf. Figure 3) attains the lowest limit of 75°. Any increase of the applied voltage does not induce increase of the maximum field at the TPL, which would decrease the contact angle, but simply causes further propagation of the conducting region. Notice that the highest field strength is lower than 10 MV/cm (the SiO2 breakdown strength) at any applied voltage, as the color scales in Figure 4b show. The boundary of the region that switched from dielectric to conductive is charged, and its charge density distribution along the solid/air interface is presented in Figure 5, for selected values of the applied voltage. The corresponding region, labeled with ABCD, is schematically shown in the inset of Figure 5. The computed charge density is now negligible at the TPL (point A); however it sharply increases close to the right boundary (the line BC) of the region. Such a distribution indicates a concentration of the trapped charge at some small distance away from the TPL. This distance, or equivalently the length, L, of the line AB,
Breakdown, Charge Trapping, Saturation in Electrowetting
Langmuir, Vol. 25, No. 1, 2009 151
Figure 6. Mean charge density and trapped charge density vs the applied voltage. Curve 1 is the mean charge density at solid/liquid interface. Curve 2 is the mean charge at the elements at which E > Ebd. (εr ) 3.8 (SiO2), δ ) 1 µm, θy ) 114°, γ ) 0.072 N/m, Ebd ) 10 MV/cm.)
increases by increasing the applied voltage as the curves in Figure 5 show. For the case studied, L (L ) R2 - R1, see the inset in Figure 5), varies between 0.01 and 0.13 (dimensionless value) for applied voltages between 68 and 80 V, respectively, indicating that it can reach the order of magnitude of d ) 0.357, the dielectric thickness. It is reminded that the locations of the TPL and of the point B depend on the applied voltage, since not only TPL moves as the droplet deforms but also the conductive region, ABCD, grows as the voltage increases. A measure of the charge trapped in the dielectric is its spatial mean value, calculated by averaging the surface charge over the boundary of the region that switched from dielectric to conductive. Figure 6 shows that the mean trapped charge significantly increases when the applied voltage exceeds the saturation onset, sharply attaining values almost half of the mean charge density of the solid/liquid interface. The predictions based on the proposed computational technique were tested against experimental measurements concerning different thicknesses of SiO2 as well as different non-silicon materials, widely used in electrowetting systems, like Parylene N and Teflon AF. Experiments with SiO2 were performed using a standard setup presented in ref 7 whereas data for Parylene N and Teflon AF (see Figure 7b) were obtained from published works.3,5 The agreement with the experimental measurements is remarkable. The computations not only capture the deviation of the experiments from the simplified Lipppman equation for the thin SiO2 (see Figure 7a) but also successfully predict the very low saturation contact angle (namely 28°) that the 50 µm thick Teflon AF sample attains. Notice that this particular low limit was considered to be thermodynamically unstable according to Quinn et al.9 a statement which is in contrast to our equilibrium computations. Since our technique is phenomenological and based on the continuum hypothesis, the detailed physics of the system cannot be captured in a scale that could elucidate the nature of the predicted trapped charge. However we believe that our work, in addition to the successful electrowetting simulation, offers a convenient base for a multiscale modeling framework, utilizing mesoscale or eventually molecular level treatment of the TPL region, coupled with the traditional finite element methodology. This way the saturation mechanism could be elucidated in a resolution that eventually could allow the design of new materials which inhibit saturation to higher voltages and lower contact angles.
Figure 7. The dependence of θv on the applied voltage for different dielectric materials and thicknesses. (a) for SiO2 samples (εr ) 3.8, θy ) 114°, γlv ) 0.072 N/m, Ebd ) 10 MV/cm; top coat, 20 nm AF Teflon) and (b) for Parylene N (εr ) 2.65, δ ) 10 µm, θy ) 119°, γlv ) 0.072 N/m, Ebd ) 4 MV/cm) and PTFE (εr ) 2.0, δ ) 50 µm, θy ) 102°, γlv ) 0.072 N/m, Ebd ) 2.2 MV/cm) samples. Published experimental measurements are also shown for Parylene N5 and PTFE.3
5. Summary and Conclusions By employing the Galerkin/finite element method, electrowetting is simulated by solving the equations of electrohydrostatics in the system of a sessile droplet. Special attention is given in the treatment of the solid dielectric in the vicinity of the TPL where high field strength is developed due to the liquid wedgelike geometry. The solid dielectric is locally treated as a conductor when the field strength exceeds its breakdown strength. Material breakdown strength is connected with dielectric charge trapping. Its spatial distribution and the dependence on the applied voltage are computed. The implementation of the simple proposed idea realized the successful simulation of the contact angle dependence on the applied voltage, even beyond the saturation limit. It is for the first time presented a computational technique where the predicted contact angle asymptotically stabilizes to a limiting value at high enough applied voltage and most importantly matching published experimental measurements by accounting only for physical properties. Limitations of the proposed technique are discussed and future analysis is suggested to illuminate molecular charge trapping or current leakage mechanisms and their effects on the geometry of the three-phase contact line. Acknowledgment. We are indebted to Dr. A. D. Tserepi at the Institute of Microelectronics of the NCSR “Demokritos”, for
152 Langmuir, Vol. 25, No. 1, 2009
Drygiannakis et al.
preparing and kindly providing the SiO2 dielectric samples. Detailed computations, by G. Pashos, of the electric field at the three-phase contact line were particularly enlightening. This work has been funded by the projects ΠENE∆ 2003 and ENTEP 2005. The projects are co-financed 80% of public expenditure through EC—European Social Fund, 20% of public expenditure through Ministry of Development—General Secretariat of Research and Technology and through private sector, under measure 8.3 of OPERATIONAL PROGRAM “COMPETITIVENESS” in the 3rd Community Support Program.
Appendix Galerkin/Finite Element Discretization of the Governing Equations The computational domain is an orthogonal parallelogram in region 1, whereas in region 2 there is a quadrant where the region occupied by the drop is subtracted, since the potential is constant therein. The unknowns are approximated in terms of finite series of the finite element basis functions
8, 9a, 9b, 10a, and 10b are directly accounted for in the line integrals of eqs 15a and 15b. Along the air/dielectric boundary, SSV, n points toward the air in region 1, whereas in region 2 the unit normal is -n. Therefore the net contribution of the line SSV to the surface integrals of eqs 15a and 15b reads
∫ (Φaiεrn · ∇ u1 - Φbjn · ∇ u2) dSSV
and vanishes due to the second interfacial condition of eq 8 and since Φai,Φbj coincide along SSV. The Dirichlet boundary conditions that prescribe the potential at parts of the boundaries are accounted for by replacing the corresponding residuals with equations setting the desired values of the potential. The Galerkin weighted residuals of the Young-Laplace equation (eq 2) are constructed by multiplying the equation by each of the basis functions, Fκ, and integrating over the liquid free surface
RYLk ≡
L
u(x, z) )
(13a)
1
in region 1
∑ ujΦbj(r, θ)
(
dSLV ) f sin θ√(f 2 + fθ2) dθ
(13b)
1
N
∑ fκF (θ) κ
(13c)
1
at the drop free surface. L and M are the number of the unknowns at the nodes created by the discretization of the regions 1 and 2, respectively (see Figure 1). N is the number of nodes at the free liquid free surface. The basis functions, Φai(y,z), Φbj(r,θ), and Fκ(θ) are linear polynomials in each independent variable. Their coefficients in eqs 13a, 13b and 13c are the unknown values of the electric potential at the nodes of the discretization and the radial coordinates of the nodes of the free surface, respectively. The Galerkin weighted residuals of eq 5 are built by multiplying the equation by each of the basis functions ΦRi, Φβj and integrating the product over the entire domain
Rip ≡
∫ ∫Φai ∇ · (εr ∇ u) dS1
(14a)
∫
π⁄2
0
∫
CFif 2 sin θ dθ )
π⁄2
0
[
(2f 2 sin2 θ + fθ2 sin2 θ)Fi
(14b)
S2
in region 2 By integrating by parts and after the application of the divergence theorem the residuals read:
I ∂S1
∇u
dλ1 - ∫∫ εr ∇ Φai · S1
∇ u dS1 (15a)
Rpj ) I Φbjn · ∇ u dλ2 - ∫∫ ∇ Φbj · ∇ u dS2 ∂S2
θ(f + fθ ) 2
2
ffθ sin2 θFθi
√sin
2
θ(f + fθ ) 2
2
P≡
ffθ sin2 θF i
√sin θ(f 2
2
+ fθ2)
S2
(15b)
The line integrals in eqs 15a and 15b are calculated at the boundaries of the regions 1 and 2 respectively. n is the unit normal to the boundary of the each of the regions, always pointing outward. The Neumann boundary conditions appearing in eqs
]
dθ - [P]π0 ⁄ 2
) F if(n · ez)|θ)π ⁄ 2
(19)
(20)
The boundary condition that sets the contact angle equal to the Young’s contact angle (eq 11) is accounted in eqs 19 and 20
P ) Fif(n · ez)|θ)π ⁄ 2 ) Fif cos θy
∫ ∫Φbj ∇ · (∇u) dS2
εrΦain ·
√sin
2
+
where
in region 1
Rpi )
(18)
The integration over the projection of the free surface facilitates the integration by parts of the curvature term and thus the reduction of the second-order derivatives in the first ones
S1
Rpj ≡
(17)
dSp)f 2 sin θ dθ
in region 2
f(θ) )
)
NeE2 F -Ngh + + C - K dSp 2 k
where dSp is the area element which is the projection of the free surface area element, dSLV is in a sphere of radius r ) f(θ). The area elements, dSp, and dSLV are respectively:
M
u(r, θ) )
∫ Sp
∑ uiΦai(x, z)
(16)
SSV
(21)
The volume constraint (eq 6) is already an integral equation and is used as is; therefore the corresponding residual is the following
Rv )
∫0π ⁄ 2 f 3sin θ dθ - 2
(22)
All the integrals at the Galerkin residuals are computed by a Gauss quadrature in a standard isoparametric element.12 By vanishing the residuals, a system of nonlinear algebraic equations arises, which is solved for the nodal unknowns of the potential, u1, u2, the location of the nodes of the free surface, f, and the reference pressure K. The set of the nonlinear algebraic equations is linearized and solved by the Newton’s method combined with a frontal solver of the sparse and structured linear system that arises. LA802551J
