ELECTROSTATIC INTERACTIONS B E T W E E N INTERFACIAL PARTICLES By Michael Peter Lyne B. Sc. (Chemistry Co-Op) University of Waterloo M. Sc. (Chemistry) University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF A P P L I E D SCIENCE in T H E FACULTY O F GRADUATE STUDIES CHEMICAL ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA Augustl989 © Michael Peter Lyne, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Chemical Engineering The University of British Columbia 6224 Agricultural Road Vancouver, Canada V6T 1W5 Date: Abstract The Finite Element Method (FEM) is employed, using Galerkin weighted residuals, to simulate the electrostatic interactions between interfacial particles in two separate models. The three phase problem of two identical, parallel, infinitely long, solid cylinders at a given separation and each embedded to the same extent in an oil/water interface is considered. It is revealed that due to the asymmetric charge distribution surrounding an interfacial particle, a force (normal to the interface) acts on the particle in the direction of the electrolyte phase. The electrostatic forces acting on the particles normal and parallel to the oil/water interface are calculated and it is found that under certain conditions they may be of comparable magnitude. The Debye-Hiickel form of the Derjaguin method is also applied to this model; this method showed good agreement with the low potential numerical results for the parallel force when the oil/water interface is assumed to be uncharged. Interactions between the cylinders on a constant potential interface are also considered and it is revealed that for some circumstances the magnitude of these interactions will be significantly different from the corresponding interactions on an uncharged interface. For practical situations, however, increases in the interaction force were not large enough to account for a 10 fcT change in the free energy associated 6 with the interaction. The more complex interaction between particles in a dense, two-dimensional structure on an oil/water interface is simulated using a cell model approach. The normal forces acting on the spheres and the electrostatic free energy change associated with forming the structure are calculated. The normal force per unit of particle perimeter in contact with the oil/water interface is shown to be comparable to that acting on cylinders while ii the free energy calculations are only accurate enough for order of magnitude estimations. Additionally, to test the effectiveness of the numerical approach, the farniliar case of two identical interacting spheres completely immersed in a single electrolyte is modeled for both constant potential and constant charge boundary conditions on the sphere surfaces. The implicit assumption inherent in the Derjaguin method when it is applied to constant charge interactions (i.e. zero potential gradient inside the particle) is examined and found to cause errors of > 10% in cases of low na and large particle dielectric constants. 111 Table of Contents Abstract ii List of Figures v List of Tables vi Acknowledgement vii 1 Introduction 1 2 Theory 4 2.1 Electric Double Layer Surrounding an Isolated Charged Solid 4 2.2 Electrostatic Interactions between Charged Surfaces in an Electrolyte Solution 11 2.3 Electrostatic Interactions between Particles at an Electrolyte Interface . . 21 2.4 Model Descriptions 23 2.4.1 Two Identical Spheres Fully Immersed in an Electrolyte 25 2.4.2 Two Identical Cylinders Embedded in an 0/W interface 27 2.4.3 Isolated Sphere Embedded in an 0/W Interface in a. Unit Cell . . 30 2.5 Finite Element Method 32 2.5.1 Introductory Remarks 32 2.5.2 Galerkin Method of Weighted Residuals 35 2.5.3 Choice of Basis Functions 36 2.5.4 Weak Formulation via Green's Theorem 37 iv 2.5.5 Application to Thesis Models 40 3 Methods 43 3.1 Finite Element Grid 43 3.2 Sensitivity Tests 47 3.2.1 Parallel Force Between Two Particles 48 3.2.2 Normal Force Acting on a Particle at an Interface 52 3.2.3 Free Energy 57 3.3 Convergence Tolerance and Relaxation Parameter 61 3.4 Orders of Integration Used 62 3.5 Solution of Equations 3.6 Software Overview 64 3.6.1 Data Input 64 3.6.2 Grid Construction 65 3.6.3 Setting up the Equations 66 3.6.4 Iterative Convergence 67 3.6.5 Force and Free Energy Calculations 68 .- 4 Results and Discussion 4.1 4.2 63 69 Two Identical Spheres Fully Immersed in an Electrolyte 70 4.1.1 Constant Potential Interactions 71 4.1.2 Constant Charge Interactions 76 4.1.3 Summary 83 Two Identical Cylinders Embedded in an Oil/Water Interface 4.2.1 84 Constant Surface Potential Cylinders Embedded in an Uncharged O/W Interface 84 v 4.2.2 Constant Surface Potential Cylinders Embedded in a Constant Potential 0/W Interface 4.2.3 4.2.4 4.3 99 Constant Surface Charge Cylinders Embedded in an Uncharged 0/W Interface 103 Summary 112 Close-Packed Two-Dimensional Swarm of Spheres at an 0/W Interface . 113 4.3.1 Normal Force Acting on a Constant Potential Sphere Embedded in an Uncharged 0/W Interface 4.3.2 4.3.3 5 6 114 Free Energy Calculations For a Constant Potential Sphere Embedded in an Uncharged 0/W Interface 119 Summary 121 Conclusions 123 5.1 Model I 123 5.2 Model II 124 5.3 Model III 125 Recommendations 126 6.1 Integrating the Interaction Force 127 6.2 Improving the Finite Element Grid 128 6.2.1 Algebraic Grid Generation 130 6.2.2 Elliptic Grid Generation 132 6.2.3 Adaptive Grids 133 6.2.4 PLTMG 134 6.3 Distortion of the Interface 135 Bibliography 137 vi Appendices 141 A Derivation of Force Equations for Three Phase System 141 B Quadratic Order Triangular Elements 145 B.l C Basis Functions and Their Derivatives . 145 B. 2 Constant Charge Boundary Condition 147 F E M Software 149 C. l 149 POWI.CYL : Model II C.2 ENERGY Subroutine 182 C.3 Typical Input File 184 D Supplemental Results 185 E Nomenclature 201 vii List of Figures 2.1 Derjaguin Method Applied to Equal, Parallel Cylinders 23 2.2 Model I: Two Identical Spheres in an Electrolyte 26 2.3 Model II: Two Identical Cylinders at a Planar 0/W Interface 28 2.4 Model III: Unit Cell Representation of a Sphere at an 0/W Interface . . 31 3.1 Typical Grids for the Three Models 45 3.2 Electrostatic Potential Distribution for One of Two Interacting Spheres in an Electrolyte 49 3.3 Effect of NEA and NED on the Parallel Force 51 3.4 Effect of ZOI and XI on the Parallel Force for Minimal Double Layer Overlap 52 3.5 Potential Distribution for One of Two Parallel, Interacting Cylinders on an Uncharged 0/W Interface 54 3.6 Back Side Normal Force as a Function of FCXI 56 3.7 Back Side Normal Force as a Function of NEA 57 3.8 Free Energy of a Sphere in an Electrolyte as a function of NEA and FCXI 59 3.9 Free Energy of a Sphere in an Electrolyte as a function of NWA 60 3.10 Free Energy of a Sphere in a Unit Cell as a function of ZOI 60 3.11 Converged Value of/* as Function of 8 61 4.1 Fraction of'Complete' to 'Incomplete' Force Solutions for Constant Charge Interactions between a Pair of Spheres in a Single Electrolyte vin 77 4.2 Comparison of 'Incomplete' Numerical Solutions and Derjaguin Approximations for the Interaction Force between a Pair of Spheres in a Single Electrolyte 4.3 82 Low Potential Derjaguin Approximation Compared with Numerical Predictions of the Force Acting between Identical, Parallel Cylinders Embedded in an Uncha rged 0/W Interface 4.4 Potential Distribution for One of Two Parallel Cylinders on an Uncharged 0/W Interface, 6=45° 4.5 89 Potential Distribution for One of Two Parallel Cylinders on an Uncharged 0/W Interface, 0=135° 4.6 90 Normal Force acting on an Isolated Constant Potential Cylinder (with variable na) Embedded in an Uncharged 0/W Interface 4.7 93 Normal Force Acting on Constant Potential Cylinders Embedded in an Uncharged 0/W Interface 4.9 91 Normal Force acting on an Isolated Constant Potential Cylinder (with variable 6) Embedded in an Uncharged 0/W Interface 4.8 86 94 Parallel Interaction Force acting between Two Parallel Cylinders Embedded in an Uncharged 0/W Interface 98 4.10 Parallel Interaction Force acting between Two Parallel Cylinders Embedded in an 0/W Interface with a Constant Potential 101 4.11 Potential Distributions for Constant Potential (top plot) and Constant Charge (bottom plot) Interactions between Pairs of Identical Cylinders Embedded in an Uncharged 0/W Interface 105 4.12 Parallel Interaction Force acting between Two Constant Charge Parallel Cylinders Embedded in an Uncharged 0/W Interface ix 107 4.13 Normal Force acting on an Isolated Constant Charge Cylinder (with variable KCL) Embedded in an Uncharged 0/W Interface 109 4.14 Normal Force acting on an Isolated Constant Charge Cylinder (with variable 6) Embedded in an Uncharged 0/W Interface 109 4.15 Normal Force acting on Constant Charge Cylinders Embedded in~an Uncharged O/W Interface 110 4.16 Normal Force acting on an Isolated Constant Potential Sphere (with variable na) Embedded in an Uncharged 0/W Interface 115 4.17 Normal Force acting on an Isolated Constant Potential Sphere (with variable 6) Embedded in an Uncharged O/W Interface 116 4.18 Normal Force Acting on a Constant Potential Sphere Embedded in an Uncharged O/W Interface 117 4.19 Change in Free Energy for a Constant Potential Sphere Embedded in an Uncharged O/W Interface 120 A.l Two Cylinders at an O/W Interface 142 A. 2 Co-ordinate System 143 B. l 146 Quadratic Order Triangular Elements D.l Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (ip — 2) Embedded in an Uncharged O/W 0 Interfa ce 187 D.2 Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (tl> 0 = 2) Embedded in an Uncharged O/W Interfa ce 188 x D.3 Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (V'o = 2) Embedded in an Uncharged 0/W Interfa ce 189 D.4 Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (V'o = 4) Embedded in an Uncharged 0/W Interfa ce 190 D.5 Normal Force Acting on Constant Potential Cylinders (V'o = 4) Embedded in an Uncharged 0/W Interface 191 D.6 Parallel Force Acting between Two Constant Potential Cylinders (V'o — 2) Embedded in a Constant Potential 0/W Interface 192 D.7 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 2) Embedded in a Constant Potential 0/W Interface 193 D.8 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 2) Embedded in a Constant "Potential 0/W Interface 194 D.9 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 2) Embedded in a Constant Potential 0/W Interface 195 D.10 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 2) Embedded in a Constant Potential 0/W Interface 196 D.ll Parallel Force Acting between Two Constant Potential Cylinders (V'o = 2) Embedded in a Constant Potential 0/W Interface 197 D.12 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 4) Embedded in a Constant Potential O/W Interface 198 D.13 Parallel Force Acting between Two Constant Potential Cylinders (V'o = 4) Embedded in a Constant Potential 0/W Interface 199 D.14 Normal Force Acting on a Constant Potential Sphere (V'o = 4) Embedded in an Uncharged 0/W Interface 200 xi List of Tables 4.1 Force of Repulsion (/"*) Between Two Spherical Particles with Equal Constant Surface Potentials in a Single Electrolyte 4.2 72 Double Layer Free Energy (G'/ATT) of Two Spheres with Equal Constant Surface Potentials in a Single Electrolyte 4.3 74 Free Energy (G") of Unit Cell with a Constant Potential Sphere (ip = 2) 119 Q D.l Force of Repulsion (/") Between Two Cylindrical Particles (/co. = 10, tpo = 2) Embedded in an Uncharged O/W Interface xn 186 Acknowledgement Many thanks are due to Dr. Bruce D. Bowen who provided valuable suggestions throughout this work and who took great care in criticizing the earlier versions of this thesis. I would also like to thank Dr. Samuel Levine for his patience and the wisdom he shared with me. Thanks are also due to David G. Taylor who provided stimulating conversation regarding the fundamentals of the finite element method. Finally, I wish to acknowledge the Alberta Oil Sands Technology and Research Authority for partial funding of this work. xiu Chapter 1 Introduction The stability of oil-in-water or water-in-oil emulsions may be greatly enhanced by the presence of fine particles in the mixture. This phenomenon is of obvious concern to those involved in any liquid-liquid-particle separation process such as oil recovery or tar sand processing. Accordingly, there has been a wealth of work done to investigate this phenomenon but the majority of it is empirical and process specific. A complete, fundamental understanding of the mechanisms responsible for the stabilization phenomenon is lacking. A comprehensive study of the mechanisms governing the stability of oil-in-water emulsions in the presence of fine particles is being conducted by Levine et. al. [1] and this thesis details the contributions of the author to that study. Particle-stabilized emulsions have been recognized and studied since the beginning of this century. The first detailed experimental study was performed by Pickering [2] in 1907 and consequently particle-stabilized emulsions are often called Pickering emulsions. Since that time, subsequent investigations have revealed that two conditions must be met for Pickering emulsions to arise : (1) the particle size must be smaller than the droplet size [3, 4, 5], and (2) the particle surface must be wettable by both the oil and the water phases. The second condition requires that the three phase contact angle, 6, (measured across the water phase by convention) must be significantly greater than 0° and significantly less than 180 . Finkle et al. [6] were the first to suggest that the type of emulsion formed depends on the contact angle; they stated that the better wetting liquid should constitute the continuous phase. From subsequent work [4, 7, 8] it is now 1 Chapter 1. Introduction 2 generally assumed that oil-in-water emulsions result when 6 is below 90° and water-in-oil emulsions result when 6 is above 90°. Levine and Sanford [9] performed a thermodynamic analysis of particle partitioning between a continuous aqueous phase and a planar oil/water (O/W) interface that indicates a spherical particle with an appropriate 6 is trapped in. a deep energy well ( O[10 6 kT]) at the surface of an oil-in-water emulsion droplet. A similar result was calculated by Pieranski [10] for polystyrene particles trapped at a water/air interface. Since the adsorption energy is so large, Levine [1] proposed that the 0 / W interface should be covered by a close-packed structure of particles. This prediction has been verified experimentally for various liquid-liquid-particle systems [1, 11, 12]. It is believed that this close-packed film of particles, because it projects predominantly into the continuous phase due to contact angle considerations, acts as a steric barrier to coalescence. If the close-packed structure of particles is to be stable at equilibrium there must be a net repulsive interaction between particles in the structure which counteracts the strong attraction of the particles to the interface. To theoretically explain the existence of this film, the contributions to the net particle interaction must be characterized and shown to balance the adsorption energy. The net interaction between particles in the structure is thought to be a result of [13, 15] electrostatic repulsion, short-range solvation and solid-elastic repulsions, van der Waals attraction and a capillary interaction. Levine et. cd. [1, 13, 14] have made preliminary investigations into all of these interactions but have yet to satisfy the energy balance implicit for stability at equilibrium. This thesis provides a numerical treatment of the equilibrium electrostatic interaction to establish whether the approximations made in the preliminary analysis were valid. The second chapter of this thesis presents the theory used to determine the electrostatic interaction between particles adsorbed at an 0/W interface. Two separate cases are considered: i Chapter 1. Introduction 3 • two identical, parallel, infinitely long, solid cylinders at a given separation where each is embedded to the same extent i n an oil/water interface, and • a spherical particle adsorbed at the 0 / W interface i n a two dimensional closepacked structure. Mathematical descriptions are given for both models and the basis for the governing equations and boundary conditions is examined. In order to test the solution methodology the popular case of two identical spheres completely immersed i n the aqueous phase is also considered. T h e electrostatic potential distribution i n the aqueous phase is described by the Poisson-Boltzmann equation ( P B E ) which has received much attention i n the field of colloid stability. A summary of relevant solutions to the P B E is given with an emphasis on the numerical work performed to date. T h e second chapter also includes a review of the finite element method ( F E M ) which is used i n the present study to solve the governing equations. Contrary to popular belief numerical analysis does involve experimentation; tests are required to determine the convergence characteristics of a given solution. T h e design of the finite element grid used to physically represent the emulsion models and the dependence of solution accuracy on grid parameters is discussed i n the third chapter. T h e results and discussion are presented i n the fourth chapter. A comparison between previous solutions to the P B E for two spherical particles is made with the results of this thesis. A n extensive study of the cylindrical and spherical particle systems is made and the dependence of the double layer interaction on model parameters (particle separation, 6, particle radius, particle surface potential, etc.) is revealed. Finally, the conclusions and recommendations for future work related to this thesis are given i n the fifth and sixth chapters, respectively. Chapter 2 Theory 2.1 Electric Double Layer Surrounding an Isolated Charged Solid Most solids acquire a surface electric charge when they are in contact with an aqueous medium; the possible charging mechanisms are • ionization of surface groups of the solid, • unequal dissolution of ions from the solid, and • preferential adsorption of ions from solution. For the purpose of this thesis it is not necessary to know the details of these mechanisms. However it is important to note that if the surface charge is a result of preferential adsorption of ions then the surface potential will remain constant; but if the surface charge is a result of ionization the surface charge will remain constant when the surface interacts with a second charged surface [16]. To preserve the electroneutrality of the system, ions of opposite charge (counter-ions) concentrate in the solution close to the surface while ions of like charge (co-ions) are repelled away from the surface. The charged surface and solution comprise the electric double layer, which is conventionally assumed to consist of an inner layer, where the ions are associated so strongly with the surface that they are unaffected by thermal agitation, and a. diffuse outer layer, where the equilibrium distribution of ions is assumed to be determined by electrostatic and thermal motion (i.e. diffusion) effects. 4 Chapter 2. Theory 5 For a rigorous mathematical treatment of electric double layer phenomena it is appropriate to disregard the various quabtatitive models of the inner layer and only consider the diffuse outer layer. The simplest quantitative description of the diffuse outer layer, and the one that is employed in this thesis, is known as the Gouy-Chapman model which is based on the following assumptions: • the solvent is assumed to influence the diffuse layer only through its dielectric constant which is assumed to be constant throughout the diffuse layer, • the ions in the diffuse layer behave as point charges distributed according to a Boltzmann distribution, and • the electrostatic potential distribution is governed by the Poisson-Boltzmann equation (PBE). The PBE is obtained by considering the differential form of Gauss' law which relates the divergence of the electric field to the charge distribution in a given volume. This relationship is given by Poisson's equation as V * = -Airp/e (2.1) 2 where $ is the electrostatic potential, p is the charge density of the medium, e is the dielectric constant of the medium and V is the Laplacian operator. Arbitrarily, the solid 2 surface charge is assumed positive; the results are applicable to negatively charged solids by reversing the charge signs. The local charge density of the medium is determined by the local ion concentrations. Thus, assuming a binary electrolyte, the charge density is given by p = e(u z + + — v^zJ) (2-2) where u and i / _ are the local concentrations (number/cm ) of the co-ions and counter3 + ions respectively, z± are the ion valences and e is the elementary charge. Chapter 2. Theory 6 The local ion concentrations are described by Boltzmann's law which assumes a statistical equilibrium based on Coulomb and thermal energies, i.e. v±=n exp(T^-) (2.3) ± where n± are the concentrations of the ions in the bulk solution, T is..,.the absolute temperature and fc is Boltzmann's constant. The bulk solution must be neutral; so the individual electrolyte concentrations are related to the overall electrolyte concentration by nz + = n_z_ = CN /1000 + (2.4) a where C is the electrolyte concentration (equivalents per litre) and N number. a is Avagadro's Thus, by making the appropriate substitutions, equation (2.1) may now be expressed in the form known as the Poisson-Boltzmann equation (PBE) „ , T V * = 2 -AireCN -z+eV. ar ^ elOOO exp( n n n 1 + kT F V ^-e* x l ) - exp( — — ) ' kT v y n (2.5) ' v The dimensionless form of equation (2.5) is obtained by invoking the reduced quantities used by Loeb et al. [17]. The reduced potential, if;, is defined as i> = (2-6) w A general reduced distance co-ordinate, a, is defined as where a is any displacement co-ordinate (x,y,z,r) and where lOOOefcr v ; and 2 A _ z +z_ + -~2zT~ (2.9) Chapter 2. Theory 7 The parameter K incorporates several properties of the double layer and the quantity 1/K is often referred to as the double layer thickness, i.e. the distance over which the potential decreases by an exponential factor at low potentials. By substituting equations (2.8) and (2.9) into (2.7) it is apparent that the reduced distance is independent of co-ion concentration; this allows results with the same C and z_ but different z + to be compared conveniently. The dimensionless form of the P B E is = *p(*-V0- *p(-*+V0 e e ( 2 . 1 0 ) where the V is expressed with respect to the reduced co-ordinates. The right-hand side 2 of equation (2.10) can be linearized by expanding the exponential terms and retaining the first order terms. This is often called the Debye-Huckel (DH) approximation and it is only valid when z±tfj <C 1 which corresponds to <C 25mV at room temperature. The electrostatic potential distribution surrounding a body immersed in an aqueous medium may be obtained by solving equation (2.10) subject to appropriate boundary conditions. If the medium can be assumed to be infinite then one type of boundary condition is the usual, arbitrary assumption of electrostatic theory: \£ has a value of zero at infinity. A second type of boundary condition applies at the solid surface where it is usual to assume either a constant surface potential or a constant surface charge. The choice of solid surface boundary condition depends primarily on the charging mechanism as noted above. More complex surface boundary conditions would have to be employed if a dynamic study was being performed, but that is beyond the scope of this thesis. A third type of boundary condition normally arises through symmetry considerations. Owing to the importance of double layer effects in the theory of colloid stability, the solution of the P B E has received constant attention since the turn of the century. Closedform analytical solutions exist for just a couple of limiting cases. Overbeek [18] provided a solution for the case of a plane of infinite extent with a constant surface potential in an Chapter 2. Theory 8 infinite aqueous medium. This result may be applied to curved surfaces where na ^> 1 (na is the reduced radius of curvature) since the particle is large compared to the double layer thickness and planar geometry may be assumed. As detailed below, this approximation is useful when employing the Derjaguin [19] method to solve for the interaction energy between two spheres. If the DH approximation applies, an exact solution to the PBE may be obtained for a spherical particle [16]. For the majority of colloidal systems, however, neither the DH approximation nor the na ^> 1 condition applies, so approximate methods must be used to solve the PBE. Most model colloidal systems consist of spherical particles, so the solution of the PBE for an isolated sphere assumes a special importance. The first numerical solution was provided by Muller [20]; but without the aid of a computer this amounted to a cursory treatment. Hoskin [21] used an integration scheme similar to Miiller's to compute the potential distribution for a limited number of cases of a sphere in a binary symmetrical electrolyte. Loeb et al. [17] improved on Hoskin's results by using an adaptive technique that adjusted the integration interval according to the slope of the potential curve which becomes very steep near the surface in cases of high surface potential. By employing equation (2.9), Loeb was also able to account for unsymmetrical electrolytes. There have also been a number of approximate analytical solutions presented based on perturbation [22],[23],[24],[25] and variational [26] methods. The most practical approximate analytical solution is provided by White et al. [27] who use a change of variable technique that assumes planar geometry very close to the particle surface but allows for the correct asymptotic behaviour far away from the surface. This solution shows excellent agreement with Loeb's [17] numerical results except for cases of small reduced radius and/or large reduced surface potential. White's method [27] yields an accurate surface charge density/surface potential relationship and it may be readily applied to other geometries such as cylindrical particles. His expression relating the reduced surface charge, cr", to the 9 Chapter 2. Theory surface potential, ipoo, for cylinders at infinite separation in a single electrolytic medium is /3~ - 1 or* = 2 s n ih ( V > o o / 2 ) , \| cosh (VW4) 2 1 + 2 ( 2 ' U ) where p = ( 2 1 2 ) In equation 2 (1 .1 ) K and K~i are the zero andfirstorder Bessel functions of the second 0 kind and a is the reduced radius of the cylinder. The corresponding expression for spheres at infinite separation in a single medium is , 81n[cosh(V.„/4)] , a* = 2 s n ih V (W 2 ) 2 1 N a cosh ('0cx)/4) a sinh('0oo/2) where a is the reduced sphere radius in this case. Both equations (2.11) and (2.13) were 2 2 2 used in this work to set constant surface charge boundary conditions. Before discussing how they were used, it is best to clarify the relationship between the surface potential and the surface charge that was assumed in this work when comparisons were made between constant potential and constant charge interactions. The reduced charge on the interface between two media (1 and 2) is given by Gauss' law as 47T On 47r on where the derivative of the potential is taken normal to the interface. If an isolated particle (say medium 2) has a constant surface potential over its whole surface, then the potential within that particle must be constant and at the same value as the surface potential (assuming that there is no charge source or sink within the particle) and thus thefirstterm of the right-hand side of equation (2.14) vanishes. This was the assumption made when equations (2.11) and (2.13) were derived. It should be noted that by specifying a constant potential boundary condition under this assumption, a constant surface Chapter 2. Theory 10 charge will also be present on the particle as long as it remains isolated (i.e. at infinite separation from all other particles). The two simplest ways of analyzing double layer interactions between two particles is to assume that either the surface potential remainsfixedor the surface charge remains fixed as the particles are brought together from infinite separation. When comparing these two scenarios it is common to specify a surface potential at infinite separation, V'oo, and to then hold either V'oo or the corresponding cr' constant for various separations of the particles. Thus equation (2.13) was used to predict the surface charge for a given ipoo when constant surface charge interactions between spheres in a single electrolytic medium were evaluated. The constant charge interactions between cylinders considered in this thesis involve cylinders at an O/W interface that are only charged on the surface immersed in the aqueous phase. Equation (2.11) cannot be applied to an isolated cylinder in this situation since the potential will now vary within the cylinder and the full form of equation (2.14) applies. We may justify the application of equation (2.14), however, if we choose to compare cylinders which are charged over their whole surface before they reach the interface; in this case equation (2.11) does apply. It is further assumed that the surface which is immersed in the oil phase becomes uncharged while the charge on the surface in contact with the aqueous phase remains fixed. In short, one assumes that the constant charge boundary condition applied to a particle in a single medium will still apply once that particle is adsorbed on the interface with the one exception that the surface in contact with the oil phase is no longer charged. As will be revealed later in this chapter, thefiniteelement method requires that a constant surface boundary condition be given in terms of equation (2.14), i.e. a difference between the charge fluxes on either side of the interface. Thus, if this difference is calculated for the situation of the isolated particle in a single medium, this same difference 11 Chapter 2. Theory may then be applied to the particle at the 0 / W interface if the above assumptions are made. This is the simplest way to compare constant charge and constant potential interactions between particles at an 0 / W interface. Unfortunately the above discussion necessitated jumping ahead in the subject matter; it is appropriate now to digress and discuss electric double layer interactions in more detail. 2.2 Electrostatic Interactions between Charged Surfaces in an Electrolyte Solution The problem of solving for the double layer interaction between two or more bodies presents further complications. As discussed below, calculating the interaction requires knowledge of the potential distribution around the bodies and hence a solution to the PBE. Unless very simple geometries and symmetric surface boundary conditions are assumed, it is impossible to obtain a closed-form. anlaytical solution to the PBE and therefore approximate methods must again be employed. Before describing the available analytical and numerical solutions to the PBE for two (or more) body systems, the methods for calculating the double layer interaction will first be discussed. There are two ways to characterize the double layer interaction between two charged bodies; by the force acting between the bodies at a given separation or by the free energy change associated with bringing the two bodies to a given separation from an infinite separation. The two quantities are related in that the free energy may be obtained by integrating the force from infinite separation to the actual particle separation, or conversely, the force may be calculated by differentiating "the free energy with respect to the separation. The free energy of a single double layer can be found by considering an imaginary 12 Chapter 2. Theory charging process. If the particle surface charge, <r, remains constant then the free energy is given by G = I" V {a)doJo (2.15) 0 where \&o is the surface potential. When the surface potential is low then cr is proportional to \£n and the double layer free energy is \o~^o [17]. For a constant surface potential boundary condition the free energy is given by (2.16) G = - [*° o-{y )dVo Jo 0 Levine [29] obtained an alternative form of equation (2.16) which is expressed in terms of the electrostatic energy and the excess osmotic pressure in the double layer and is also valid under conditions equivalent to a constant surface potential boundary condition, G = -W - f AUdv (2.17) Jv where W — ^ J \ V^\ dv is the ordinary field energy in electrostatics, All is the dif2 v ference between the osmotic pressure in the double layer and the bulk, and v is the electrolyte volume being considered. The W term represents the work done by distributing the charged ions in the potential field while the osmotic term represents the work done against thermal motion by concentrating the ions. Equation (2.17) is more easily handled by the numerical method considered in this thesis because it only requires that the potential distribution be known for the system. Equations (2.15) and (2.16) require that either an analytical relationship between tpo and a be available or that relationship be calculated at each stage in the charging process. The former requirement is unavailable for the systems under consideration here. The latter would entail solving the governing equations for each stage in the charging process, which is grossly impractical. For two identical particles separated by a distance 2d, the energy of repulsion, AG, is the difference between the free energies per particle, G and G^, at 2d and at infinity d Chapter 2. Theory 13 respectively, i.e. AG = 2(G -G ) d (2.18) 00 where the factor of 2 arises because two particles are being considered. Calculating A G from equation (2.18) involves taking the difference of two very large numbers to yield a much smaller number. The energy of repulsion is usually less than 10% of 2Gd for typical interactions [28]. Thus, A G calculated from equation (2.18) will be very sensitive to the error associated with 2Gd and 2Goo. Hoskin and Levine [28] developed a direct method for calculating the force between two particles based on the electrostatic stress and the osmotic pressure associated with concentrating the ions between the particles. Using the general considerations of electrostatics and employing the PBE one obtains v - r = P E = — v * v * = v n (2.19) 47T where T is the Maxwell stress tensor, E = — V $ is the electric field and II is the osmotic pressure. The Maxwell stress tensor is of rank two; each of its rows can be considered a vector whose divergence gives the force density due to electrostatic forces. A useful notation to express the components of T is [30] r = -^(EE - y J) (2.20) where EE is the dyadic product of the electric field, I is the unit tensor, and E is the magnitude of the electric field. Since the electrostatic force is given as a divergence of a vector, one can make use of the divergence theorem to calculate the net electrostatic force over a large volume. This is the advantage of using the Maxwell stress tensor because it allows the force acting on a charged volume to be calculated in terms of the electric field at the surface bounding the volume, irrespective of the behaviour of the electric field within the volume. 14 Chapter 2. Theory Therefore, by using the divergence theorem and equation (2.19) we have 0 = / ( V • T - VU)dv = I \-^{{E Jv Js4.iv • n)E - \ E U } - (II - U )n]dS 2 (2.21) 2 0 where S is a closed surface bounding the volume v, n is the unit outward normal vector to the surface 5, and TJ. is the osmotic pressure of the bulk solution. The I I / ndS = 0 0 0 5 term is added onto equation (2.21) so the osmotic contribution goes to zero in the bulk. To calculate the force between two spheres, Hoskin and Levine consider a hemispherical volume which is bounded by the following surfaces S is the surface surrounding one of the particles, p Sm is the median plane between the two particles, and Sh is the surface of the hemisphere. As the radius of the hemisphere goes to infinity the integral over Sh will tend to zero since the potential at this surface approaches zero. If the integration of equation (2.21) is carried out over the separate surfaces S and 5 each result yields a vector that points p m along the line of centres joining.the particles / (4~E J S M 2 + AIl)dS = I S7T JS,, {-^-E 2 + AU)ndS (2.22) 07T and since ideal behaviour is being assumed for the electrolyte, the osmotic pressure term is given by [16] An = fcT[(i/ -n ) + + + (i/_-ii_)] (2.23) which depends on the overall electrolyte concentration and the electrostatic potential. Note that this problem is axisymmetric and consequently no forces act on the particles which are perpendicular to the line of centres joining the two spheres. The left-hand side of equation (2.22) is the total force acting across the median plane due to the electrostatic stress and osmotic pressure and the right-hand side is the Chapter 2. Theory 15 force acting across the particle surface. Hoskin and Levine proved that equation (2.22) does actually represent the force acting between particles by differentiating afreeenergy expression with respect to separation to obtain the same result. Both the free energy and force approaches mentioned above have been used to calculate double layer interactions between particles. The majority of the analytical solutions attempt to solve for the free energy directly while the numerical treatments usually find the force and then integrate to get the free energy. The reasons for this situation will be revealed as the various analytical and numerical treatments for interacting particles are reviewed. Verwey and Overbeek [18] derived an exact solution to equation (2.18) for the case of two parallel plates with constant surface potentials and a symmetrical electrolyte. The exact expression contains elliptic integrals, the values of which need to be obtained from function tables, so Levine and Suddaby [31] derived approximate expressions valid for moderate potentials and small and large separations and they also treated the case of asymmetric electrolytes [32]. Derjaguin [33] has shown how the interaction of two spheres can be calculated from the interaction of two parallel plates. In the Derjaguin approximation for spheres, concentric, parallel rings are paired off from each sphere and the parallel plate solution is used to obtain the double layer interaction between the pairs of rings. The contributions from all the rings are summed to give the total double layer interaction between the spheres. A more detailed description of the method will not be given here as it is generally available in any comprehensive text on colloid chemistry (the Derjaguin approximation applied to cylinders embedded in an 0 / W interface is provided later on in this chapter). This approximation is only valid when the double layer thickness and the particle separation are similar but are both small compared to the particle sizes. Applying the DH approximation, McCartney and Levine [34] used a new approach to solve for the interaction energy between two identical spheres with constant surface Chapter 2. Theory 16 potentials which involved integrating an expression that related the surface potential to a surface distribution of electric dipoles. This technique yielded improved solutions over those predicted by the Derjaguin formula and their general formula actually reduced to the Derjaguin expression when the separation between particles was much less than the particle radii. Bell et al. [35] extended this approach to cases involving dissimilar spheres with unequal, constant surface potentials. For larger separations, Bell et al. [35] developed an approximation to the interaction energy based on the linear superposition of the potentials of the overlapping double layers. Healy et al. [36, 37] used the DH form of the PBE to derive energy expressions for asymmetric interactions (different sphere sizes) and for both constant surface charge and constant surface potential cases (surface charge remains constant but not necessarily equal on both spheres or surface potential remains constant but not necessarily equal on both spheres during the interaction). Bell and Peterson [38] studied the effects of either constant potential or constant charge interactions between dissimilar particles using the Derjaguin approximation at small separations and the superposition approximation at larger separations. While much work has been published on approximate analytical solutions for the free energy of interaction of different particles with different boundary conditions, this work is beyond the scope of this thesis but an excellent review of the current methods is provided by Overbeek [39]. To complement, the wealth of analytical approximations to the double layer interaction between two particles there has also been a steadyflowof numerical solutions presented; in fact it is common practice to test the validity of approximate analytical expressions against the more exact numerical solutions. Since the test case of two identical spherical particles is numerically evaluated in this thesis, a more detailed review will be given of previous numerical solutions to this problem. Hoskin and Levine [28] (HL) were the first to publish a thorough numerical analysis of Chapter 2. Theory 17 the double layer interactions between identical spheres with constant surface potentials in a single electrolyte. They solved the P B E by the finite difference method using a coarse grid (144 points) and a dipolar co-ordinate transformation. The dipolar transformation concentrates the nodal points in the region between the spheres but the nodal separations on the 'outer' side of each sphere are~so great as to severely limit the accuracy of the solutions in that region. HL only retained first order finite difference approximations and they assumed their solutions had converged when the iterative corrections to the potentials did not affect the third decimal place. The force between the particles was calculated using both the right-hand and left-hand sides of equation (2.22). It was found that the results obtained by integrating over the particle surface were very inaccurate because of the steep potential gradient that existed there. H L also calculated the double layer interaction energy numerically by four separate methods. Two involved numerically integrating, with respect to separation, the force results just mentioned. The superposition principle was used to predict the potential distribution and hence the force between the spheres at large separations thereby obviating the need for computations in this region. Again, the free energy calculated using the median plane force results were superior to those using the particle surface force. The other two methods involved expressing the free energy of the double layer as: 1) a volume integral over the complete double layer, equation (2.17) and 2) a volume integral over the particle surface, equation (2.16). The second method suffers from the same problems as the force calculation performed over the particle surface; numerical differentiation of steep gradients is prone to large errors. The dipolar transformation introduces inaccuracies into both methods as well since the nodes are so unevenly spaced about the particle. Another factor that limits the accuracy of the direct free energy calculations is that equation (2.18) has to be used to calculate the net interaction energy. This involves taking the difference of two large quantities to get a much smaller quantity. One of the main conclusions of the HL paper was that Chapter 2. Theory 18 integrating the force that was obtained from the left-hand side of equation (2.22) yields the most accurate value for the free energy of interaction. McCartney and Levine [34] (ML) repeated HL's calculations using a finer grid (1200 points) and using the dipolar transformation with an added transformation that served to concentrate more nodes on the outer side of the sphere where the double layer interaction is least and the potential gradient is steepest. This added transformation made the thirdorder differences sufficiently small on the outer side that three to four figure accuracy was claimed for the potential values. ML calculated the force of the interaction using the midplane method and integrated this force in the same manner as HL to get the interaction energy. In nearly all cases, ML's force and energy results were 2-5% lower than those of HL. ML used these figures to check the validity of their new analytical technique described above. Honig and Mul [40] presented tables of the energy of repulsion for parallel plates and two spheres both at constant equal potential and at constant equal surface charge. The Derjaguin method was used to treat the case of the interacting spheres so these results only apply to spheres with radii much greater than the double layer thickness. Philip [41] used a cell model to obtain double layer interaction energies for a sphere with a constant surface potential in two and three-dimensional homogeneous particle swarms. Philip used the DH approximation and extended it to cases of higher surface potential using a 'matched approximation'. This involved expressing the double layer potential for a given separation as the DH result multiplied by the ratio of the exact/DH results for a single particle at infinite separation. The numerical results of Loeb et al. [17] were used as the exact values. A comparison of the three-dimensional results with HL's values for two interacting spheres was made by assuming twelve two-particle interactions, superimposed, represented the interaction between a particle and its enviroment. As could be expected the comparison was poor due to the drastic assumptions made. 19 Chapter 2. Theory McQuarrie et al. [42] solved for the double layer force between spheres that had equal constant surface potentials or surface charges. Rather than solving directly for the total electrostatic potential, McQuarrie et al. solved a differential equation for the correction to an approximate solution for the PBE. The approximate solution was obtained using the superposition principle and the one particle approximation to the PBE given by Brenner and Roberts [26]. The differential equation was solved by finite differences (2500 points) in dipolar co-ordinates. The force was calculated from the left-hand side of equation (2.22). The results for constant potential spheres agreed to within 1% of ML's results. The important advantage of McQuarrie's method is that it allowed the force between spheres of equal constant charge to be calculated. When implementing the constant charge boundary condition, McQuarrie ignored the potential distribution inside the sphere. This is an acceptable approximation if one is considering two, interacting, semi-infinite plane walls which is the assumption made in the Derjaguin approximation. McQuarrie found that his constant charge results for large spheres showed good agreement with the Derjaguin approximation providing the proper choice of surface potential at infinite separation was made. The Derjaguin approximation for the reduced double layer interaction force, /", between two spheres is given [38] by (2.24) /" = 2-KKO,V (S) + where a is the sphere radius, S is the closest approach of their surfaces and V (S) is + the parallel plate interaction energy per unit area which is given by equation (5.15) in reference [38]. The y (S") term is a function of the surface potential at infinite separation, + and the surface potential at separation S, tjj. The surface charge, cr*, is related to ipoo through the relationship cosh V'oo = 1 + cr" /2 2 (2.25) and ip is implicitly dependent on ip^ and S through a complicated expression involving Chapter 2. Theory 20 elliptic integrals (see equation (6.4) in [38]). Usually, when using the Derjaguin formula for constant charge interactions, the surface potential at infinite separation is calculated for a plane wall; McQuarrie revealed that the surface potential calculated for a plane wall at infinite separation, rather than a sphere, indeed yields a better approximation, i.e. closer to his numerical results. Constant charge interactions are considered in this thesis where the interior of the particle is ignored (as McQuarrie did) and also where the potential distribution in the interior of the particle is simultaneously solved with the aqueous phase potential distribution; the results of these calculations are given in Chapter 4. Chan and Chan [43] used the finite element method to obtain the force of repulsion between two equal spheres. They considerd two cases: 1) constant equal surface potentials and 2) constant and equal surface potentials which are opposite in sign. Again, the force of the interaction was calculated using the left-hand side of equation (2.22). The motivation for this work was to confirm the validity of the approximate expressions derived by Healy et al. [36] using the Derjaguin method which were challenged by Barouch et al. [44]. The results of the numerical calculations supported Healy's expressions. Unfortunately no details of the numerical approach were given in the paper and the force results were just presented in graphical form. To check that thefiniteelement method developed in this thesis is correct, the case of two equal spheres with constant, equal surface potentials or surface charges is considered. The force of repulsion between the spheres is calculated using the left-hand side of equation (2.22). Specific details of the model used will be given later in this chapter in Section 2.4. Chapter 2. 2.3 21 Theory Electrostatic Interactions between Particles at an Electrolyte Interface When particles are trapped at an interface between an aqueous electrolyte and a second non-polar fluid there will be an electrostatic interaction acting through the non-polar phase in addition to the double layer interaction acting through the water phase. The particles are assumed to be charged only on the part of their surfaces immersed in the water phase. The distribution of ions in the water phase is still governed by the PBE so the interaction from the double layer overlap will decrease approximately exponentially with increasing separation between particles. There is no similar screening effect in the non-polar phase however so one would expect the interaction to be more long range. A particle at such an interface has an asymmetric charge distribution and therefore a dipole moment. Pieranski [10] was the first to suggest that the particles could interact through the non-polar phase according to a dipolar force law. Nichols and Pratt [45] derived the interaction energy for two point charges located at an equal distance from an air/water interface that decayed as r~ where r is the charge separation, which is large 3 compared to the double layer thickness. Hurd [46] considered charged particles at an air/water interface using an integral expression obtained by Stillinger [47] based on the linearized PBE. In addition to the exponentially decaying double layer interaction, a long range dipole-dipole interaction was found that was described by an expression similar to Nichols'. Jancovici [48] proved that pairwise interactions, parallel to the interface, between charged particles decay as r~ if the particles are near the interface. 3 Levine et. al. [1] used a modified version of Nichols' expression to estimate the long range dipolar interaction energy of a charged particle in a close-packed hexagonal array. They found that the dipole repulsion was less than 10~ of the ordinary double layer 3 repulsion. Solving for the double layer interaction between two spheres at an 0/W interface is 22 Chapter 2. Theory a three-dimensional problem since there is no longer axial symmetry about the line of centres joining the two spheres. The problem may be reduced to a simpler integral form by using the Derjaguin method. Recall that in this approximation parallel, concentric rings are paired off from each sphere and the force is calculated using the analytical parallel plate solution. This approximation is valid when the particle radius is large compared to: 1) the separation of the particles and 2) the double layer thickness. Levine et. al. [1] used the Derjaguin approximation to estimate the double layer interaction between particles with equal, constant surface potentials at an O/W interface where only the portion of the particle immersed in the water phase contributes to the interaction and the interface is assumed to be uncharged. The additional assumptions implicit in their solution were that the interface was planar (i.e. oil droplet diameter particle diameter) and that the degree of immersion of the particle was determined solely by its contact angle (i.e. gravitational effects were ignored). The Derjaguin method can also be applied to two equal, parallel cylinders. Instead of pairing off parallel rings, one pairs off parallel strips from each cylinder. The free energy of the double layer is given by V (H )= e 0 f J— a V (H)dS p (2.26) cos 8 where V (H) is the energy expression corresponding to two parallel plates, Ho is the P smallest separation between the cylinder surfaces whose axes are separated by R, a is the cylinder radius, 8 is the contact angle, and H is the separation of the surfaces at a height h from the line joining the two cylinder axes and is given by H — R — 2y/a — h 2 (see Fig. 2.1). In the limiting case of low surface potential and a 1:1 electrolyte lp is given as [19] V (H) = ^^5[l-tanh(l/2/tfr)] r P 2 23 Chapter 2. Theory Figure 2.1: Derjaguin Method Applied to Equal, Parallel Cylinders = ^^{[1 - tanh(l/2 (# - 2vV - h )}} 2 K (2.27) 47T To get an expression for the force between the two cylinders we differentiate the energy expression with respect to separation and obtain r f (H) = c 87T ch [l/2 (R 2 se K 2y/a?-W))dh (2.28) J-a cos 8 where f (H) is the force per unit length of the cylinder. Equation (2.28) must be inc tegrated numerically and upon doing so it is found that most of the interaction comes from the double layers overlapping in the region about the plane joining the axes of the particles. Note that this method completely ignores any interaction through the oil phase or the possibility of having a charge on the 0/W interface. 2.4 Model Descriptions As noted in the introduction, three separate double layer interaction models are considered in this thesis. The familiar case of two spheres fully immersed in an electrolyte Chapter 2. Theory 24 (Model I) is an axisymmetric problem where only the PBE needs to be solved to calculate the interaction. The two equal cylinders partially immersed at an 0/W interface (Model II) represent a two-dimensional problem in Cartesian space. A lone sphere partially immersed in an 0/W interface surrounded by a cylindrical bounding envelope perpendicular to'the interface (Model III) is an axisymmetric problem. In Models II and III, three separate phases are being considered; an aqueous phase, the solid dielectric particle and the oil phase. Laplace's equation governs the potential distribution in the non-polar solid and oil phases while the PBE applies in the aqueous phase. In Model I only the solid and aqueous phases are present. When considering Model I with a constant potential boundary condition, it is not necessary to solve for the potential distribution in the particle. This is because we assume Laplace's equation holds inside the sphere and the whole sphere surface is fixed at the same potential. However, for a constant charge boundary condition, the interior of the sphere must be considered since the constant surface charge is related to the potential distribution in both phases through Gauss's law. Thus the potential distribution is solved simultaneously in both phases subject to the fixed charge boundary condition. For Models II and III we again have situations where surface charge boundary conditions are being applied, e.g. a zero charge boundary condition is always assumed on the particle/oil interface. This necessitates that the potential distribution in the oil and particle phases must always be solved simultaneously. If there is a constant charge boundary condition on either the O/W interface or the particle/water interface then the potential distribution must be solved in all three phases simultaneously, subject to the fixed charge boundary conditions. Before defining the specifics of each model the applicable forms of the PBE will be given; Laplace's equation is obtained by setting theright-handside of the PBE to zero. 25 Chapter 2. Theory The two-dimensional form of equation (2.10) in Cartesian co-ordinates is \ _L 9 9 (, ^ \ d ex P ( - V Q - exp(-2 V>) z + and the axisymmetric form is 1 o dip d dip r or Or oz Oz exp(z-ip) - exp(-2 V>) + Lz_ The reader will notice that the dielectric constant, e, which is customarily included in the definition of K, equation (2.8), now appears in equations (2.29) and (2.30): this is necessary so that the material discontinuities inherent in a multi-phase problem can be accounted for in the numerical treatment. Accordingly, a modified reduced distance co-ordinate is defined by a =^ (2.31) where K = Ky/e^ (2.32) For comparison purposes all distances will still be expressed in terms of double layer equivalents, equation (2.7), so the conversion between the modified reduced co-ordinate used in the governing equations (2.29) and (2.30) and that defined by equation (2.7) is simply a = a/^/i^ where e is the dielectric constant of water. w 2.4.1 T w o Identical Spheres Fully Immersed in an Electrolyte The centres of two identical spheres (radius=a) are fully immersed in an electrolyte solution, and are separated by a distance R. The spheres may have equal, constant surface potentials or surface charges. Consequently, the median plane separating the spheres is a mirror plane of symmetry so we need only consider one of the particles. There is also axial symmetry about the line joining the centres of the spheres so only the upper hemisphere needs to be considered (see Fig.2.2). Chapter 2. 26 Theory r N. V ( Figure 2.2: Model I: Two Identical Spheres in an Electrolyte The smallest distance between a sphere surface and the median plane is d = (R — 2a)/2. This is an axisymmetric problem so equation (2.30) applies. The boundary conditions are 1. ., $ —• 0 as z —» oo and r —> co; Sh- 2. constant potential (9 = or charge (a = a ) on the particle surface; S . 0 p 3. ^ = 0 on the median plane; z = — (d -f a), (S ). m 4. = 0 at r = 0. The force, / , acting between the spheres is calculated using the left-hand side of equation (2.22). It is convenient to define a reduced force which has a similar form to that used by Hoskin and Levine [28]: _ f ~ /c 1000 2 kTCN € a w f (2.33) Chapter 2. Theory 27 In terms of reduced quantities the expression for the reduced force is / . = 2ix(z_ + z ) * 5 + i f°° £„, <9V + i z e,„ ± ) p ^ , Jo ) + A T 1 H l I or Z (2 , 34) where AT is the reduced osmotic pressure given by A T = ^-[exp(-z V>) - 1] + z [exp(z_-0) - 1] + + 2z and where I/J, ^ ^ + and AT are evaluated along the median plane. The free energy of the double layer is calculated using equation (2.17). Again, it is convenient to express this quantity as a dimensionless reduced free energy, G", which is defined as G - = - * ™ „ kTCNjJ B (2.36) 2 In terms of the reduced quantities the expression for the reduced free energy of the double layer is c-=-' '-;' "r i 1 2.4.2 < r&W + xrv**; (2.37) Two Identical Cylinders Embedded in an O / W interface Two infinitely long, parallel, identical cylinders of radius a are situated at a planar 0/W interface, separated by a distance R. Each is immersed to the same depth at the interface. It is assumed that the gravity effects acting on the particle are negligible and the immersion depth is determined solely by the contact angle 8. The part of both cylinders in contact with the aqueous phase may have equal constant potentials or charges. The median plane separating the two cylinders represents a mirror plane of symmetry so only one cylinder needs to be considered. The immersion angle, 8, measures the extent of immersion of the cylinders into the oil phase and is given as (r — 0), where 8 is the conventional three phase contact angle. Again the smallest separation between a cylinder and the median plane is d = (R — 2a)/2 (see Fig.2.3). Figure 2.3: Model II: Two Identical Cylinders at a Planar O/W Interface This problem may be solved in Cartesian co-ordinates so equation (2.29) applies in the water phase and the Laplace form of equation (2.29), right-hand side equal to zero, applies in the oil and solid phases. Note that the appropriate dielectric constant must be used when applying equation (2.29) in each of the three separate phases. The boundary conditions for this case are more complicated since the interfaces between the phases must be considered. These are 1. constant potential (ty = $ ) 0 o r constant charge (a — o- = Q — f^^fS where w and p refer to the water and particle phases respectively) on the particle surface in contact with the aqueous phase; 5,^. 2. | ^ = 0 on the median plane; x = —(ci + a), (S^ + Som)3. no charge on the particle surface in contact with the oil phase; (i.e. e^(lf = 0 Chapter 2. pltr e Theory 29 ° ^ P refer to the oil and particle phases respectively and tp repre- w n e r e a n n sents the potential distribution in the non-polar oil phase). 4. either no charge (e ^f = e ^-) or a constant potential (9 — * ) ° the O/W n w 0 0 interface; S^. 5. if there is no charge on S ow and 5 ^ as z —> ±oo then f -» 0 on the surfaces and x —> oo, but if there is a constant potential on then | ^ = 0 as x — > oo. For the case of a constant potential on the O/W interface, the second part of boundary condition 5 applies because the potential distribution normal to the interface will resemble the distribution for a charged plane wall far from the particle. The force acting on the particle maj be found using the method of Hoskin and Levine r but since three discontinuous media are now present, all interfacial boundaries must be explicitly considered. An interesting result of this model is that there is a force acting normal to the O/W interface (the z direction) on the cylinders as well as a force acting parallel to the line of centres joining the cylinders. This normal force is present even in the case of infinite separation. The force equations are presented below and their full derivation is given in Appendix A. The force acting parallel to the line of centres will have contributions from the aqueous as well as the oil phase. The contribution from the aqueous phase is . Z j j ^ j t y & r (2.38) + LTW and from the oil phase £ = ( e w ! z ± J 2 dz i ± ) ( 2 . 3 9 ) where e is the dielectric constant of the oil phase. The above equations represent the D force per unit length acting on one of the cylinders. Chapter 2. Theory 30 The normal force acting on the particle is obtained by integrating over the O/W interface 5„ 'OW f; = where ^I ^^ i sr {j ^jp is the normal derivative of the { +c (? f]+r}di (24D) potential into the oil phase. If there is a constant potential boundary condition on the O/W interface then the normal force is given by 2.4.3 Isolated Sphere Embedded in an O/W Interface in a U n i t Cell This model. attempts to simulate the three-dimensional problem of an isolated sphere in a close-packed, two-dimensional array of identical spheres by using a two-dimensional unit cell technique. The unit cell technique is based on the concept that a collection of particles can be divided into a number of identical cells, one particle occupying each cell [49]. The boundary value problem is solved for a single cell and, for this application, the resulting energy is used to approximate the double layer energ}' per particle in the array. In general, a unit cell consists of a single particle surrounded by a bounding envelope. In the case of our two-dimensional array of particles, a cylindrical bounding envelope is employed which is perpendicular to the plane of the two-dimensional array and whose axis penetrates the centre of the cell particle. The bounding envelope represents the point where ^ is zero , i.e. the midpoint between the cell particle and an imaginary particle outside of the cell. By assuming symmetry about the axis of the envelope, i.e. that there is an imaginary particle at every angular displacement about the cell axis, the fully three-dimensional problem is simplified to a much more amenable axisymmetric one. There is a minimum radius for the bounding cylinder that corresponds to the closest packing allowed in the two-dimensional array. Assuming a close-packed structure where 31 Chapter 2. Theory WATER Figure 2.4: Model III: Unit Cell Representation of a Sphere at an 0/W Interface each spherical particle has 6 nearest neighbors this minimum radius is 1.0513a, where a is the radius of the particles in the array. The details of the cell model employed in this thesis are given below. Consider a lone sphere of radius a at a planar O/W interface. The sphere is centred on the axis of an axisymmetric, bounding cylinder. The nearest distance from the sphere surface to the cylinder wall is d and at this point ^ = 0. The immersion of the sphere at the interface is again defined by f3 = 7r — 8 (see Fig.2.4). Since this is an axisymmetric problem, equation (2.30) applies in the water phase and the Laplace form of equation (2.30) applies in the oil and solid phases. Again, note that the appropriate dielectric constant must be used when applying equation (2.29) in each of the three separate phases. The boundary conditions for this case are 1. $ —» 0 on the surfaces S wz and S ; as z —» ±oo. oz 32 Chapter 2. Theory 2. constant potential = \Po) on the particle surface in contact with the aqueous phase; S^. 3. ^ = 0 on the imaginary surface; S = d + a, (Sam + Swm). 4. zero charge ( | £ ^ f — f^^f = 0) on the oil/particle interface; S^. 5. zero charge (f^^f - = 0) on the oil/water interface; Sou,. There is no net force acting parallel to the interface in this model but there is still a force acting normal to the interface. The expressions for the reduced normal forces are similar to equations (2.40) and (2.41), i.e. for a zero charge boundary condition on S ow 2ir(z-+*+) / e J = { ^Z5. 2 w [ ( «W), Ox + r}rdr «.(»£)>) e Oz + (2.42) w and for a constant potential boundary condition on S^ = ar(»- + «+) / Z 5 . ( ^ . ) . e„, J { 2 oz + + 2 rdr r) oz (2.43) The reduced free energy of the double layer is given by equation (2.37) where the integration limits are changed according to the unit cell geometry, i.e. r has the limits J£ where S — a + d is the radius of the bounding envelope and z has the limits JT^,s When applying equation (2.37) in the non-polar phases the osmotic term vanishes and the sign of the electrostatic field energy term, W, must be changed [29]. 2.5 2.5.1 Finite Element M e t h o d Introductory Remarks The finite element method (FEM) is a general technique used to approximate solutions to complicated continuous systems. The origin of the FEM cannot be traced to any one person or period in time; rather, it has evolved from various approximation techniques Chapter 2. Theory 33 used by mathematicians and engineers. As would be expected from its varied origins, the FEM encompasses many approximation methods and allows for a great deal of choice when tackling a problem. A brief introduction to the philosophy behind the FEM is provided followed by a discussion of the techniques used for the analyses presented in this thesis. There are many comprehensive texts on the FEM available should the reader require more specifics; the text by Zienkiewicz [50] is particularly recommended. When posed with a continuum problem that has no tractable solution the investigator must resort to approximation methods to obtain a result. These approximation techniques usually involve discretizing the domain of the problem to model the continuous system.The FEM is such an approximation technique that involves subdividing the domain of the problem into individual subdomains or elements. The behaviour of the continuum (or field) variable is approximated in each element by a parameter function. In a thermal analysis this parameter function would represent the temperature function; in an electrostatic analysis the potential function. This process allows the overall behaviour of the system to be approximated by solving for the parameter function within each of the elements and then reassembling the elements to form the original domain. Each element is defined by a set of nodes that represent discrete locations in the domain. The value of the field variable at each node is the unknown quantity for which we are solving. The continuum within each element is approximated by the elemental parameter function which is, in effect, a function that interpolates between the nodal values of the field variable for that element. Each elemental parameter function, $ , is e represented as a sum of the products of the field variable, fa, and the so-called basis function, Ni, associated with each of the i nodes in that element, that is r = Y^Ni<fn (2.44) where M is the number of nodes associated with the given element. The basis functions Chapter 2. Theory 34 are in turn functions of the displacement co-ordinates of the element and are most often chosen such that equation (2.44) has the form of a polynomial. As will be revealed shortly, integrals of the products of the N and products of the spatial derivatives of the N often arise in the FEM. These integral expressions can be very complicated so it is advantageous to transform each element to a standard co-ordinate system such that a common numerical integration scheme, such as Gauss-Legendre quadrature, can be applied to each element. If the basis functions are expressed in terms of this local co-ordinate system then the resulting set of equations can be solved in a very efficient manner. When the integrals associated with each element are evaluated the effect of the transformation from the global domain to the standard elemental domain is accounted for by the usual methods used in calculus. This routine transformation procedure leads to one of the most popular features of the FEM: an ability to handle irregular element shapes and sizes in a routine manner. If the geometry of the element is denned by the same nodes that are used to represent the field variable then the element is said to be isoparametric. In this case any position within the element domain can also be expressed in terms of the basis functions, Ni, and the global co-ordinates of the nodes. Thus a position in the element may be expressed in terms of the local co-ordinate system and then readily transformed to the corresponding position in the global domain by using an expression similar to equation (2.44) where the (pi are replaced by the co-ordinates of the nodes in the global system. This has obvious advantages when integrating expressions that contain displacement terms. The premise used to obtain the approximation has not yet been discussed and it is appropriate to do so now. There are several approaches to take when establishing a criterion for the approximation; this should not come as a surprise since, as mentioned above, the FEM encompasses a general family of approximation techniques. The most common approach taken to solve problems described by second order differential equations is the Chapter 2. Theory 35 Galerkin method of weighted residuals. 2.5.2 Galerkin Method of Weighted Residuals The true behaviour of the field variable, <3>, is only being approximated by the parameter function, so when $ is substituted into a governing equation, = / , a residual Rn results R n = L{$)-f (2.45) where L is some general differential operator, / is some known function that represents the boundary conditions of the problem, and R is considered over the whole domain, fi. A weighted residual method involves minimizing RQ by integrating a biased residual expression over the domain of the problem. The bias is provided by introducing a set of independent weighting functions, Wi. Thus, the overall residual error is reduced by requiring that a sufficient number of integrals of the error over fi, each weighted independently, be zero [w R dn l n =0 ;Z = 1,2,... M ) (2.46) where it is convenient to use M functions; one associated with each node in the element. If the domain is represented by only one element then $ would be given by an expression of the form of equation (2.44). Substituting $ for $ in the residual term of equation (2.46) results in a set of M simultaneous equations and unknowns (<pi). The system of equations can be efficiently represented in matrix notation as K $ =f e (2.47) e where $ = (01,^2,^3, •• Am) (2.48) 36 Chapter 2. Theory Kf = / WiL(Ni)dn; 1 < l,i < M { ft = / Wtfdtt; Jo (2.49) 1< I< M (2.50) By selecting a suitable set of weighting functions and defining the boundary conditions, the unknowns can be solved to obtain the approximating function given by equation (2.44). If the domain is divided into J elements then there will J equations having the form of equation (2.47). If the boundary conditions between the elements can be specified, which is usually a trivial matter based on continuity requirements, then the complete set of unknowns can be obtained by assembling the individual K matrices and f vectors e e for each element into a global matrix K and vector / respectively. This assemblage step is a common feature of all FEMs. As a legacy to the FEM's beginnings in structural engineering K is customarily referred to as the stiffness matrix and / as the load vector. In the Galerkin method of weighted residuals the basis functions are used as the weighting functions. This choice leads to a significant computational advantage in that it frequently leads to symmetric stiffness matrices. 2.5.3 Choice of Basis Functions As noted above the basis functions, JV,-, are the functions used to interpolate the field variable between the nodes of the element. Naturally it is desirable to select a set of basis functions that best describes the behaviour of the continuum as a more accurate result will be obtained with fewer elements. For instance, if a one-dimensional continuum is described exactly by a quadratic relationship then the exact answer could be obtained by using a single element with second order basis functions but, if linear basis functions were being used, the domain would have to be successively discretized until the piecewise linear approximations sufficiently represented the curve of the quadratic. Generally, the Chapter 2. Theory 37 true behaviour of the continuum cannot be described simply and so one has to weigh the advantages and disadvantages of using a more complex set of basis functions. More complex parameter functions require more discrete variables per element to be solved so unless the parameter functions almost (or exactly) describe the 'true' behaviour of the solution within the element, there is no great saving to be made when employing complex parameter functions in an effort to reduce the number of discrete variables. The geometry of the elements has to also be considered. Recall from above that if the element is isoparametric, the same basis functions are used to transform the element between the global and local domains as are used in equation (2.44). If curved boundaries are to be represented then a set of basis functions of at least quadratic order will have to be used. The alternative is to approximate the curved boundary by a collection of linear boundaries but this approach is prone to error as the radius of curvature of the boundary decreases. Quadratic order basis functions were selected for the applications described in this thesis. Another property of the basis functions that has to be considered is their continuity between elements. This factor is best left until the weak form of the governing equation is discussed. 2.5.4 Weak Formulation via Green's Theorem The concept of the weak formulation is best explained by an example. Consider Laplace's equation for heat transfer in two-dimensions where T is the temperature (field variable) andfcis a thermal conductivity. Using the Galerkin method of weighted residuals we must evaluate the elemental stiffness matrices 38 Chapter 2. Theory which have the form =L 'l<*s> £<*f»** r NT + (252) where the integration is performed over the element area A and N is the set of basis e functions associated with the element, i.e. one basis function per node. The second order derivatives of the basis functions in equation (2.52) imply that the temperature as well as the the first derivative of the temperature must be continuous along the boundary of the element; this represents a C -continuous problem [50]. This 1 means that there must be slope continuity between all elements in our domain and that k must also be continuous. It is not possible to treat a composite material using the original form of equation (2.52) but we can allow for material discontinuities if we use the weak formulation. The weak formulation to a problem is an integral formulation that contains in its integrand derivatives of the field variable that are of a lower order than those in the original governing differential equation [51]. The weak form of equation (2.52) is obtained using Green's Theorem which is essentialy a multi-dimensional version of integration by parts [52]. Green's Theorem is related to the divergence theorem in that it connects an integral over a set of points in some n-dimensional space with another integral over the boundary of that set of points. If we consider the identity V • (N pT) = N V -p + V N -p T T (2.53) T then the two-dimensional form of Green's Theorem is / N V -pdA= I N p- ndC - [ V N • pdA T JA T T JC JA (2.54) where C is the line enclosing the area A and n is the unit outward normal vector from C. It is customary to perform the integration around C in a counter-clockwise direction. Chapter 2. Theory 39 Identifying kVT as p, equation (2.52) becomes (2.55) where n and n are the direction cosines between the normal vector and the co-ordinate x y axes. Equation (2.55) is the weak form of equation (2.51) and represents a C°-continuous problem. Putting the governing equation in this form removes a degree of smoothness from the problem since there is no longer a requirement for the first derivative terms to be continuous and material discontinuities are now permissable. This form of the equation is certainly more physically realistic than the original equation if it is to be applied to a discontinuous problem. The integral over the boundary C in equation (2.55) can be identified as the heat flux leaving the element. Rather than solve for the boundary flux we may specify it, that is each elemental boundary integral represents the function f in equation (2.50). e When assembling the load vector, thefluxterms associated with each node will be added together to get the net flux crossing the corresponding boundary. Unless there is some sort of heat source (or sink) associated with an element boundary then by continuity the heatfluxentering the boundary must equal the flux leaving the boundary. Thus all the boundary integrals for the nodes from the interior of the global domain, Q, will cancel each other out and the load vector will have zero values for the interior nodes. The flux terms associated with the global boundary are determined by the boundary conditions of the domain. If one boundary is subject to a heatfluxq then the result of the boundary integral (first term in the right-hand side of equation (2.55)) is added into the load vector elements corresponding to the nodes on that boundary; an example of this procedure is given in Appendix B. If another boundary has a constant temperature, T , boundary condition then T„ is added into the appropriate elements in the load vector e 40 Chapter 2. Theory and the stiffness matrix K is conditioned to ensure that the unknowns c6; associated with the constant temperature boundary equal T, when equation (2.47) is solved. 2.5.5 Application to Thesis Models This section presents the forms of equations (2.29) and (2.30) that were used in the'~ programs for determining the electrostatic potential distributions. It is convenient to consider a single elemental stiffness matrix and to separate it into two parts: K = K + K% e (2.56) e L where K\ is the contribution from the Laplacian operator and K B is the contribution from the charge density term in equations (2.29) and (2.30), i.e. the left-hand and right-hand sides of the equations. The contribution from a Laplacian operator in Cartesian co-ordinates was given in the last section, equation (2.55), so the (i,j) element of if £ is K (i,j) = e L e JJ—X-) + (— )( — )}dxdz (2.57) where e,, is the value of the dielectric constant in the element and A? is the area of the element. The weak form of the Laplacian operator in cylindrical polar co-ordinates contains a factor of r, the radial displacement co-ordinate. Recall that for isoparameteric elements, a co-ordinate within an element is given as a linear combination of the basis functions and the nodal co-ordinates. Thus the element of the cylindrical form of the Laplacian operator is =^ t e / ^ M ^ X ^ ) + i^§)(^f)¥rdi} (2.58) where the sum is over the M nodes in the element and r is radial co-ordinate of the k th k node in the global domain. The advantage in using isoparametric elements is apparent 41 Chapter 2. Theory in equation (2.58); there is economy in using the same basis functions to represent the field variable and the elemental co-ordinates. The non-linear charge term has to be linearized (i.e. a factor of xp has to be extracted) so that it may be accounted for in the system of linear equations (2.47). The exponential terms from the righthand side of equation (2.29) can be expanded to yield [exp(z_V) - exp(z V)] = + *+) + ^(zl-z ) 2 + + + ^(z _ 3 + *») + ••• (2.59) This may then be linearized by extracting a power of rp [exp(z_V>) - exp^+VO] = ip • A (2.60) where O + ^-(zi+zl) A = (*_ + z ) + + (2.61) + .- and where rp is the unknown quantity and rp is a guess value of rp. For the purposes of e this work the first twenty terms were retained in the expansion; this results in an error of <.1% for a maximum reduced potential of 10 for a 1:1 electrolyte. This problem requires iteratively solving for rpi until the change from the previous value rp\ is less than a selected tolerance 8, that is the solution is assumed to have converged when \tiZjtl\< • S (2.62) for all values of i. The initial guess potentials are obtained by solving the DH form of the PBE where only the first order term in the expansion is retained. In this case the K matrix element B is K (i,j) e B = X[ NiNjdxdy 2 (2.63) J Ae where A was defined in equation (2.9). On subsequent iterations we have 2 K (hJ) B = ^ - I lZ_ JAe NiNjdzdy (2.64) Chapter 2. The new Theory 42 guess value of xjj* is calculated using an under/over relaxation parameter ct #=# + aty,--#) (2-65) For the special case of a = 1 this amounts to using the calculated value from the previous iteration. Chapter 3 Methods This chapter contains details of the methods used to obtain the results presented in the next chapter. The finite element grid design used for the numerical analysis has gone through several evolutionary stages. Thus the reasons for the final design choice and the advantages and disadvantages of this design are discussed. Sensisitivity tests were performed to determine solution dependence on element densities and domain sizes. The force and free energy calculations involving the electric field can be strongly dependent on the grid density so special consideration is given to these calculations. The elements of the stiffness matrix, the forces andfreeenergies are all calculated via numerical integration. Thus, the integration method and the number of integration points chosen for performing the various numerical integrations are discussed. Thefiniteelement analysis yielded very large, sparse sets of simultaneous equations so the method used to solve them is also discussed. Finally, an overview of the computer program is presented. Unless notedc otherwise, all calculations presented in "this chapter were performed with e = 4, e = 2, p D and e — 80. w 3.1 Finite Element Grid As described in the previous chapter, the new models considered in this thesis are twodimensional and involve a particle with a circular cross-section embedded in a planar interface. To characterize the dependence of the electric double layer interaction on the system parameters (e.g. particle size, contact angle, etc.) a series of test cases had to be 43 Chapter 3. Methods 44 evaluated. Rather than design a separate grid for each test case, it was desirable to have a general routine that created a grid for any combination of independent variables. The geometric parameters that must be variable are: • the contact angle (hence the immersion depth) of the particle, • the separation between particles (the closest distance from the particle surface to the median plane), • the particle size, • the distance of the surface S (representing infinite distance) from the interacting T particles, and • the element density anywhere in the grid. While separate grids were used for each of the three models considered in this thesis, they all had the same basic design (see Fig.3.1). The particle is always represented by a circle that is divided into a number of equal inner 'pie' sections and outer 'spoked' sections. In the case of the cylinder model, the area surrounding the particle is divided into eight sections. There are four sections in each phase with two sections on the front side (side closer to the median plane) and two on the back (side furthest from the median plane). In Model I there is axial symmetry about the line of centres joining the particles so only four of the eight sections surrounding the particle need to be mapped. If there is a constant potential boundary condition on the particle surface, the interior of the particle may be ignored. If there is a constant charge boundary condition on the particle surface then half of the particle is mapped since the potential distribution inside the particle is axisymmetric. In Model III there is symmetry about an axis that passes through the centre of the particle and is perpendicular to the interface so only half the area needs to be mapped, i.e. four sections plus half the particle. 45 Chapter 3. Methods Model II Model III Figure 3.1: Typical Grids for the Three Models Chapter 3. Methods 46 In the present scheme the particle may only be sectioned into equal fractions of TT and thus the contact angle is limited to those fractions. The radius of the inner 'pie' region is variable but it was found that it didn't affect the answer unless it was much smaller than or almost as great as the particle radius; it was customarily set to half the particle radius. The sections surrounding the particle are all quadrilateral so they may be treated by a common set of subroutines to calculate the nodal co-ordinates, spacings, etc. When using a standard quadrilateral grid one usually varies the number of elements in the X and Y directions to see how the solution varies. For the present grid design one therefore has to consider variations in all eight sections plus the particle. Fortunately, however, a basic requirement of the FEM is that whenever any two sections join, the nodal locations must match. This places restrictions on the grid, for instance each of the four sections on the back side of the particle in Model II must have the same number of elements in the X direction. This matching restriction reduces the number of variables that have to be tested. The node spacing in the X and Y directions is also variable in each of the eight sections; these spacings are also subject to the matching requirement. A spacing capablility was incorporated so that more elements could be packed around the particle where the potential varies most strongly. The major disadvantage of this design is that there isn't a smooth distribution of elements surrounding the whole particle; the elements located at the top and bottom of the particle are much more distorted than those around the middle portion. This isn't very important when solving for the potential but it is a factor when calculating the electric field around the particle. Even with the above restrictions there is still a lot of choice when selecting the grid parameters; this is advantageous since it allows a good deal of flexibilty but disadvantageous because it makesrigoroussensitivity analyses impractical. The procedures used to 47 Chapter 3. Methods evaluate solution sensitivities to grid variations are discussed below. 3.2 Sensitivity Tests The general strategy employed for the sensistivity tests was to use a dense grid to calculate a solution and to then decrease the density of the grid in certain regions to see what effect this had on the solution. The advantage of using this strategy is that a dominant parameter is readily identified. An obvious disadvantage is that the computing costs are substantial. The upper limits on the number of nodes and elements allowed in the present scheme is 5000 and 2500 respectively. In most cases it was found that the solution had converged (i.e. maintain a constant value for further increases in grid density) for node and element numbers of ~ 2500 and ~ 1500 respectively. The reader will notice that the bounds of the finite element grids are finite whereas the boundary conditions call for an infinite boundary. Rather than map the elements to infinity, which introduces added complexities [50], it was decided to set the infinite boundary far enough away from the particle that the potential could be assumed to be zero. This is straightforward for the water phase if the Debye-Hiickel expression for the potential distribution at a plane wall is considered. In this case, $(z) = ^exp ** - (3.66) where ty(x) is the potential at a distance x from the plane wall which has a surface potential of \& . Equation(3.66) predicts the potential will have decreased by a factor of ~ 10 -7 0 at a distance of 15 double layers ( « ) from the surface. For surface potentials greater _ 1 than the DH limit and/or curved surfaces the potential drops even more rapidly. Thus setting the water phase boundaries at 15/c from the particle surface is a conservative -1 estimate for the zero potential assumption. The potential drops off much more slowly in the oil phase than in the water phase because there is no screening effect due to the Chapter 3. Methods 48 presence of ions. The 'infinite' boundary in the oil phase was determined by performing a series of sensitivity tests to see how the interaction quantities were affected. These tests will be discussed below. As shown in the last chapter, the force and free energy expressions depend directly on the electrostatic potential through the osmotic term and also on the electricfield,which is expressed in terms of spatial derivatives of the electrostatic potential. The potential distribution is obtained directly by the FEM but since this problem must be treated as C°-continuous, the potential derivatives can only be approximated from the results. Thus, under conditions where the local potential undergoes rapid changes with position, approximating derivatives will lead to large errors which are sensitive to grid parameters. It is best to examine separately the effects that grid parameters have on each of the three calculated quantities; the force parallel to the 0/W interface acting between two particles, the force normal to the O/W interface acting on the particles and the free energy calculations. Each of these quantities involve differentiating and integrating over different regions of the grid and thus exhibit different sensitivities. 3.2.1 Parallel Force Between Two Particles The force acting parallel to thefineof centres joining two particles was calculated using the left-hand side of equation (2.22), i.e. by integrating the electrostatic and osmotic stresses over the median plane between the two particles. The electrostatic stress term only involves the potential derivative normal to the fine of centres. It turns out that the effect of the double layer interaction is to cause the electric field to vary relatively slowly in the region of overlap. The potential distribution for a typical case for Model I is shown (as a set of equipotential contourfines)in Fig. 3.2. The effect of the double layer interaction is starkly obvious. On the back side of the sphere, the potential drops off very rapidly due to the screening effect of the ions. On the front side, the overlapping Chapter 3. Methods 49 O in Figure 3.2: Electrostatic Potential Distribution for One of Two Interacting Spheres in an Electrolyte double layers serve to moderate the potential distribution. The influential grid parameters for calculating the parallel force are: 1. the number of elements in the overlap region, 2. the extent of the median plane required to approximate an infinite boundary, 3. the number of elements extending out along the median plane, and 4. the number of integration points per element along the median plane. The number of elements in the overlap region is determined by the number of radial elements in the particle (NEA) and the number of elements between the particle and the median plane (NED). The dependence of the parallel force on each of these parameters is shown in Fig. 3.3 for a few representative cases of Model I. The behaviour of the solution as a function of NEA is plotted for several cases of tta with I{J = 2 and K[R— 2a] = 1. Note 0 50 Chapter 3. Methods that the effect of increasing NEA is insignificant above NEA=8 and that cases with larger na are more sensitive to changing NEA than cases with smaller na. The dependence of the solution on NED is plotted for Model I for ip — 2 and na = 10 at various separations. 0 Note that smaller separations require fewer elements for convergence. For appreciable double layer overlap, the figures indicate that the parallel force calculation is insensitive to changes in NEA and NED above the coarsest settings. This observation is not surprising if the contributions to the integral over the median plane are examined. When significant double layer overlap occurs (i.e. K[R — 2a] < 10) the osmotic term contributes ~ 95% of the total force with the electrostatic stress term accounting for the remainder. For two cylinders at an O/W interface with significant double layer overlap, the electrostatic stress contribution from the oil phase accounts for <1% of the total force. Since the osmotic term accounts for the majority of the parallel force in cases of significant double layer overlap it is not surprising that the parallel force is so insensitive to grid changes for these cases. The reason for this is that the osmotic term is calculated directly from the electrostatic potential which is in turn being directly solved by the FEM. The potential distribution is not very sensitive to changes in the grid density above a moderate setting. When double layer overlap is minimal (i.e. n[R — 2a] > 10 or contact angles ^> 90° for cases with large KO) the parallel force resulting from the electrostatic stress in the oil phase dominates. The effect of NEA and NED is not much different from that shown above because the potential varies relatively slowly with displacement in the oil phase. The dominant parameters for this circumstance are the infinite boundary dimensions in the oil phase. Thus sensitivity tests were performed to determine the effect of setting the 'infinite' boundaries at ZOI in the oil phase (i.e. the depth of the oil phase) and at XI in the X-direction on the backside of the cylinders. Both of these quantities are measured from the particle surface to the appropriate infinite boundary. The results of these tests for cylinders with na=10, 9 = 135°, n[R — 2a] = 5, rp = 2, O/W interface a 51 Chapter 3. Methods % Difference from Converged Result NEA NED Figure 3.3: Effect of NEA and NED on the Parallel Force uncharged) are shown in Fig.3.4. In the ZOI tests XI was held constant at IOOAC and -1 in the XI tests ZOI was held constant at 100/c . -1 Note that the solution has essentially converged on the final result in both tests. For the purposes of this work, XI and ZOI were fixed at 80/c when calculating longer range -1 forces (K[R — 2a] > 10). For short range calculations XI was set to 20K . _ 1 The number of elements extending out along the median plane will influence the parallel force calculation differently in the water and oil phases. For the aqueous phase, an examination of the contributions to the force integrand reveals that >99% of the contributions come from the region below the height of the particle above the O/W interface. The greatest double layer overlap occurs within this region (see Fig.3.2). The number of elements extending out along the median plane in this region is determined by NEA whose effect has already been addressed. As expected, varying the element density in the region above the particle had a minimal effect on the final result. In the oil phase the interaction results only from the electrostatic stress which has its maximum close Chapter 3. Methods -*H 0 i 25 52 1 1 1 50 75 100 ~i 40 f - 70 /cXI 1 100 1 130 1 160 *ZOi Figure 3.4: Effect of ZOI and XI on the Parallel Force for Minimal Double Layer Overlap to the interface and then falls off slowly, asymptotically approaching zero towards the 'infinite' bound. This smooth behaviour allows the electrostatic force term to converge rapidly as a function of the number of elements extending out along the median plane. The number of integration points per element necessary for a satisfactory result are discussed in Section 3.6. 3.2.2 N o r m a l Force A c t i n g on a Particle at an Interface As shown in Appendix A, there is also a force acting normal to the interface on a particle at an O/W interface. This force is calculated by evaluating, over the interface, surface integrals which involve osmotic and electrostatic stress terms. The normal force was only calculated for particles on an uncharged interface for reasons given in Chapter 4. Before discussing the details of the sensitivity tests a few aspects regarding the normal force should be clarified. When discussing the normal force below, a distinction is made between the 'front' side and 'back' side normal forces. The 'front' side force refers to 53 Chapter 3. Methods the force calculated by integrating over the 0/W interface between the clyinder surfaces and the median plane separating the pair of cylinders. The 'back' side force refers to the force calculated by integrating over the O/W interface from the cylinder surface to the 'infinite' boundary. The sum of these two forces yields the total normal force acting on the particle. The normal force still exists at infinite particle separation and in this case the 'front' and 'back' side forces are equal. In general, the effect of the double layer interaction is to reduce the 'front' side force. In the case of an uncharged interface, potential derivatives normal and parallel to the surface and the osmotic pressure at the interface are integrated. Figure 3.5 is an equipotential contour plot for one of a pair of cylinders under the conditions: na = 10,rpo — 2,K[R — 2a] — 5,6 = 90°, on an uncharged interface. The double layer overlap interaction dominates at this separation and thus the electric field varies slowly on the 'front' side of the cylinder. The osmotic contribution accounts for ~ 85%, the (^) term 2 for ~ 15% and the (f^) term for ~ 1% of the total front side normal force when double 2 layer overlap is significant. Because the major contribution isfromthe osmotic term and thefieldvaries slowlyfromthe particle to the median plane, thefrontside normal force is not very sensitive to changes in grid parameters when significant overlap occurs. On the 'back' side of the cylinder, the potential varies much more sharply and the (f^) term assumes a greater contribution; osmotic ~ 50%, (^) 2 2 ~ 50%, (f^) ~ 1% 2 for 6 = 90°. The node spacing in the X-direction on the back side of the cylinder is determined by three variables • FCXI which is a ratio of the length of the elements furthest awayfromthe cylinder to the length of the elements nearest the cylinder, • XI which is the distancefromthe cylinder surface to the infinite boundary, and Chapter 3. Methods 54 Figure 3.5: Potential Distribution for One of Two Parallel, Interacting Cylinders on an Uncharged 0/W Interface 55 Chapter 3. Methods • NEI which is the number of elements, in the X-direction, extending from the cylinder surface to the infinite boundary. For short range (K[R — 2a] < 10) interactions, XI is always set to 20 and because of memory limitations NEI is always set to 15. Thus the node spacing, in the X-direction, on the back side of the cylinder was varied using FCXI to determine the effect the spacing had on the normal force. It was found that increasing FCXI from 1 had a minimal effect on the back side normal force when 6 = 90°. When 6 ^ 90°, however, there is a noticeable effect which is shown in Fig.3.6 where the back side force is plotted as a function of FCXI for cylinders with «a = 10 and = 2 for the maximum arid minimum contact angles considered in this thesis. Both graphs show that the force has essentially reached a constant value for FCXI = 12. FCXI was set to 12 for cases where 6 ^ 90°. On the front side of the cylinder, the node spacing in the X-direction is determined by FCXD. As the separation between the cylinders increases, the double layer overlap decreases and the electric field will vary more rapidly in the X-direction (i.e., it will resemble the back side) and thus the calculated normal force on the front side of the cylinder becomes more sensitive to the node spacing. When there is appreciable double layer overlap, varying FCXD does not affect the calculated forces (parallel or normal) so FCXD was set to 12 for cases with 0 f 90°. When there is no surface charge on the interface the normal force can be calculated using the oil side or water side (|^) term (see Appendix A). This provides a useful check 2 as to whether the normal force is being accurately calculated. When 6 = 90° the (f^) 2 term accounts for ~ 1% of the total force and so the check is rather meaningless. For contact angles other than 90° the normal derivative term becomes more important as the electric fieldfinesare no longer parallel to the interface. In general, the normal forces 56 Chapter 3. Methods o 6 = 157.5 0} O 2 O o D E -B O S B B B Q 1 1 1 10 15 20 1 10 FCXI —T— 15 FCXI 1 20 Figure 3.6: Back Side Normal Force as a Function of FCXI calculated using the oil and water side derivatives showed poor agreement so a series of tests was performed to determine which of the two results was more reliable. The dominant parameter affecting the (^) 2 term is NEA, which determines the number of elements distributed radially about the cylinder. Recall that the radial spacing of the nodes must be in equal fractions of TT (by design) and so the height of the elements on the 0/W interface adjacent to the cjdinder will be determined by NEA. The back side normal forces, calculated for various values of NEA using the oil and water side derivative terms, are shown in see Fig.3.7. The forces were calculated for contact angles of 45° and 135°. For both cases the force calculated using the oil side derivative term appears to equal the true asymptotic result whereas the water side result is more sensitive to NEA. This observation is not surprising since the potential varies much more slowly in the oil phase than in the water phase. The fact that the water side result is more stable for the 45° case than the 135° case may be explained as follows; at 45° the double layer is projected onto the 0/W interface, which has a moderating effect on the electric field, Chapter 3. Methods 57 Figure 3.7: Back Side Normal Force as a Function of NEA but at 135° the double layer is projected away from the 0/W interface. The net effect is that at 45° the potential varies more slowly, in the X-direction, along the O/W interface than at 135°. 3.2.3 Free Energy The free energy associated with the electric double layer and electrostatic stress in the oil phase is calculated by integrating the square of the electric field over the whole system and the osmotic pressure over the water phase. As will be shown in the next chapter this calculation yielded poor results when compared with the published results for Model I. Since the electric field is integrated over the whole domain, the calculated free energy will be sensitive to grid densities over the whole domain, especially in the water phase in the region close to the particle surface where the potential gradients are very steep. The parameters governing the element density in the region near the particle surface are: 1. the number of elements distributed radially about the particle, 58 Chapter 3. Methods 2. the number of elements between the particle and the 'infinite' boundary, 3. the spacing of the nodes close to the particle, and 4. the extent of the infinite boundaries. As mentioned above, the number of elements distributed radially about the particle is determined by the grid variable NEA while the node spacing (in the X-direction) between the particle surface and the infinite boundary is determined by the grid variable FCXI. The reduced free energy calculated for Model I with na = 5,ip = 2, and at 0 'infinite' separation as a function of NEA and FCXI is shown in Fig. 3.8 (for the NEA tests FCXI was fixed at 10, for the FCXI tests NEA was fixed at 32). Note that the solution is more sensitive to the FCXI parameter than NEA; this is reasonable since the potential gradients in the angular direction about the particle are very small compared to the gradients normal to the particle surface. Neither of the solutions in Fig. 3.8 have converged, in fact they don't even show a smooth asymptotic approach to a final value. These calculations were repeated several times, varying the convergence tolerance (6), and the same behaviour was reproduced each time. This behaviour clearly indicates the inability of the present grid design to accomodate the free energy calculations. It should be noted, however, that the overall trend when increasing the grid density is that the free energy value decreases. The number of elements in the Y-direction in the regions above and below the particle will also affect the number of nodes near the particle surface. The number of elements in the Y-direction in the region above the particle in the water phase is determined by the grid parameter NWA. The result of varying NWA is shown in Fig 3.9 for Model I with na = 5,^0 = 2, and at 'infinite' separation. From the figure, it is apparent that the solution has not yet converged for the range of NWA shown, and again the solution doesn't show a smooth asymptotic approach to a final answer. NWA=25 is the maximum 59 Chapter 3. Methods d "H 12 —. 21 i 30 1 39 —:n 48 ci "M 0 , 5 : NEA r- 10 1 1 15 20 FCXI Figure 3.8: Free Energy of a Sphere in an Electrolyte as a function of NEA and FCXI value allowed due to memory constraints (not to mention costs) while keeping all other grid parameters at sufficient settings. The extent of the 'infinite' boundaries will also affect the free energy calculation since the integration is supposed to be performed over an infinite space. In the water phase the potential falls off rapidly so there is virtually no change in the free energy once the 'infinite' boundaries in the water phase are more than 15K away from the particle -1 surface. Model III, however, contains a non-polar oil phase and hence the potential drops off much more slowly than in the water phase. Figure 3.10 shows the reduced free energy calculated as a function of the 'infinite' boundary in the oil phase, ZOI, for Model III with na = 1 0 , = 2,6 = 90°. As can be seen in thefigure,the free energy has approached the asymptotic result by ZOI=40/e . An examination of the contributions to the total -1 free energy of the unit cell reveals that the oil phase accounts for less than 1% of the total. This is explains why the total free energy is relatively insensitive to the 'infinite' boundary in the oil phase. Chapter 3. Methods SI Figure 3.9: Free Energ)' of a Sphere in an Electrolyte as a function of NWA *1 <N Figure 3.10: Free Energy of a Sphere in a Unit Cell as a function of ZOI 61 Chapter 3. Methods d 0) o l- o o , I.I.I.UIIII mill . . I.I.IJJIIII i.I.I.IJIIII . l.i.i.uml • i.i.t.uiiit . 1 . 1 . 1 . U I 1 1 I . I.I.UJIIII 10"° 10" 10" 10" 10" 10" 10" 10"' 10° 7 6 5 4 3 2 Figure 3.11: Converged Value of/* as Function of 8 3.3 Convergence Tolerance and Relaxation Parameter The non-linear term in the PBE forces us to solve for the potential distribution iteratively. The solution is assumed to be converged when all nodal values of ip change by less than the tolerance, 8, from one iteration to the next. The definition of 8 was given in equation (2.62). A test was performed to evaluate the rate of convergence as a function of 8 for a Model I case: na = 10, K[R — 2a] = l ^ o = 2. Fig. 3.11 clearly shows that the solution attains the true asymptotic result by 8 = 10 ; the force values were equal for -4 five significant digits for values of 8 < 10 . -3 All results presented in this thesis were obtained with a tolerance of 10~ . All of the comparable data in the literature are either 4 reported to this accuracy or less. The values of tp are altered from one iteration to the next using equation (2.65). The relaxation parameter, a, was incorporated because it was found that under certain conditions the solution became unstable and, without under-relaxation, could not converge. Therefore, for constant potential surfaces with ip > 4, a had to be reduced below 1. In 0 62 Chapter 3. Methods general a was reduced by .025 for each unit increase i n ipo above 4 (the highest value of ipo tested was 10). W h e n two surfaces with equal, constant surface charges are brought together, the surface potentials must increase (for a positive charge). Thus it was found that a had to be decreased from 1 as the separation between the two particles decreased. For very small separations (K[R — 2a] =.25), a had to decreased to ~ . 6 , depending on the surface charge. 3.4 Orders of Integration Used T h e contributions to the stiffness matrices discussed i n the last chapter are such that numerical integration is required for them to be evaluated practically. The boundary integrals used to calculate the force expressions and the area integrals used to calculate the double layer free energy also require numerical integration. T h e additional error which is introduced into the approximation through numerical integration is a function of the distribution of the integration points which is i n turn a function of the order of the integration used and the number of elements over which the integration is performed. Accuracy and computing costs have to be considered when the order of integration for evaluating the stiffness matrix elements is chosen. T h e first inclination is to sacrifice computer costs for the sake of a higher integration order to increase the accuracy of the solution. It has been shown, however, that using higher integration orders may sometimes decrease the accuracy of the solution because the fortuitous 'cancellation of errors' inherent to numerical integration may be diminished [50]. Zienkiewicz [50] states that for a C°-continuous problem where quadratic order triangular elements are being used, a three point quadrature scheme should be adequate. W h e n the three point scheme was used with the axisymmetric models, anomolous solution behaviour occurred in the 63 Chapter 3. Methods region of the median plane, i.e. the potential oscillated slightly. When a four point integration scheme was employed this behaviour vanished and, in the case of Model I, the force calculations yielded results which showed excellent agreement with published data. The four point scheme was also used for Model II as slightly better results were obtained (~ .2% difference when limiting cases were compared with the Derjaguin approximation) compared with the three point scheme. The sampling points and weights for quadratic, triangular elements are found in most finite element texts [50, 51]. The force calculations discussed in the last chapter involve boundary integrals along the sides of the elements on the global boundaries. These integrals were performed over each element using Gauss quadrature by transforming the basis functions from the standard 0—»1 domain (area co-ordinates) to a -1—>1 domain. The results from each element were summed to give the total force. A number of different integration orders were evaluated for both the parallel and normal force calculations. It was found that a three point scheme was sufficient. The free energy calculations involve an area integral (or a volume integral for axisymmetric models) over the domain. Three, four and five point integration schemes were tested for Model I. A moderate improvement was realized by using the four point over the three point scheme but no further improvement was obtained with the five point scheme. The four point scheme was used for all the results in this thesis. 3.5 Solution of Equations The order of the stiffness matrix is determined by the number of nodes in the global domain. If the domain is regularly shaped the nodes may be numbered in such a fashion as to minimize the bandwidth of the stiffness matrix. Unfortunately, for the irregular geometries considered in this thesis, the resulting bandwidths were too large to use equation Chapter 3. Methods 64 solvers specifically designed for banded matrices. SPARSPAK [53] is a collection of subroutines for solving a system of linear equations where the coefficient matrix is large and sparse. This package uses Gaussian elimination to solve the system of equations but takes advantage of the sparseness of the coefficient matrix by reordering it to reduce the number of zeroes. There are five separate ordering choices available within SPARSPAK. Both the oneway dissection and the nested dissection methods are ideally suited for finite element applications. For non-linear problems where the structure, and hence the ordering, of the coefficient matrix remains constant but the numerical values do not, the nested dissection method realizes lower computer costs [54]. The nested dissection method for symmetric coefficient matrices was chosen for this application. 3.6 Software Overview This section gives a brief overview of the structure of the programs used to solve the governing equations by the FEM. A complete copy of the program used for analyzing Model II is provided in Appendix C. Three separate programs were written, one for each model, but they all have the same basic design. (While it would have been ideal to write one general program that could be used for all models the additional effort was deemed unwarranted.) The general structure of the software is outlined in the following subsections. 3.6.1 Data Input All the model variables are supplied via an input data file. An example of this file is given in Appendix C. This file contains the dimensions of the grid, the number of elements in the various regions of the grid, boundary condition values and flags (e.g., for constant Chapter 3. Methods 65 charge or constant potential), and the non-linear convergence tolerance and iteration limit. The same input file can be used for all three models although some quantities are meaningless for different models (e.g., the oil phase parameters aren't used in Model I). The grid dimensions are supplied in units of « _ 1 so the first task performed is to scale them by a factor of ^/e^; the reason for doing so was explained in the previous chapter. A few preliminary calculations of grid variables that are required by the grid construction subroutine are then performed. 3.6.2 Grid Construction There are two distinct stages to the grid construction • assigning the global co-ordinates to the nodes of the grid, and • assigning the appropriate nodes to the elements in the grid. As discussed earlier in this chapter, the grid design is such that the particle is divided up into 'pie-shaped' elements having equal fractions of 7r while the remaining regions are all quadrilateral. Assigning co-ordinates to the nodes in the particle is trivial. The node co-ordinates for the separate quadrilateral regions are all calculated by a common set of subroutines. All that needs to be supplied to this set of subroutines are the global co-ordinates of the boundary nodes of the regions and the spacing factors for the nodes in the X and Y directions. The node numbering is performed as the global co-ordinates are assigned. An essential part of the grid construction is assigning the global node numbers to the elements in the grid. As shown in Appendix B, each quadratic-order triangular element has six nodes associated with it which are numbered locally from 1 to 6. This local numbering scheme is used when calculating the contributions to the 6x6 elemental stiffness matrix. When the individual elemental stiffness matrices are assembled into Chapter 3. Methods 66 the MxM globed stiffness matrix the global node numbers (which range from 1 to M) are used to address the appropriate locations in the global stiffness matrix. The global node number assignments of the nodes in each element are stored in a K E Y array which has dimensions NELx6 where NEL is the number of elements. Once all the nodes have been assigned their global numbers (and co-ordinates) the KEY array is filled in the subroutine FILKEY. 3.6.3 Setting up the Equations Now that the nodes are numbered and assigned to the elements, the global stiffness matrix may be constructed for the first iteration. Before calculating the matrix elements it is necessary to perform a few SPARSPAK preliminaries. SPARSPAK requires the non-zero structure of the stiffness matrix so that it may reorder it. This structure may be input as a collection of submatrices (our elemental matrices) which are readily available using the K E Y array. The ordering choice must be explicitly made by calling the appropriate subroutine. The elemental submatrices are calculated in the subroutine SETMAT. The elemental stiffness matrix elements are calculated using Gauss quadrature. Efficiency is maximized by calculating the contribution for the same quadrature point for every element in the grid, before moving on to the next quadrature point. Since this is a symmetric problem, only the lower left-hand triangle of each elemental matrix is calculated and stored. Recall that the elements are transformed from the global domain to a normalized domain so that the integration may be performed efficiently. The spatial derivative terms and the integration intervals must account for this transformation; this is accomplished using the standard procedures in calculus that require evaluation of the Jacobian matrix [52]. When integrating numerically, the Jacobian must be evaluated at each integration point; this is performed in the subroutine JACOB. There are limits to the degree of distortion 67 Chapter 3. Methods an element can be subject to; the general rule [50] is that if the sign of the Jacobian changes at any point within an element then that element is too distorted. Too much distortion indicates that one-to-one mapping is no longer possible between the global and local domains, i.e. a point in the local domain may represent more than one point in the global domain. Thus, the sign of the Jacobian is checked from point to point (in JACOB) during the numerical integration on the first iteration. The elemental submatrices are then input into the SPARSPAK routine one by one. SPARSPAK automatically adds multiple contributions to the same global matrix element. The boundary conditions are now ready to be incorporated into the system of equations. This step involves assigning values to the right-hand side or load vector and conditioning the global stiffness matrix. For a constant potential boundary condition the boundary value, multiplyed by a very large number (10 in this program), is assigned to 25 the load vector elements representing the boundary nodes. The stiffness matrix is then conditioned by adding that same large number to the appropriate diagonal elements so that when the system of equations is solved, the boundary nodes will have the assigned boundary values. For a constant charge boundary condition the result of the boundary flux integral term is simply added to the load vector elements representing the boundary nodes at constant charge. The boundary flux integral term was discussed in the last chapter and a worked example of how to calculate it is provided in Appendix B. 3.6.4 Iterative Convergence The set of equations is solved iteratively until all the ip meet the tolerance requirement discussed above. Since the structure of the matrix is not changing it is not necessary to reorder the stiffness matrix from iteration to iteration. Only the contributions from the elements in the water phase need to be recalculated after the first iteration because the Chapter 3. Methods 68 non-linear term is only present in that phase. The whole matrix does have to be reloaded into SPARSPAK however, since the stiffness matrix is automatically zeroed after each solving step. 3.6.5 Force and Free Energy Calculations Once the solution has converged, the force and free energy calculations are performed. Boundary elements are collected for the appropriate force calculations and are passed to the subroutine BNDINT. This subroutine is used to calculate both the parallel and normal forces (where applicable). The force expressions are numerically integrated for each element on the boundary and summed to give the total force. The free energy calculation involves integrating over the whole grid and is identical to the integration procedure in SETMAT. The free energy calculations were not performed for the cylinder model (because they are so inaccurate) so the ENERGY subroutine used for Models I and III is also included in Appendix C. Chapter 4 Results and Discussion This chapter presents the numerical results obtained for the three models described in Chapter 2. The results for each model are presented and discussed separately. Before examining the results in detail, it is beneficial to review what order of magnitude the reduced quantities used in this thesis represent in practical situations. All distances in this chapter are referred to in units of K ; this is the standard 'yardstick' _ 1 used in double layer theory. Recall that / c _1 is the so-called double layer thickness and it represents the distance required for the potential to drop by an exponential factor at low potentials. For a 1:1 electrolyte concentration of 10~ molar, at room temperature; 4 rc - 1 ~ 3 X 10~ am. The double layer cases considered in this thesis range from na = 1 2 to 20; the experimental cases studied by Levine et. al. [1] had na % 10. The reduced potential, ip = ^ , represents the ratio of the electrostatic energy to thermal energy at a point in the electric field. A potential of ~ 25 mV at room temperature is equivalent to a unit of the reduced potential. The reduced force, / * , is measured in units of ^ Q * T • For a 1:1 electrolyte at room temperature this corresponds to a unit force of approximately 2 x l 0 - 1 3 i V . For a silica sphere with a radius of .3/im this unit force is approximately 80 times the gravitational force (gravity alone, no bouyancy effects) acting on that particle! It is common practice for workers in double layer theory to report a scaled reduced force for comparison purposes. The reduced force is usually divided by ip ^ and a factor which arises in the 2 Derjaguin approximation. For spheres, the scaled reduced force is given by 69 Chapter 4. Results and Discussion 70 oo and for cylinders r = ~ (4.2) T O O The force expressions for cyhnders are given as per unit length of the cylinder, where a unit length is / c . -1 The reduced energy, G", is measured in units of ^ 5 5 ^ • For a 1:1 electrolyte concentration of 1 0 -4 imately 7 x 1 0 molar at room temperature this corresponds to a unit energy of approx-21 J which is ~ 2kT. Again, a scaled reduced energy is often employed which for spheres is given as 47r/ca0 , 2 X) ^ ^ Free energy calculations were not performed for cyhnders for reasons given later in the text. Unless noted otherwise, all calculations presented in this chapter were performed with e = 4, e = 2, and e = 80. p 4.1 0 w Two Identical Spheres Fully Immersed in an Electrolyte This model was evaluated for two purposes: • to test the performance of this F E M analysis against previous work, and • to examine the validity of the common approximation employed when considering constant charge interactions between particles in a single electrolyte, i.e. assuming that | ^ is zero on the particle side of the particle/electrolyte interface. 71 Chapter 4. Results and Discussion The bulk of the previous numerical solutions for this model only consider the case where there is a constant potential on the spheres. Only McQuarrie et. cd. [42] have presented numerical results for the constant charge case by directly solving the PBE. Other numerical solutions based on the Derjaguin [37, 38, 40] and superposition [35, 38] approximations have been presented for the constant charge case. It is convenient to consider constant potential and constant charge interactions in separate sections. 4.1.1 Constant Potential Interactions As was detailed in Chapter 2 several numerical solutions to this model have been presented previously, the most comprehensive being those of Hoskin and Levine (HL) [28] and McCartney and Levine (ML) [34]. Both HL and ML report the force of repulsion and the free energy of interaction for two identical spheres, for a variety of sphere sizes, separations and constant surface potentials. The interaction force results calculated by the present FEM are listed in Table 4.1 along with selected results from ML and HL. Where available the ML results are shown as they are the more accurate of the two. The greatest difference between the FEM and ML values is 1.3 % (at larger separations, i.e. K[R — 2a] = 3) and in general the difference is less than .5 %. For the ip = 6 cases, only the HL results were available and 0 the comparison with the FEM results is poorer. This is not surprising since the HL data were obtained with a relatively coarse grid (144 points). ML reported 2-5% differences between their results and HL's; in nearly all cases the results of ML were lower than those of HL. This same trend is evident when comparing the FEM results with those of HL. The FEM results, however, are all marginally larger than ML's results so they may be slightly less accurate. McQuarrie et. cd. [42] considered constant potential cases with a much finer finite difference grid than ML and reported agreement to within 1%. Unfortunately McQuarrie did not publish any numerical values nor did he indicate if his 72 Chapter 4. Results and Discussion Table 4.1: Force of Repulsion (/"*) Between Two Spherical Particles with Equal Constant Surface Potentials in a Single Electrolyte V'o na method K{R 0.5 - 2a) 3.0 1.5 2 FEM 0.7772 0.3136 0.0691 5 ML 0.7737 0.3109 0.0681 5 FEM 10 0.7667 0.3168 0.0740 10 ML 0.7634 0.3142 0.0730 FEM 0.7632 0.3178 0.0757 15 ML 15 0.7601 0.3152 0.0747 4 FEM 0.7332 0.2286 0.0484 5 ML 5 0.7287 0.2264 0.0477 FEM 10 0.7228 0.2278 0.0504 ML 10 0.7185 0.2257 0.0497 FEM 15 0.7194 0.2274 0.0503 ML 15 0.7150 0.2253 0.0503 6 FEM 0.5992 0.1517 0.0312 5 HL 5 0.6098 0.1583 0.0343 FEM 0.5929 0.1500 0.0319 10 FEM 15 0.5936 0.1496 0.0322 HL 15 0.5939 0.1557 0.0353 FEM : Present Finite Element Method. ML : McCartney and Levine [34]. HL : Hoskin and Levine [28]. Chapter 4. Results and Discussion 73 results were higher or lower than those of ML. In any case, Table 4.1 indicates that the FEM employed in this work shows good agreement with the previously published results for the interaction force between two spheres with constant equal surface potentials in a single electrolyte. The electric double layer free energy may also be directly calculated using the potential distribution obtained using the present FEM. Table 4.2 shows a comparison of the reduced free energy of two lone particles as calculated by the FEM using equation (2.17) and as calculated by Hoskin [21] and Hoskin and Levine (HL) [28]. Note that the values in Table 4.2 represent the double layer free energies for two spheres. The values for the particles at infinite separation were calculated by Hoskin while the remaining values in the HL column were calculated by HL. Hoskin calculated the free energy for an isolated particle using a similar expression as equation (2.17) but he solved the one-dimensional integral problem. In the case of two interacting particles, the double layer free energy values in the HL column of Table 4.2 were obtained by taking the sum of the change in free energy due to the interaction (which was calculated by integrating the force results) and the free energy of the isolated particles calculated by Hoskin. HL also calculated the free energy of the double layer surrounding the interacting spheres by integrating equation (2.17) over their finite difference grid, but reported that their results were unsatisfactory. The failure of this method in the HL treatment was attributed to the uneven node distribution around the particles (because of the dipolar transformation employed) and the error incurred in approximating the steep potential gradients near the particle surface. In this work, the double layer free energy due to a lone sphere was calculated using the two-dimensional finite element grid in order test the method before extending it to more complex cases involving interfacial particles. As can be seen in Table 4.2, there is very poor agreement between the FEM results and the previously published results for the double layer free energy of the lone and interacting spheres. The obvious conclusion Chapter 4. Results and Discussion 74 Table 4.2: Double Layer Free Energy (G'/4ir) of Two Spheres with Equal Constant Surface Potentials in a Single Electrolyte FEM % Diff. HL 47.8 2 oo 24.37 1 16.49 41.7 21.39 0.5 15.1 42.4 22.35 1.0 15.7 22.4 oo 73.2 4 89.60 17.5 68.4 80.40 0.5 83.17 18.0 1.0 70.5 5 2 oo 255.8 336.7 31.6 31.2 324. 0.5 247. 328.4 30.8 1.0 251. 1340.2 7.63 oo 4 1245. 7.94 1315.8 0.5 1219. 7.97 1328.6 1.0 1230. oo 2070. 2681.4 29.5 15 2 2645.4 29.4 .5 2045 oo 4 10,370 10,951 5.61 6.43 10,291 10,952 0.5 FEM : Present Finite Element Method. HL : Hoskin [21], Hoskin and Levine [28] (see text for details). Kd V'o /c[i2 - 2a] Chapter 4. Results and Discussion 75 is that the present FEM suffers from the same problems as HL's finite difference method did when equation (2.17) is used to calculate the double layer free energy. The present FEM uses a much higher node density than that used by HL, and the nodes are certainly distributed more evenly throughout the double layer in the present FEM scheme. The node distribution around the particle surface is not uniform--in the present design (see Fig. 3.1), however, and this is most likely the cause for large discrepancy in Table 4.2. The inability of the present method to adequately determine the double layer energy may be overcome with a better finite element grid design; this possibility is discussed in Chapter 6. Recall that the change in the double layer free energy (AG) due to the interaction is the quantity that is of real interest. AG may be calculated using equation (2.18), i.e. taking the difference between the double layer free energy at a given separation (G<f) and the energy for two isolated particles (GQO). Since Gd and G ^ are similar (~ 10% difference for small separations), any errors occurring in the separate determination of Gj and Goo have their effect increased substantially when AG is calculated. Thus, the large errors in the Gd and G ^ values evident in Table 4.2 cause enormous errors (>100%) in AG when compared to the most recent literature values (ML [34]). Note that ML calculated AG by integrating the interaction force with respect to sphere separation and subtracting the free energy of the isolated spheres; they did not publish the separate Gd and Goo values. The present method of calculating AG for the double layer interaction (using equations (2.17) and (2.18)) is inadequate and at best provides an order of magnitude approximation to the true value. Integrating the interaction force with respect to particle separation has proved [28, 34] to be the most reliable numerical method for obtaining the interaction free energy. The present method was attempted because it is the only method available to approximate the interaction free energy for Model III (recall from Section Chapter 4. Results and Discussion 76 2.4.3 that there is no net parallel force acting on the particle iri the unit cell representation). The AC? results obtained by the free energy approach for Model III are discussed in Section 4.3. The double layer interaction energy for Model II was not evaluated using this method since the interaction force is available and hence the interaction energy can be determined by integrating the force. Due to time constraints, the force integration method for calculating the interaction energy was not attempted in this thesis for Models I or II, but its application is discussed in Chapter 6. 4.1.2 Constant Charge Interactions Constant charge interactions, by nature of the boundary condition on the particle surface, are more difficult to model than constant potential interactions. Recall from equation (2.14) that a surface charge is given by Gauss' law as the difference between the electric displacements, ^:(§f), on either side of the surface. Unless the normal potential derivative is zero on the particle side of the particle/electrolyte interface, then the potential distribution in the particle must be solved simultaneously with the potential distribution in the electrolyte. When the Derjaguin method is used to approximate a constant charge interaction between particles, the differential pairs of parallel, concentric rings (for spherical particles) on the two particles are assumed to act as portions of parallel, semi-infinite plates. The potential gradients in a semi-infinite plate (with no charge source or sink, i.e., when Laplace's equation governs the potential distribution in the solid) with a uniform constant surface charge (or potential) will always be zero. Thus, when employing the Derjaguin method to approximate constant charge interactions betweenfiniteparticles, an implicit assumption is made that the component of the potential gradient normal to the particle surface is zero. Whether or not this implicit assumption is a reasonable approximation is examined below. A series of constant charge interactions between identical spheres in an electrolyte Chapter 4. Results and Discussion 77 Figure 4.1: Fraction of 'Complete' to 'Incomplete' Force Solutions for Constant Charge Interactions between a Pair of Spheres in a Single Electrolyte were evaluated where the interior of the particle was ignored (as in the Derjaguin method) and where the potential distribution was simultaneously solved along with the distribution in the electrolyte. For the latter cases, the dielectric constant of the particle was varied to see what effect this would have on the interaction. Various values of KCL were also considered. For the sake of brevity, the solutions obtained by considering the interior of the particle will be referred to as 'complete' solutions while the solutions obtained by ignoring the interior of the particle will be referred to as 'incomplete' solutions. Fig. 4.1 plots the force acting between spheres (V>oo = 2) with different dielectric constants as a fraction of the result obtained when the force is calculated by ignoring the interior of the particle. In all three cases, the 'complete' solution produces results which are less than or equal to the 'incomplete' solution. The effect of the dielectric constant of the particle is apparent; a larger dielectric constant causes a greater deviation between the 'complete' and 'incomplete' solutions especially at small particle separations. Smaller Chapter 4. Results and Discussion 78 values of na also cause a greater difference between the 'complete' and 'incomplete' solutions at small particle separations. These discrepancies between the 'complete' and 'incomplete' solutions are consistent with Gauss' law as is discussed below. The magnitude of the force between the spheres increases with the magnitude of the electrostatic potential in the double layer overlap region; recall from Section 3.2.1 that the osmotic term in equation (2.34) accounts for ~95% of the total interaction force and the osmotic term in turn is directly dependent on the electrostatic potential. Thus, as shown in Fig. 4.1 the 'complete' solution results in a reduction of the electrostatic potential in the double layer overlap region when compared to the 'incomplete' solution and this potential reduction increases with the dielectric constant of the sphere; This behaviour is consistent if we examine Gauss' law; equation (2.14). The surface charge is given as the difference between the electric displacements on either side of the sphere surface. When the spheres are brought together the potential on the sphere surface must increase (for a positive charge). For the 'incomplete' solution, the surface charge is held constant by maintaining the potential gradient normal to the particle surface on the electrolyte side of the particle surface only. For the 'complete' solution, the surface charge is held constant by maintaining the difference of the normal potential gradients on either side of the particle surface. The potential gradients on either side of the particle surface are opposite in sign (i.e. the potential is decreasing away from the surface on both sides) so the surface charge is actually expressed as the sum of the magnitudes of the potential gradients on either side of the particle surface. Since there is a non-zero gradient within the particles (at small separations) the potential gradient on the electrolyte side of the particle surface will decrease from the value it had when the particles were at infinite separation. The net result is that the potential on the surface of the particles does not increase as much for the 'complete' solution as for the 'incomplete' solution. Another interesting feature of Fig. 4.1 is that the 'complete' and 'incomplete' solutions Chapter 4. Results and Discussion 79 are essentially equal at a separation of K[R — 2a] = 5 except for the na = l,e p = 80 case. This is consistent with the superposition approximation which assumes that at large separations, the interacting particles do not influence each other in any way and the potential distribution in the double layer overlap is given simply as the linear superposition of the potential distributions of two isolated spheres. If the particles do not influence each other (significantly) in a constant charge interaction then the surface potential over the whole particle surfaces should be constant, and thus the 'complete' and 'incomplete' solutions should be equivalent. Bell et. al. [34] found that the superposition approximation is valid when K[R — 2a] > 4 for constant potential interactions. The results in Fig. 4.1 indicate that this same general rule should apply to constant charge interactions providing na > 1 and e < 20. p It is interesting to note that McQuarrie et. al. [42] assumed that their model particle had a dielectric constant of 80 (same as water) which meant that electrostatic image effects could be ignored between the electrolyte and particle phases. Nullifying image effects however, does not preclude the applicability of Gauss' law which McQuarrie has done. Ironically, as the results presented above show, assuming a dielectric constant of 80 for the particle results in a greater error than assuming the typical dielectric constant values associated with most solids (i.e., 2-25). Since McQuarrie never solves for the potential distribution within the particle, his choice of particle dielectric constant is a moot point, but his assuming that e = e is the appropriate choice to make when p w employing the Derjaguin approximation is incorrect. The results of the present method were also compared with the Derjaguin approximation for constant charge interactions between equal spheres in a single electrolyte. Recall from Chapter 2 that the Derjaguin approximation for constant charge interactions between spheres is given by equation (2.24). The parallel plate interaction energy term, 80 Chapter 4. Results and Discussion V (5'), in equation (2.24) depends on the particle surface potential at infinite separar+ tion, V'oo, and the surface potential at separation S, xp (see equation (5.15) in [38]). The quantity xp is dependent on S, V'oo an£ f the particle surface charge, cr*, through a compli- cated, implicit expression involving elliptic integrals (see equation (6.4) in [38]). Thus, the Derjaguin approximation requires that xp be determined for a given set of 5,'0 ,cr" oo by finding the root of equation (6.4) in [38]. This was done for the results discussed below using the Newton-Raphson technique with a tolerance of 10~ . The repulsive force 6 acting between the particles is obtained by calculating V (S) using equation (5.15) in + [38] and then substituting this value into equation (2.24). When employing the Derjaguin approximation for constant charge interactions between spheres, one must choose which value of V'oo to use. In the Derjaguin approximation xpoo is related implicitly to <x" by equation (2.25), which governs the charge/potential relationship at the surface of a semi-infinite flat plate. The surface charge/potential relationship of a sphere is dependent on the size of the sphere and is given by equation (2.13). When constant potential : constant charge interactions are compared it is common practice to assume that the spheres in each interaction will have identical surface potentials at infinite separation, i.e. xp = 0 xpoo,c (where xp is the surface potential of the 0 spheres in a constant surface potential interaction and V'oo.c is. the surface potential of the spheres at infinite separation in a constant charge interaction). The question when employing the Derjaguin approximation is whether to use V'oo = Vo ' in the approximation (which corresponds to a lower cr" on the differential parallel plate segments than on the sphere) or to use V'oo = V'oo,*, where V'oo,. is calculated from equation (2.25), using the a" value of the spherical particle? McQuarrie et. al. [42] revealed that the Derjaguin approximation showed better agreement with their numerical results when V'oo = Vo ' was chosen. McQuarrie presented his results in graphical form showing the Derjaguin approximations obtained with -0°° = 81 Chapter 4. Results and Discussion ip 0 and ipoo = V'oo,* as separate curves on a graph of /** vs K[R — 2a] and his numerical results as points on the graph for the case of na = 5, cr* = 2.78. All of the numerical results showed good agreement with the V'oo = V'o = 2 curve (i.e. < 1% when read off his graph). There is one peculiarity in McQuarrie's results; his choice of a'. For a sphere with na = 5 and ip = 2, a* should equal 2.72 not 2.78. The value of 2.72 is obtained 0 from equation (2.13) for a sphere with na = 5 and V ' o o / : = 2; this value also agrees with the numerical results published by Loeb et. al. [17]. From equation (2.13), a" = 2.78 corresponds to a reduced surface potential of 2.033 for an isolated sphere. The Derjaguin approximation with V'oo = 2.033 yields results within .5% of the results obtained when tpoo = 2.0 so this error by McQuarrie is not very significant. This same numerical experiment was performed for this thesis and the results for the 'incomplete' solution (analogous to McQuarrie's numerical approach) are shown in Fig. 4.2 for two particle sizes. The solid lines are the Derjaguin solutions obtained with V'oo = V'o = 2 and the dashed lines are the Derjaguin solutions obtained with V'oo These results correspond to spheres with V'oo.c = V'oo,* • = 2 so for the na = 5 case; cr* = 2.723 and V'oo,* = 2.231, and for the na = 20 case; cr* = 2.443 and V'oo,* = 2.059. The circles in Fig. 4.2 represent the numerical results obtained in this work. As can been seen in Fig. 4.2, the numerical values give marginally closer agreement with the Derjaguin approximation when the choice of V'oo — V'o is made. The na = 5 case, when compared with Fig. 4 in [42], does not show as good agreement as McQuarrie's numerical results did with the Derjaguin approximation with V'oo = V'o- I* this work, the 1 numerical results are slightly greater than the Derjaguin approximation with V'oo = V'o (solid line) in the range 1 < n\R — 2a] < 2. The na = 20 case in Fig. 4.2 shows this same behaviour as well. The greatest difference between the FEM results and Derjaguin (with V'oo = V'o) approximation is ~4%. A series of tests were performed where the element densities were varied in the region close to the particle surface but no improvement on the Chapter 4. Results and Discussion 82 Figure 4.2: Comparison of 'Incomplete' Numerical Solutions and Derjaguin Approximations for the Interaction Force between a Pair of Spheres in a Single Electrolyte 83 Chapter 4. Results and Discussion results in Fig. 4.2 was obtained. Nonetheless, the results found in this thesis, although they appear to be of inferior quality, support McQuarrie's conclusions. It is interesting to note that for the na = 20 case the Derjaguin approximations using either choice of V'oo (2.0 or 2.059) yield force values that are very similar. Thus, for cases where na > 20, the question of which V'oo to use is meaningless since the flat plate and sphere results become identical. 4.1.3 Summary The results of this section are summarized below; all but the last point are based on the force calculations: • the constant surface potential results show good agreement with the previously published numerical data for constant potential interactions between identical spheres, which indicates the present FEM is acceptable for simulating constant potential interactions, • the constant surface charge results show adequate agreement with the Derjaguin approximation for cases where the approximation applies, further indicating that the present FEM is acceptable for simulating constant charge interactions, • the constant charge results support the findings of McQuarrie et, al. [42] as to which choice of V'oo is correct when applying the Derjaguin approximation to constant charge interactions between spheres. It was also revealed that this choice only makes a substantial difference to the predicted force results when na < 20, which is approaching the limit (na ~ 5) where the Derjaguin approximation is valid, • the implicit assumption — 0 in the interior of the particle) made in the Der- jaguin approximation for constant charge interactions introduces substantial errors 84 Chapter 4. Results and Discussion which increase as na and the particle separation decrease and as the dielectric constant of the particles increases, • the 'complete' and 'incomplete' solutions for the constant charge interactions become identical at large particle separations (except for na = 1 and e = 80), p indicating that the linear superposition approximation should apply to constant charge interactions between identical spheres for the same particle separations as it does for constant potential interactions, i.e. n[R — 2a] > 4, and • the double layer free energies obtained using equation (2.17) are very inaccurate probably due to the present grid design. 4.2 T w o Identical Cylinders E m b e d d e d in an O i l / W a t e r Interface This section presents the results obtained for Model II which was described in Section 2.4.2. A number of cases were studied: constant particle surface potential interactions where the oil/water (O/W) interface remained uncharged, interactions where both the particle surface and the O/W interface potentials were held constant, and constant particle surface charge interactions where the O/W interface remained uncharged. The parallel force acting on the cylinders was evaluted for all these cases but, for reasons explained later, the normal (to the O/W interface) force acting on the cylinders was only evaluated for those cases where the O/W interface was uncharged. It is convenient to discuss these various cases in separate sections. 4.2.1 Constant Surface Potential Cylinders E m b e d d e d in an Uncharged O / W Interface This section presents the results for the cases of identical, interacting cylinders with constant surface potentials which are embedded in an uncharged O/W interface. The Chapter 4. Results and Discussion 85 material presented in this section serves three purposes: • to examine the ability of the Derjaguin approximation to estimate the interaction force between particles embedded at an 0/W interface where only the portion of the particle immersed in the aqueous phase is considered, • to compare the interaction forces acting parallel and normal to the interacting cyhnders, and • to examine the course of the interaction over a wide range of particle separations to see whether the behaviour agrees with previous predictions (see Section 2.3). Accuracy of the Derjaguin Method The Derjaguin approximation for calculating the repulsive force between identical, parallel cyhnders at an 0/W interface was discussed in Section 2.3. The method presented in Section 2.3 is only valid for low potentials, i.e. when the Debye-Hiickel (DH) approximation applies. High potential solutions involving Derjaguin's method were not attempted as they involve numerically differentiating, then integrating elliptic integrals as part of an involved root finding procedure; this is not impossible but it was deemed unnecessary since the low potential approximation should indicate whether or not Derjaguin's method provides an adequate estimate of the interaction when compared with the low potential numerical results. It should be noted that Levine et. cd. [1] only used the DH approximation when applying the Derjaguin approximation to pairs of equal spheres embedded in an O/W interface. A comparison of the low potential Derjaguin approximation (solid lines) and the numerical results (symbols) obtained with the present FEM for the parallel force acting between cyhnders embedded in an uncharged 0/W interface is shown in Fig. 4.3 for Chapter 4. Results and Discussion 86 Figure 4.3: Low Potential Derjaguin Approximation Compared with Numerical Predictions of the Force Acting between Identical, Parallel Cylinders Embedded in an Uncharged 0 / W Interface Chapter 4. Results and Discussion 87 four values of /ca, all with V'o = -5. Results are plotted for four different contact angles. The comparison between the FEM results and the Derjaguin predictions is good for all cases except for /ca = 5 where there is ~10% difference at the smallest separations. This same trend is present to a lesser degree as /ca increases, i.e. the FEM results are slightly larger than the Derjaguin predictions at small separations. The probable explanation for this discrepency is that at small separations the electrostatic potential in the overlap region will be very similar to V'o and thus the neglected terms in the DH approximation could be significant. Linearizing the right-hand side of the Poisson-Boltzmann equation (PBE), equation (2.10), and retaining only the first order term results in an error of ~ 4% when ip = .5. The most probable reason the /ca — 5 numerical results show such poor agreement with the Derjaguin approximation is that at this small value of /ca the assumptions inherent to this method may no longer be valid. Recall that the limit when applying the Derjaguin's approximation to interacting spheres is na >5 [28] and since the double layer does not decay as quickly for cylinders as it does for spheres it is plausible that the /ca limit for applying the method to cylinders would be > 5. This argument is qualitatively supported by the results of Hoskin and Levine [28] who showed that the Derjaguin approximation tmder-estimates the interaction force for spheres for small /ca and separations, K[R — 2a] < 1. Perusing Fig. 4.3, one appreciates that the majority of the interaction comes from the overlapping double layers in the region about the plane joining the axes of the cylindrical particles. This feature is most prominent in the na = 30 case where the results for contact angles (6) of 45° and 67.5 0 are almost identical. It is also evident that for 6 > 112.5° the effect of the interaction is negligible, i.e the double layers do not have an opportunity to overlap. This effect is seen even more clearly in Figs. 4.4 and 4.5 which show the equipotential contours for two interacting cylinders with constant surface potentials on an uncharged interface at 45° and 135° respectively; both figures are for cylinders with Chapter 4. Results and Discussion 88 na = 10, rpo — 2 and K[R — 2a] = 5. Thesefiguresmay also be compared with Fig. 3.5 which is a contour plot of the same interaction but with 6 = 90° . Recall that the interaction force is calculated by integrating the electric field and osmotic pressure over the median plane separating the two interacting cyhnders; this median plane corresponds to the left-hand border of the plots in Figs. 3.5, 4.4, and 4.5. As seen in the equipotential plots, the double layer overlap (when possible) is most prominent in the region about the plane joining the axes of the two cylindrical particles so it is obvious why the majority of the interaction is contributed from this region. The Derjaguin method employed by Levine et. al. [1], and used in this thesis, completely ignores any interaction acting through the oil phase and any effects due to the boundary condition governing the potential distribution on the O/W interface. An examination of the contributions to the overall interaction force calculated by the present FEM reveals that for small separations (n[R— 2a] < 5) the major contribution arises from the aqueous phase (> 99%). Fig. 4.3 clearly indicates that (for low surface potentials) the DH form of the Derjaguin method provides a good approximation for the interaction of a pair of parallel cyhnders embedded in an uncharged O/W interface providing the conditions associated with the Derjaguin approximation apply, i.e. the double layer thickness and particle separation are both significantly smaller than the particle size. Comparison of the Parallel and N o r m a l Forces The forces acting parallel and normal to a pair of interacting cyhnders embedded in an uncharged O/W interface were calculated for several different cases. The dependencies of the normal and parallel forces on model parameters are completely different. Some discussion on the parallel force was given in the previous section, so only a brief commentary on the parallel forces for the cases considered in this section will be provided. The effect of various model parameters on the normal force are examined below. Finally, the Chapter 4. Results and Discussion 89 EQUIPOTENTIAL CONTOURS o Figure 4.4: Potential Distribution for One of Two Parallel Cylinders on an Uncharged 0/W Interface, 0=45° Chapter 4. Results and Discussion 90 EQUIPOTENTIAL C O N T O U R S o iri >- XX Figure 4.5: Potential Distribution for One of Two Parallel Cylinders on an Uncharged O/W Interface, 0=135° 91 Chapter 4. Results and Discussion "S—B—B—B •e B- B- €1 12 15 CM - o CN 0 3 6 9 ICQ Figure 4.6: Normal Force acting on an Isolated Constant Potential Cylinder (with variable no) Embedded in an Uncharged 0/W Interface effect of the interaction on the relative magnitudes of the two forces will be compared. The parallel force results are presented in Appendix D in Figs. D.l to D.4. Three values of na (2,5,10) are considered, each at three separate contact angles (45°,90°,135°). The parallel forces for all three fta cases were obtained with particle surface potentials of 2 and the force for the na = 10 case was also evaluated with tp =4. These results display 0 the same trends evident in the low potential results which were discussed in the previous section. It was revealed in Chapter 3 that the normal force acting on a charged cylinder at an O/W interface is present even when the cylinder is isolated. Before discussing the interaction effects it is sensible to examine the factors that cause variations in the normal force acting on isolated cyhnders. The magnitude of the normal force, for an isolated particle, is dependent on two geometric parameters; na and the contact angle, 6. Fig. 4.6 shows the variation of Chapter 4. Results and Discussion 92 the normal force for cylinders as a function of na where the cylinders are embedded (6 = 90°) in an uncharged O/W interface and the cylinder surface in contact with the aqueous phase has V'o=2. For KO > 6 the normal force is virtually independent of na. Recall that the normal force is calculated by integrating equation (2.40) over the O/W interface. The potential at the O/W interface decreases rapidly as a function of distance from the particle (see contour plots in Figs. 3.5, 4.4, 4.5) and consequently > 95% of the total contribution to the integral comes from the region of the interface within a couple of double layers of the particle surface. Fig. 4.6 indicates that this integral will be virtually equivalent for all cases with na > 6. This occurs because the thickness of the double layeris so much smaller than the radius of curvature of the particle that the intersection of the double layer (extending from the particle surface) with the O/W interface resembles the intersection of two planes. In other words, for na > 6 the normal force acting on the cylindrical particle is equivalent to the normal force acting on a semi-infinite planar charged surface (with the same 6 as the cylinder) immersed in the same O/W interface. It is not surprising that this phenomenon occurs at approximately the same value of na that is the lower limit of the Derjaguin approximation. Fig. 4.7 shows the variation of the normal force as a function of 6 for an isolated cylinder with na = 10 and = 2. There is a minimum in the force at 6 — 90° and it is substantially greater for 6 < 90° than for 6 > 90°. This behaviour can be rationalized by examining the contributions of the separate terms in equation (2.40) and the contour plots in Figs. 3.5, 4.4, 4.5. Recall from Chapter 3 that the relative contributions to equation (2.40) for a case with na = 10 and 6 = 90° are; osmotic term ~ 50%, (§^) term ~ 50%, and (§j) term ~ 1%. The osmotic term in equation (2.40) accounts for 2 ~ 92% of the total normal force when 9 = 22.5° but only ~ 30% of the total when 6 = 157.5°; the absolute change in the osmotic term from 157.5° to 22.5° is ~9 times. At low contact angles, the double layer, which extends in a direction normal to the particle 2 93 Chapter 4. Results and Discussion CO o •e- o 30 60 90 120 Contact Angle 6 [deg] 150 180 Figure 4.7: Normal Force acting on an Isolated Constant Potential Cylinder (with variable 6) Embedded in an Uncharged 0/W Interface surface, is projected onto the O/W interface. As seen in Fig. 4.4, the net effect is that the potential on the 0/W interface does not decay as rapidly as a function of distance from the particle surface as it does when 8 = 90°. At high contact angles, the double layer is projected away from the O/W interface. As seen in Fig 4.5, the net effect is that the potential on the 0/W interface decays more rapidly as a function of distance from the particle surface for high contact angles than it does when 6 = 90°. This reduces the contribution of the osmotic term but increases the contribution of the (^) 2 term; for the same case discussed above, when 6 = 157.5° , the (§£) term accounts for 67% of the 2 total normal force. Fig. 4.8 contains plots of the normal force acting on one of a pair of parallel, identical cyhnders (with ip = 2) for various separations, contact angles and values of na. Fig. D.5 0 in Appendix D plots the same quantity for the case of na = 10 and xp — 4. In almost 0 every case the interaction causes the normal force to decrease. While Fig. 4.8 plots Chapter 4. Results and Discussion KQ 0.0 1.0 2.0 = 2. 3.0 c[R-2a] 94 o 4.0 5.0 ICQ 0.0 1.0 2.0 = 5. 3.0 e[R-2a] o 4.0 ' 5-0 0.0 KO 1.0 2.0 = tO. 3.0 «:[R-2a] 4.0 3.0 Figure 4.8: Normal Force Acting on Constant Potential Cylinders Embedded in an Uncharged 0/W Interface the total normal force ('front side' plus 'back side'), only the 'front side' contribution is affected in these constant potential interactions. The potential distribution on the 'back side' of the cylinder remains unaffected by the interaction; this may be observed in all of the contour plots for constant potential interactions presented in this thesis. The only cases where the normal force is not appreciably affected by the interaction are when na = 10 and 0 ^ 90°. Again, for this value of na, the double layer is very small compared to the radius of the cylinder so the interaction has little effect at the O/W interface (see Figs. 4.4 and 4.5). Examination of the separate contributions to equation (2.40) revealed that the osmotic term on the 'front side' of the cylinders increased as the separation between the particles decreased, but the reason that the overall normal force is reduced is that the 'front side' (^) 2 term (the electrostatic stress) decreases by up to 95% for small separations. This may be appreciated by examining the double layer overlap region in Fig. 3.5, i.e. the potential varies very slowly along the interface which 95 Chapter 4. Results and Discussion is in the overlap region. Another interesting feature in Fig. 4.8 may be seen in the na — 2 plot; the force is greater at 135° than at 45°, which is opposite to the trend shown in Fig. 4.7. The separate contributions to the total integral in equation (2.40) were examined and it was found that when na < 5 the (f^) term dominates over the osmotic term (the opposite 2 is true at larger KCX'S). For the case of na = 2 and 8 = 135°, the (§^) term accounts for 2 ~ 80% of the total 'back side' force. Thus, the electrostatic stress term dominates over the osmotic term when the double layer is comparable in size to the radius of curvature of the particle. All the normal force values calculated in this work are positive. Thus, according to the derivation given in Appendix A, there is a force normal to the O/W interface acting on the cylinder that is directed into the water phase. In other words, when a charged particle is embedded in an O/W interface there will be an electrostatic/double layer interaction that tends to force the particle into the electrolyte phase. For all the constant potential cases studied in this thesis, the normal force per unit length (K~ ) of the cylinder (/**) varied between 2-3 (except for the case at 22.5° in X Fig. 4.7 which had / " ~ 4). Recall from the beginning of this chapter that one unit of the reduced force, /*, is equivalent to ~ 2 X 10 iV. Assuming that the cylinder is _13 composed of amorphous silica with a density of 2200 kgrn~ [55], and is embedded in 3 an O/W interface where: the electrolyte concentration is 10 -4 molar, the system is at room temperature (/c % 3 x 10 m), the cylinder has a radius of .Zpm (so KCL = 10), _1 _8 and V'o = 2(~ SOmV^), then the gravitational force acting per unit length on the cylinder is ' 2 10 N. From equation (4.2) (/ x 16 = 4h) the normal force due to electrostatic effects will be ~ 2 x 10 iV or ~ 10 times the gravitational force! Since the normal —14 4 force remains constant, except for very small particles (assuming n remains constant), the ratio of the gravitational/electrostatic forces will increase with a , where a is the 2 96 Chapter 4. Results and Discussion cylinder radius. If the above example is reworked using a cylinder radius of 30/im and it is assumed that the normal force remains constant, then the gravitational force is approximately equal to the electrostatic normal force. For larger particles or a more concentrated electrolyte the gravitational force quickly becomes the dominant normal force acting on the cylindrical particle. Recall (Section 2.4) that this analysis is meant to apply to particles so small that gravity effects can be ignored and the interface is assumed to be flat. The above results indicate that the electrostatic normal force may be large enough to affect the curvature of the interface in the vicinity of the particle. Whether or not this occurs depends on a number of factors such as the relative densities of the three phases and the 0/W interfacial tension. This possiblity will be addressed further in Chapter 6. The parallel force acting between two interacting cylinders depends strongly on the particle separation, the particle size and the particle contact angle; see Figs. D.1-D.4 in Appendix D. Assuming the same conditions as above, the reduced force per unit length acting between two cylinders with na = 10, 6 = 45°, ipo — 2 and at a separation of K[R — 2a] = .25 is ~ 9 (= 7 x 10 iV~) or about 4 times the electrostatic normal force -12 acting on one of those cylinders. The quantity that is of more interest is the free energy change associated with the interaction, which is obtained by integrating the change in the force as a function of particle separation. Unfortunately this was not performed for this thesis, but perusing the results for the normal and parallel forces as a function of particle separation, it is obvious that the parallel interaction (i.e. the double layer repulsion) will have a much greater free energy change than the normal interaction. The normal force is only slightly affected by the interaction and it quickly reaches a constant value at small particle separations, thus its integral will be very small. The parallel force, however, declines from a maximum to zero as the particles are separated so the free energy associated with this interaction will Chapter 4. Results and Discussion 97 be substantially greater. Effect of Particle Separation on the Parallel Force It is believed [10, 45, 46] that charged particles trapped at an interface between an electrolyte and a non-polar phase will interact through both phases. There is a double layer interaction acting through the electrolyte that exponentially decays with particle separation. The particles are thought to interact through the non-polar phase according to a dipolar force law which decays algebraically' as a function of the particle separation. The parallel force was calculated over a wide range of particle separations to discover the nature of the long range interaction between charged cyhnders at an O/W interface. The six plots in Fig. 4.9 show the interaction between two cyhnders with na = 10 and tpo = 2 over a wide range of particle separations. Semi-log plots are given on the left-hand side of the figure and log-log plots are on the right-hand side. The force contributions from the oil and water phases to the total force are plotted separately for three different contact angles. The data in Fig. 4.9 are also presented in numerical form in table D.l in Appendix D. The plots in Fig. 4.9 clearly indicate a short range double layer interaction acting through the electrolyte and a weaker, but much longer ranging interaction acting through the non-polar phase. The contribution from the water phase dominates the short range force (except for the 6 = 135° case where double layer overlap is negligible) while the contribution from the oil phase dominates the long range force. The double layer interaction decays exponentially with separation; the semi-log plots of the interaction yield straight lines at small separations (except for 6 = 135°). The long range interaction is not purely algebraic as is predicted from previous theoretical treatments (see Section 2.3), i.e. instead of a straight line, the log-log plots of the Chapter 4. Results and Discussion 98 Contact Angle: 45. Figure 4.9: Parallel Interaction Force acting between Two Parallel Cylinders Embedded in an Uncharged 0/W Interface Chapter 4. Results and Discussion 99 long range interaction show a mild curvature. Grid sensitivity tests ensured that this effect was not a function of the element density between the interacting particles. Another interesting feature of the log-log plots is that the long range interaction acting through the electrolyte shows a transition from an exponential decay to a pseudo-algebraic decay at na % 50. These two features are not predicted from the current [10, 45, 46] mathematical treatments which all employ simplifying assumptions (such as the Debye-Hiickel approximation) to yield tractable solutions. It is possible that terms discarded in these simplifying assumptions are responsible for this unpredicted behaviour or the assumptions-are not valid in the case of infinite interacting cyhnders. Further speculation on this point is beyond the scope of this thesis, but thesefindingsmay warrant further investigation. 4.2.2 Constant Surface Potential Cylinders E m b e d d e d in a Constant Potential O / W Interface It is more likely than not that an O/W interface will be charged. The charging mechanism will undoubtedly involve an adsorption of some species from either phase and thus, assuming kinetic equilibrium, a constant potential boundary condition will govern the potential distribution about the interface. A series of simulations were performed where both the O/W interface and the cylinder surface in contact with the electrolyte had constant potential boundary conditions. Before discussing the results, the manner in which these boundary conditions were imposed will be outlined. The cylinder surface constant potential boundary condition was imposed in the same manner for these simulations as for all previous simulations, i.e. all nodes on the cylinder surface in contact with the electrolyte (including the node at the point where the three phases meet) were assigned a constant potential, tp . The O/W interface boundary 0 condition was imposed in the same manner, i.e. all the nodes on the O/W interface Chapter 4. Results and Discussion 100 (excluding the node at the point where the three phases meet) were assigned a constant potential, V'om- The surface potentials, V'o and V'oto, are not necessarily equal so there may be a step jump in the nodal potential values between the point on the particle where the three phases meet and the node on the interface closest to the particle. This situation is not realistic as the presence of the particles on the interface would perturb ^ou; in the vicinity of the particles, but it does provide a first approximation for modeling the parallel force between the cylinders. Recall that the normal force calculation is very dependent on the potential distribution about the O/W interface in the vicinity of the particles and thus the force calculation strongly depends on the interfacial boundary condition imposed in this region. For this reason, the normal force calculations were not performed for the interactions described in this section. Recall from Chapter 2 that to calculate the parallel force between the two cylinders, one must integrate the left-hand side of equation (2.22) over: the median plane in both fluid phases, the O/W interface and, the surfaces at infinite distance from the interacting particles. When the O/W interface is uncharged the contributions to the parallel force from the infinite surfaces vanish. When the O/W interface has a constant potential boundary condition, however, the contribution at the infinite surface, in the vicinity of the interface, must be subtracted from the force calculated at the midplane. Interactions were modeled with the following parameters; na = 2,5,10, 6 = 45°, 90°, 135°, ipcnu = -2,-1,0,1,2, /c[R-2a] = .25-5., and V'o = 2. Additional simulations were performed for the na = 5 case where Vw = -4,-3,3,4. The bulk of the results are presented in Figs. D.6-D.13 in Appendix D. The general trends in the results are discussed below. The major differences between the V'o ~ Tpow a n d the V'ouncharged interface inter- actions studied in the previous section are: at large separations there is no long range interaction acting through the oil phase, and at small separations, if the interface influences the overlapping double layers, the interaction force is dependent on ipow a n d may 101 Chapter 4. Results and Discussion *[R-2a] «[R-2a] /c{R-2a] Figure 4.10: Parallel Interaction Force acting between Two Parallel Cylinders Embedded in an O/W Interface with a Constant Potential be larger or smaller than the corresponding interaction on an uncharged O/W interface. The former difference arises because at large separations the potential distribution in the oil phase is unaffected by the presence of the particles. As Laplace's equation governs the potential distribution in the oil phase, which is assumed to be semi-infinite, then the potential in the oil phase far from the particles will be constant throughout at ipow- Any long range dipolar interaction acting through the oil phase is swamped out by the presence of the charged O/W interface. The latter difference requires more elaboration and is discussed below. The O/W interface will not always influence the short range magnitude of the parallel force acting between the pair of cyhnders. Fig. 4.10 shows the interaction forces for three different pairs of cyhnders (na — 2,5,15) all with a contact angle of 45°. As na increases, the force values for the various ip^ cases become equal (in fact they are actually all Chapter 4. Results and Discussion 102 approaching results found for the uncharged interface interactions). This is the familiar effect of na already referred to above in the discussion on normal forces; when the double layer thickness is small compared to the radius of curvature of the particle, then at low contact angles, the interface will not influence the overlapping double layers. This explains why the force values for the na = 15 case are all equal, independent of VwWhen the interface does influence the interaction, the magnitude of the interaction will depend on Vw- Assuming V'o is positive, then when Vw is equal to V'o the interaction is minimized and it gradually increases as is decreased. This behaviour may be explained if the contributions to the total parallel force are examined, i.e. equation (2.38) (equation (2.39) is not considered because the oil phase contribution to the parallel force is negligible). When Vw — V'o, the osmotic term in equation (2.38) accounts for 95% of the total parallel force acting between the two cylinders. As Vw is decreased from V'o, the ( § j ) term increases: for the case of tp = 2, the (f^) term increases by a factor of ~40 2 2 0 when •Vow as decreased from 2 to -2. This increase in the electrostatic stress component w of the force is not surprising considering that at small separations the potential along the midplane must vary from ip ow at the interface to a value close to ip in the region where 0 the particle double layers overlap. It is also useful to compare these force results with those in Figs. D.1-D.2 for the uncharged interface case. For na = 2, n[R — 2a] = .25, and 9 = 45°, the force from Fig. D.l is 3.85 while the forces from Fig. 4.10 range from 2.88 to 8.11 depending on V'ow This same trend was observed in all cases where the interface influenced the overlapping particle double layers; the effect of the constant potential interface is to either lower or raise the interaction force, depending on the relative magnitudes and signs of V'ow a n d V'o, when compared to the corresponding interactions occurring on an uncharged O/W interface. When ip ow is greater than V'o the interaction is minimized further and when V'cu, is decreased below a negative value with the same magnitude as V'O) the force is 103 Chapter 4. Results and Discussion further increased; compare Figs. 4.10 and D.8 with D.12 and D.13. It is worthwhile to compare the V'o —V'o™ and the V'o—uncharged interface interactions in Fig. D.2 to D.13 as a guide to speculate on the free energy change associated with the V'o V w interactions. For the na = 2, 6 = 90° case, the force for the constant O / W — interface (tpow — — 2) interaction is ~9 times as great as for the same interaction on an uncharged interface: this is the largest proportional difference found in this data for the same interaction occurring on the two types of interfaces. Recall that the free energy is obtained by integrating the area under the force vs separation curves. The force curves for the V'o — V'ow interactions are in general much steeper and than the corresponding curves for the uncharged interface interactions; this is especially true for those cases where V w = — 1,-2, which result in the greatest interaction forces. Considering this, and the fact that there is no long range interaction when the O / W interface is charged, then it is reasonable to assume that for the cases simulated in this thesis, the free energy change associated with bringing two cyhnders together on a charged interface (with V>ow = -V'o) will at the most be 10 times greater than the same interaction occurring on an uncharged interface. Levine et. al. [1] have calculated the free energy change associated with forming a close-packed structure of spheres on an uncharged interface and made an upper estimate of 2x 10 kT per sphere. Multiplying this result by a factor of 3 10 (i.e. assuming that V'o™ = ~V'o) still leaves the free energy change due to electrostatic interactions two orders of magnitude below 10 kT. 6 4.2.3 Constant Surface Charge Cylinders E m b e d d e d in an Uncharged O / W Interface Interactions between constant charge cyhnders embedded in an uncharged O / W interface were evaluated for comparison with the constant potential interactions. Interactions involving constant charge cyhnders embedded on a constant potential O / W interface Chapter 4. Results and Discussion 104 were not performed because the boundary conditions imposed in the vicinity of the point where the O/W interface meets the particle surface lead to unstable solutions, i.e. a constant charge node on the particle surface and an adjacent constant potential node on the O/W interface are not compatible. Both the parallel and normal forces acting on the constant charge cyhnders embedded in an uncharged O/W interface were calculated and the results for both of these forces are presented and discussed separately below. Before discussing the force data it is worthwhile to peruse the potential distributions associated with each type of interaction to appreciate the differences between them. Fig. 4.11 shows the equipotential contour plots for one of pair of interacting cyhnders (na = 10, K\R — 2a] = 1, 6 = 90° and V'oo = 2) where the top plot corresponds to a constant potential interaction and the bottom plot corresponds to a constant charge interaction. The potential distributions in the water phase look virtually identical for both types of interactions at this plot resolution. Close inspection of the 'front side' contour lines in the water phase of both plots reveals some subtle differences; i.e the xp = 1.5 and xp = 1 contours in the constant charge plot intersect the median plane at slightly higher values of KY than the same lines on the constant potential plot. The contour hnes look similar in the aqueous phase for the two different interactions because the potential varies exponentially. The differences in the potential distributions between the two interactions are more noticeable in the solid and oil phases. Every one of the contour hnes in the sohd and oil phases in the constant charge plot is shifted away from the region between the interacting cyhnders when compared with the same contours in the constant potential plot. This indicates that overall, the potential in the region between the interacting cyhnders is, as expected, greater in the case of constant charge interactions than constant potential interactions. 105 Chapter 4. Results and Discussion ICHflPOlLNIIAt CONIOUWS XX Figure 4.11: Potential Distributions for Constant Potential (top plot) and Constant Charge (bottom plot) Interactions between Pairs of Identical Cyhnders Embedded in an Uncharged O/W Interface Chapter 4. Results and Discussion 106 Parallel Interaction Force The parallel force acting between pairs of identical, constant charge cylinders with V'oo = 2, /ca = 2,5,10 (and therefore cr* = 2.781, 2.529, 2.441, respectively), 6 = 45°, 90°, 135°, and K[R — 2a] = .25-5 are shown in the plots in Fig. 4.12. The plots in Fig. 4.12 may be compared with the corresponding constant potential interactions shown in Figs. D.1-D.3 in Appendix D. The same trends in the magnitude of the parallel force as function of /ca and 6 seen in the constant potential results are present in the constant charge results. The major difference between the two types of interactions is that the magnitude of the constant charge interaction forces at small separations (n[R — 2a] < 2) are much larger than the corresponding constant potential results. The reasons for this have already been discussed, i.e. the surface potential of interacting surfaces must increase as two surfaces with positive constant charges approach each other. At the smallest separations evaluated in this work (K[R — 2a] = -25), the constant charge interaction force is never more than four times greater than the corresponding constant potential force. The slope of the constant charge force curves are all steeper than the corresponding constant potential force curves at small separations. At large separations (K[R — 2a] > 4), the parallel forces for the constant charge and constant potential interactions become similar; this occurs because at large separations the interacting particles do not influence each other's double layers, (i.e. this is the common assumption used in the linear superposition approximation). Noting these facts, one may safely speculate that the free energy change for constant charge interactions will not be more than 4 times greater than the corresponding energy change associated with the constant potential interaction, for the cases studied in this thesis. 107 Chapter 4. Results and Discussion Figure 4.12: Parallel Interaction Force acting between Two Constant Charge Parallel C5linders Embedded in an Uncharged O/W Interface r N o r m a l Force Just as with the constant potential interactions, it is sensible to examine the factors that cause variations in the normal force acting on isolated cyhnders with constant surface charges before discussing the effect of the interaction. Recall that for an isolated particle the normal force is dependent on two geometric parameters; na and the contact angle, 6. Fig. 4.13 shows the variation of the normal force for constant charge cyhnders as a function of na when the cyhnders are embedded in an uncharged O/W interface (6 = 90°) and the constant surface charge, <x*, is calculated using equation (2.11) with V'oo = 2. Fig. 4.13 may be compared with the corresponding constant potential case, Fig. 4.6. Unlike the constant potential case, the normal force in Fig. 4.13 increases with na until na ~ 15. Once again the contributing terms to the total force, given by equation (2.40), were examined. For na > 3, the (§^) term remained relatively constant: its absolute 2 Chapter 4. Results and Discussion 108 value increased by less than 1% in the range of na = 3 to na = 20. The osmotic contribution to equation (2.40), however, increases by ~ 50% over the range /ca = 1 to /ca = 15; no such increase was noted in the constant potential case. An increase in the osmotic term implies an increase in the potential in the region where the cylinder and interface meet. Recall that the surface charge is calculated using equation (2.11) which relates a' of an isolated cylinder, completely immersed in a single electrolyte, to /ca and V'oo- When this same surface charge is imposed on a cylinder embedded in an O/W interface, it is not surprising to find that the surface potential varies over the surface of the cylinder. Examination of the surface potentials revealed that the potential was at a maximum at the points where the interface meets the cylinder. This suggests that the oil phase has an insulating effect on the electrostatic potential distribution which in turn causes the potential distribution along the O/W interface to be more sensitive to /ca than is found when the cylinder has a constant surface potential. Fig. 4.14 shows the variation of the normal force as a function of 6 for an isolated cylinder with /ca = 2 and a' = 2.441. This figure is also markedly different from the corresponding Fig. 4.7 for the constant potential case. Instead of a minimum at 90°, the normal force increases substantially as the contact angle is decreased from to 22.5°. The magnitude of the force at 22.5° is ~5 157.5° times greater when the constant charge case is compared to the constant potential case. Once again, examination of the contributions to equation (2.40) helps to explain this behaviour. When 6 = 157.5° percent contributions of the individual terms to equation (2.40) are: osmotic term ~ (ft ) term ~ 2 65%, and ( f ^ ) term ~ 2 5%, when 6 = 22.5° the 30%, the percent contributions are: osmotic term ~ 86%, (§£) term ~ 15%, and ( f £ ) term ~ 1%. The absolute value of 2 the (^) 2 term increases by a factor of ~ 2 4 from the 157.5° to the 22.5° case, but the absolute value of the osmotic term increases ~50 times over this same range! Thus, the constant charge boundary condition on the cylinder surface causes the potential in the Chapter 4. Results and Discussion 109 Figure 4.14: Normal Force acting on an Isolated Constant Charge Cylinder (with variable 6) Embedded in an Uncharged O/W Interface 110 Chapter 4. Results and Discussion «ro = 2- KO O 9 = 10. "5 ° E - Contact Angfe • = 45 o = 90 <• = 135 2.0 3.0 <e(R-2o] Contact Angle o = 45 o = 90 » = 135 2.0 3.0 «(R-2a] Contact Angle • = 45 o = 90 A = 135 2.0 3.0 <{R-2o] Figure 4.15: Normal Force acting on Constant Charge Cyhnders Embedded in an Uncharged O/W Interface region between the particle surface and the O/W interface to increase significantly and the increased osmotic pressure in this region produces a large force which acts to push the particle into the water phase i.e. the force is positive. Fig. 4.15 contains plots of the total normal force ('front' plus 'backside' forces) acting on one of a pair of parallel, identical cyhnders that have constant surface charges calculated using equation (2.11) for ^ o o = 2. Fig. 4.15 may be compared with Fig. 4.8 which shows the corresponding constant potential interactions. The most striking difference between Figs. 4.15 and 4.8 is that the constant charge interactions have a minimal effect on the magnitude of the normal force whereas there is a significant effect in the constant potential interactions. Additionally, the interactions for 8 = 45° actually show a shght increase in the normal force at small separations for the constant charge interactions while the constant potential interactions always showed either a decrease or no change at all. Chapter 4. Results and Discussion 111 As shown above, the normal force increases as the contact angle decreases. Examination of the separate 'front' and 'back' side contributions to the total normal force for the 8 = 45° cases revealed that both the 'front' and 'back' side forces increased when the cylinders were at small separations; e.g. for na = 2, and n[R — 2a] = .25, the 'front' side force increased by 15% compared to the value at infinite separation, while the 'back' side increased by 2%. When 8 — 90° the 'front' side force decreased slightly but the 'back' side force increased slightly. Examination of the separate contributions to equation (2.40) revealed that the osmotic term increased for both the 'front' and 'back' side forces ('front' side increase was always greater than the 'back' side increase), but the (§i£) term for the 'front' side force decreased substantially (~ 90%) due to the double 2 layer overlap. The constant charge boundary condition on the cylinder surface results in an increase of the cylinder surface potential (and therefore the whole particle) when the cylinders are at small separations. The results presented above for the 8 — 90° case indicate that the increased osmotic pressure associated with increased particle potential counter-balances the decrease in the electrostatic stress contribution (on the 'front' side) to the total force so that the total normal force acting on the cylinder remains relatively constant irrespective of particle separation. For the 8 — 45° cases, the decrease in the electrostatic term does not counterbalance the increase in the osmotic term as well as in the 8 = 90° cases so the overall normal force acting on the cylinders increases slightly at small separations. For 8 > 90°, the normal forces for the constant charge and constant potential cylinders are comparable, i.e. all between / " = 1 and /*" = 3. For 8 = 45°, the constant charge normal forces are approximate^ 2.5 times greater than the corresponding constant potential results. Since the normal forces for the constant charge situations show minor changes with particle separations, then the free energy change associated with the normal component of the force for interacting cylinders will be negligible compared to the free 112 Chapter 4. Results and Discussion energy change associated with the parallel force. 4.2.4 Summary The results of this section are summarized below in point form: • the numerical results for cyhnders with constant surface potentials (in the DebyeHiickel regime) which are embedded in an uncharged 0/W interface show good agreement with the corresponding Derjaguin method predictions. The results indicate that the Derjaguin method is acceptable for cases where Ka > 5, • the normal force acting on a constant surface potential cyhnder embedded in an uncharged 0/W interface exists when the particle is isolated and the effect of bringing two cyhnders together is to reduce the normal force. The normal force acts on the cyhnder in the direction of the electrolyte phase, • the normal force acting on constant surface charge cyhnders embedded in an uncharged 0/W interface shows markedly different dependencies on Ka and 6 from the corresponding constant potential interactions, • for constant surface potential interactions between cyhnders on an uncharged interface, the normal force varies with na up until Ka ^ 5 whereafter it remains constant; this is the same Ka limit found to apply in the Derjaguin analysis, • for constant surface charge interactions between cyhnders on an uncharged interface, the normal force varies with Ka up until na ^ 15 whereafter it remains constant, • an example case was given where it was shown that the electrostatic normal force per unit length acting on a constant surface potential (50 mV) silica cyhnder can Chapter 4. Results and Discussion 113 be ~ 10 times greater than the force of gravity acting on the cylinder; this example 4 was for a silica cylinder with a submicron radius embedded in an O/W interface with a dilute (10 ) electrolyte phase; for larger particles (> 30pm) and/or more -4 concentrated electrolytes the gravitational force quickly becomes the dominant normal force, • at small separations (/c [R-2a]=.25), the parallel force can be ~ 4 times greater than the normal force per unit length acting on the cylinders, • the force results indicate that the free energy change associated with the variation of the normal force (as a function of particle separation) will be negligible compared to the corresponding change in free energy associated with the parallel force, • examination of the parallel interaction force between constant potential cylinders at an uncharged O/W interface over a wide range of particle separations indicated that a strong, exponentially decaying force due to the overlapping double layers exists at small separations while a weaker, pseudo-algebraically decaying force, acting through the oil phase, exists over a much greater separation range, and • when the constant potential O/W interface affects the double layer interaction (depends on na and 8), the parallel interaction force depends strongly on the relative magnitudes and signs of ip and tpow- When tp^ is opposite in sign to V'o and of com0 parable magnitude, the parallel interaction force may be 10 times the magnitude of the force for the corresponding interaction on an uncharged O/W interface. 4.3 Close-Packed Two-Dimensional Swarm of Spheres at an O / W Interface This section presents the results obtained for Model III which was described in Section 2.4.3. Only interactions involving a constant potential sphere embedded in an uncharged 114 Chapter 4. Results and Discussion 0/W interface were considered for this model. Since the sphere is surrounded by an axisymmetric bounding envelope, there can be no net force in a direction parallel to the interface acting on the sphere. Thus, only the normal force was evaluated for this model. It was shown in Section 4.1 that the present free energy calculations yield results that are of the same order of magnitude when compared with literature values of the free energy of interaction between two spheres completely immersed in a single electrolyte. Thefreeenergy calculations were performed for this model to get an order of magnitude estimate of the free energy change associated with surrounding a single sphere on an 0/W interface with a dense two dimensional swarm of identical spheres. 4.3.1 N o r m a l Force A c t i n g on a Constant Potential Sphere E m b e d d e d in an Uncharged O / W Interface The reduced normal force (/") acting on a sphere in the unit cell model is calculated using equation (2.42). The results in this section are reported in terms of the scaled reduced force (/"*) which is obtained using equation (4.1); / " — / ^, • The normal 2ir a 2 forces reported in this section may be compared with the scaled reduced normal forces given for Model II (see Figs. 4.8 and D.5), but care should be taken when making the comparison. The most sensible way to compare the normal forces calculated for Models II and III would be to compare the force acting on the particle per unit of the particle perimeter in contact with the 0/W interface (hereafter called the 0/W perimeter). For cyhnders, the O/W perimeter is constant irrespective of the contact angle so we just consider the force per unit length values (measured in double layer thicknesses), divided by two (to account for the 'front' and 'back' sides), as reported in Section 4.2. For spheres, this 0/W perimeter varies with contact angle. The scaled reduced normal forces reported in this section are divided by the 0/W perimeter of sphere when 6 = 90°, i.e. 2TK,O.. Thus, 115 Chapter 4. Results and Discussion o - _L 10 15 20 25 30 ACQ Figure 4.16: Normal Force acting on an Isolated Constant Potential Sphere (with variable na) Embedded in an Uncharged 0/W Interface when 8 ^ 90° one divides /"* by sin 6? to get the normal force acting on the sphere per unit of the actual O/W perimeter. Again it is sensible to examine the effects that na and 8 have on the normal force of an isolated sphere before considering the effect of placing the sphere in a dense twodimensional structure. Fig. 4.16 shows the variation of the normal force for spheres as a function of na where the spheres are embedded (8 = 90°) in an uncharged O/W interface and the sphere surface in contact with the aqueous phase has tp — 2. Fig. 4.16 is 0 qualitatively similar Fig. 4.6 which plots the same variables for an isolated cylinder, i.e. as na increases the normal force asymptotically approaches a constant value. The sphere results, however, show a slower approach to the asymptote; the double layer decay along the O/W interface is more sensitive to na for spheres than cylinders. This behaviour is easily explained if one remembers the concept of the 'plane wall' result alluded to in Section 4.2.1, i.e. as na of the cylinder increases, the intersection of the double layer 116 Chapter 4. Results and Discussion Eo 30 60 90 120 Contact Angle 8 [deg] 150 180 30 60 90 120 150 180 Contact Angle 8 [deg] Figure 4.17: Normal Force acting on an Isolated Constant Potential Sphere (with variable 6) Embedded in an Uncharged 0/W Interface with the O/W interface will resemble that of a two planes. For a sphere, the 'plane wall' limit for Ka is greater than that for a cyhnder because the particle surface curves in two dimensions rather than one. For cyhnders the normal force per unit of O/W perimeter at large na is 1.08; for spheres the largest case calculated was for Ka = 28 where the force per unit of O/W perimeter is 1.10, i.e. a difference of < 2% compared to the cyhnder case. This is reassuring as one would expect the two results to be equivalent in the 'plane wall' limit. The plots in Fig. 4.17 show the variation of the normal force as a function of 6 (left plot) and the variation of the O/W perimeter-weighted normal force (right plot) for an isolated sphere with Ka = 10 and V'o = 2. The left plot in Fig. 4.17 shows markedly different behaviour from Fig. 4.7 which plots the same variables for the case of isolated cyhnders, but as expected, the plot of the O/W perimeter-weighted force shows similar behaviour to the cyhnder results. The separate contributions to the total normal force Chapter 4. Results and Discussion 117 Figure 4.18: Normal Force Acting on a Constant Potential Sphere Embedded in an Uncharged 0/W Interface given by equation (2.42) revealed that at low contact angles the osmotic term dominated and at high contact angles the (fj*r) 2 term dominates. This is the same trend that was noted and discussed for cylinders in Section 2.4.2. The important trend to note from Fig. 4.17 is that for 8 < 90° the normal force acting on the sphere is independent of 6. It is expected that this general rule should apply for na > 10. Fig. 4.18 contains plots of the normal force acting on a sphere (V'o = 2) for various sizes of the bounding envelope, i.e. the nearest distance from the sphere surface to the bounding envelope is measured by nd and this quantity has a minimum value of ~ 1.05 x na. Fig. D.14 in Appendix D contain analogous plots but with V'o = 4. Qualitatively, the results in Fig. 4.18 show the same trends as in Fig. 4.8, i.e. the force is affected the most when 8 — 90° and when na is small, the normal force is significantly affected even when 8 ^ 90°. The relative magnitudes of the normal force as a function of 8, however, show marked differences from the corresponding results for cylinders and 118 Chapter 4. Results and Discussion these differences can be attributed to the variation of the 0/W perimeter with 8. The changes in the force as a function of particle separation are greater in the case of the cell model than for cyhnders, e.g. for a sphere with na = 2 and 8 = 90° the normal force is reduced by a factor of ~ 85% when nd is decreased from 2.5 to .125, but for a cyhnder with na = 2 and 8 — 90° the normal force is only reduced by a factor of ~ 40% when K[R — 2a] is decreasedfrom5 to .25. The reason for the difference is the whole surface of the sphere is affected when the bounding envelope is shrunk, but only the 'front' side normal force is affected when two parallel cyhnders are brought together (the 'front' side normal force is reduced by ~ 85% for the interaction described above). For the same interaction as described above but with V'o = 4, see Fig. D.14, the normal force is reduced - by a factor of ~ 75%. All of the normal forces (/"") reported in this section varied .between .5 and 1.5. Assuming that a sphere composed of sihca is embedded in an O/W interface where: the electrolyte concentration is 10 -4 molar, the system is at room temperature (/c % -1 3 x 10 m), and the sphere has a radius of .3pm (so Ka = 10), then the gravitational force _8 acting on the sphere will be ~ 2.5 x 10 A . Recallfromthe beginning of this chapter _15 r that one unit of /" w 2 x 10~ N and that /" = 13 27r/eaV'L/*"- Then, if V'o = 2 (~ 50m V), the normal force due to electrostatic effects will be ~ 5 X 10 iV or about 2 X 10 times _11 4 the gravitational force! The normal force remains constant as the particle size increases (assuming K remains constant), so the ratio of the gravitational/electrostatic forces will increase with a , where a is the sphere radius. If we consider a sphere with a = 8pm and 3 assume all other factors remain constant then the electrostatic and gravitational forces will be almost equivalent. For larger particles and/or more concentrated electrolytes, the gravitational force quickly becomes the dominant normal force. 119 Chapter 4. Results and Discussion Table 4.3: Free Energy (G**) of Unit Cell with a Constant Potential Sphere (ip = 2) 0 e «d 45° 90° 135° 4.3.2 na 2 5 .5 3.22 7.45 1. 3.89 8.30 oo 4.33 8.62 .5 2.11 4.61 1. 2.48 5.03 oo 2.73 5.20 .5 .927 1.68 1. .983 1.69 CO 1.03 1.69 10 14.6 15.6 15.9 8.79 9.32 9.45 2.92 2.92 2.91 Free E n e r g y Calculations For a Constant Potential Sphere E m b e d d e d in an Uncharged O / W Interface This section presents the free energy results calculated for spheres with constant surface potentials in the unit cell model. The reduced electrostatic free energy of the system (solid, oil, water) is calculated using equation (2.37), (recall (Section 2.4.3) that in the non-polar phases the osmotic term is zero and the sign of the electrostatic field energy, W, is positive). The change in the free energy of the unit cell , AG, associated with forming the close-packed structure on the interface is calculated by taking the difference of the free energy of the unit cell when d is finite and the free energy when the particles are infinitely separated (d = co). The free energy of the unit cell was calculated for the same cases described in the last section and select values are presented in Table 4.3. Note that the energies are reported as G**, which was denned by equation (4.3). Conversion of the values in Table 4.3 into multiples of kT is straightforward, e.g. for 1.8fcr4ir/eaV>L **i G! s o f o r t h e c a s e o f K a = 1 0 > * d = - 5 a n K d 6 _ 1 = .OZpm, G ~ 1.8kTG' = 9 0 ° : G = 7 - 9 5 x 10 fcT. 3 = Chapter 4. Results and Discussion 120 Figure 4.19: Change in Free Energy for a Constant Potential Sphere Embedded in an Uncharged O/W Interface Fig. 4.19 shows the AG** values that were calculated as a function of nd and 8 for cases where na= 2,5,10 and i/'o = 2. The plots all show the expected trend, i.e. an increase in the repulsive energy as the size of the unit cell is reduced. As to be expected, the effects of the contact angle and na on AG" are very similar to the same effects seen for the force results of Model II. At 8 — 135° there is negligible change in AG** as a function of nd except at small separations for tea — 2. Noting the magnitudes of the G** (Table 4.3) and AG** (Fig. 4.19) values as a function of contact angle, it is apparent that the majority of the repulsive energy is contributed from the aqueous phase. Examination of the separate contributions to the total G " of the unit cell revealed that the aqueous phase contribution accounts for > 99% of the total. The results for the case of na — 10 in Fig. 4.19 are clearly erroneous since they show an attractive (negative) interaction energy for nd > 1.5 when 8 = 45° and 90°. Recall that in Section 4.1 it was shown that the free energy calculations for Model I resulted Chapter 4. Results and Discussion 121 in errors of up to 50% . The difference between the free energy of the system when the sphere is in a close-packed structure and when it is isolated is small, i.e. from Table 4.3 for na = 10 and 6 = 45° the difference between G** when d = 1. and oo is ~ 10%. Considering that these separate G** values may errors of up to 50% it is understandable, but disappointing, that substantial errors in AG** arise. As was discussed in Chapter 3, it is believed that a better grid design would reduce the degree of error; this subject will be addressed in Chapter 6. Nonetheless we may still use the data in Fig. 4.19 to give an order of magnitude estimate of the interaction energy. The greatest interaction energy available from Fig. 4.19 is for the case of na = 10, 6 = 45° and ud = 1: AG** = 1.30, so AG = 1.18 x 10 kT. 3 This is of the same order of magnitude as the prediction made by Levine et. al. [1] for the interaction energy for a sphere in a close-packed two dimensional structure on an O/W interface. The interaction energy due to the dipolar interaction through the oil phase cannot possibly be estimated using the present FEM. Since the double layer accounts for > 99% of the total energy of the system then it is essential that this quantity be accurately calculated for all sizes of the unit cell. If it is not calculated accurately as a function of d, then the relatively minute changes in G " due to the oil phase interaction will be masked by the error in the double layer energy. Considering that the differences in G" at large d and at oo will be < 1%, then G" itself must be known accurately to <C 1%. The chance of a better grid design, or any other numerical method for that matter, being able to yield this kind of accuracy is highly unlikely because of the very steep gradients that must be evaluated. 4.3.3 Summary The results of this section are summarized below in point form: apier 4. Results and Discussion 122 • the normal force acting on the constant surface potential sphere shows the same qualitative behaviour as a function of na, 6 and particle separation as the normal force acting on cyhnders with constant surface potentials (note that the respective normal forces must be 'perimeter-weighted' if this similar behaviour is to be borne out), • an example case was given where it was shown that the normal force acting on a submicron silica sphere, embedded in an O/W interface with a dilute ( 1 0 molar) -4 electrolyte phase, can be ~ 2 x 10 times greater than the gravitational force acting 4 on the sphere; for larger particles (> 8fim) and/or more concentrated electrolytes the gravitational force quickly becomes the dominant normal force, • the free energy calculations yielded results which were of the same order of magnitude as the predictions of Levine et. al. [1], and • the error associated with calculating the free energy by the present technique make it impossible to characterize the long range interaction acting through the oil phase. Chapter 5 Conclusions The conclusions drawn from the work presented in this thesis can be conveniently summarized for each of the models studied. The conclusions presented below are in point form; the reader is referred to Chapter 4 for details. 5.1 Model I 1. The Finite Element Method (FEM) developed for this thesis is capable of accurately simulating constant potential and constant charge interactions for a variety cases involving identical spherical particles in a single electrolyte. 2. The implicit assumption employed in the Derjaguin approximation for constant charge interactions (i.e. | £ = 0 inside the particle) will result in substantial errors (> 10%) as na and the particle separation decrease and as the dielectric constants of the interacting particles increase. 3. The hnear superposition approximation applies to constant charge interactions between identical spheres at approximately the same limit as it does for constant potential interactions (K[R — 2a] > 4) providing Ka > 1 or the particle dielectric constant is less than 80. 123 Chapter 5. 5.2 Conclusions 124 M o d e l II 1. The Derjaguin method, where the integration is performed only over the portion of the particle in contact with the aqueous phase, provides a good approximation to the parallel interaction force between constant surface potential cylinders embedded on an uncharged O/W interface when na > 5 and K[R — 2a] < 5. 2. For both constant potential and constant charge cylinders on an uncharged 0/W interface, a force in the direction normal to the 0/W interface exists which always acts to push the particle into the aqueous phase and for cases of dilute electrolytes (< 10 molar) and cylinders with submicron radii, the electrostatic normal force -4 can be substantially (~ 10 times) larger than the gravitational force acting on the 4 particle. 3. The normal force acting on isolated cylinders with constant surface potentials at an O/W interface is independent of na when na > 6, i.e. approximately the same Ka limit used when applying the Derjaguin method. 4. The major contribution to the free energy change associated with bringing two cylinders together from infinite separation will come from the integration of the parallel component of the interaction force. 5. Cylindrical particles embedded in an uncharged O/W interface exhibit a strong, short range double layer interaction which decays exponentially, and a weaker, long range interaction acting through the non-polar phase. 6. The effect of a charged 0/W interface on the parallel interaction force between constant surface potential particles embedded in the interface is highly unlikely to cause the free energy associated with the interaction to approach 10 fcT. 6 Chapter 5. 5.3 Conclusions 125 Model III 1. The normal force acting on constant potential spheres at an uncharged 0/W interface is very similar to the normal force acting on analogous cyhnders. 2. The double layer contributes more than 99% of the total electrostatic energy for the three phase unit cell model. 3. While the error incurred in calculating the electrostatic energy of the system is substantial, order of magnitude estimates of the interaction energy indicate that is well below 10 fcT. 6 Chapter 6 Recommendations This chapter presents recommendations for future work on this project. The recommendations are itemized below and the first three are then discussed in further detail for the benefit of future workers. 1. The force results for Model II should be integrated to yield the change in free energy for interacting clinders at an 0/W interface. 2. An effort should be made to either improve the grid design or obtain access to a 'state-of-the-art' commercial finite element analysis package so that the electrostatic free energy of the system can be directly and accurately calculated; this is essential for Model III. 3. For certain conditions (see Section 6.3) the electrostatic force acting on an interfacial particle in a direction normal to the interface may be sufficiently strong to distort the profile of the interface (which contravenes the initial assumption of a flat interface) ; if this occurs the interactions between interfacial particles would have to be solved iteratively. 4. The full three-dimensional model of the close-packed structure of particles at an 0/W interface could be attempted and the results compared with those of Model III. 126 Chapter 6. 127 Recommendations 5. An analytical treatment of the interaction of cylinders at an 0 / W interface could reveal the source of the pseudo-algebraic decay of the long range interaction predicted in this thesis. 6. The normal force calculated for a particle at an O/W interface is extremely dependent on the boundary conditions applied in the region where all three phases meet; an effort should be made to discover what is the most appropriate and realistic form of these boundary conditions. 6.1 Integrating the Interaction Force Integrating the interaction force has proven [28, 34, 35] to be the most accurate method to date of calculating the free energy change associated with double layer interactions between charged particles in a single electrolyte. The interaction energy between two particles at a given separation is obtained by integrating the interaction force acting between the particles from an infinite separation to the given separation. Previous authors [28, 34, 35] who have performed these calculations obtain the interaction force at discrete separation intervals and then integrate these results (no details of the integration method are given in any of the literature). The interaction forces are calculated up to a maximum, finite separation. At very large separations the integral may be evaluated analytically by employing the linear superposition approximation (LSA). Hoskin and Levine [28] (HL) give a lucid account of how to employ this technique for the case of two spheres fully immersed in a single electrolyte. Levine et. cd. [35] give a simple, more convenient expression for the interaction energy based on the LSA obviating the need for any integration at large separations. The next worker(s) may want to first test their integration procedure by integrating the force results for Model I; the resultant energies may be compared with previously Chapter 6. Recommendations 128 published data [28, 34, 35]. Evaluating the integral for Model II at large separations will introduce added complexities; the LSA applied to interacting cyhnders at an 0/W interface has not yet been derived. This obstacle may be circumvented if the force is evaluated for sufficiently large separations that at any greater separations the contributions to the integral are insignificant. The force was evaluated for a large separation range for three separate cases (see Section 4.2.1) and the numerical results are presented in Table D.l in Appendix D. The integration procedure could be applied to this data to get an idea of the magnitude of the long range interactions to the overall energy integral. Developing the integration procedure, testing it against literature values for Model I and then applying it to Model II may be suitable for a fourth year thesis project. 6.2 Improving the Finite Element Grid As was noted in Chapters 3 and 4, the present grid design is inadequate for directly calculating the free energy of the system. The main source of inaccuracy arises in the evaluation of the components of the electricfieldwhich are expressed as spatial derivatives of the electrostatic potential. A better grid design would improve the accuracy of the free energy calculations. In the past 10 or 15 years, grid generation theory ha£ exploded and really become a science unto itself. The number of present methods to choose from and the complexity of those methods can be quite daunting to the uninitiated numerical analyst. This section does not attempt to give a full review of current grid generation techniques, rather it serves to stimulate the next worker towards selecting an appropriate grid generation scheme. Before discussing possible alternatives, it is worthwhile to review the present grid design. The disadvantages of the present design are: Chapter 6. Recommendations 129 • the node distribution is not uniform around the circumference of the particle, especially close to the particle surface on the top and bottom (relative to the O / W interface) of the particle, • as a result of the non-uniform node distribution, the grid lines are badly skewed near the top and bottom of the particle, • the node spacing does not increase smoothly in a direction normal to the particle surface; the spacing can only be varied in the X and Y directions, • the element density varies unevenly from subregion to subregion in the grid, and • the element density far from the particle, in the water phase, is unnecessarily high. The present design, while riddled with faults, is not without some redeeming features: • a minor computing effort ( C P U time) is required to create the grid, • since we are dealing with a well-defined system (particles with circular cross-sections on a planar interface, separated by symmetry planes which are normal to the interface), variations of the geometric variables (i.e. particle size, contact angle, etc.) are performed easily, and • since the grid is constructed in a predefined manner for every variation of the same model, the assignment of node and element numbers follows a set pattern and thus the boundary elements required for the force calculations are easily referenced. The most promising pathways for establishing a better grid design are presented below. The majority of the information below was gleaned from a collection of review articles on grid generation, i.e. Eiseman [56], Turkel [57], Thompson et. cd. [58] and Thompson [59]. Most of the methods referenced in these reviews relate to computational fluid dynamics but these same methods are applicable to any field problem. Chapter 6. Recommendations 6.2.1 130 Algebraic G r i d Generation Algebraic grid generation schemes are the simplest to employ and they give the analyst explicit control over the node placements; the method used in this thesis is a crude version of an algebraic scheme. The basis of an algebraic scheme is that the nodal co-ordinates are obtained by interpolating between the boundaries of the domain. For severly distorted regions, for which interpolation between the global boundaries would result in unsatisfactory grids, it is necessary to interpolate between intermediate surfaces in the field; this is the concept of creating subregions within the overall domain. There is much choice when selecting the method of interpolating between surfaces; Eiseman [56] gives a good review of the more popular methods. The present grid design is unsatisfactory because its subregions have not been designed well. If an algebraic system is chosen to improve the current grid design, particluar care should be taken to: • make sure that the co-ordinate lines of the grid are normal and perpendicular to the constant potential (charge) surfaces when they are close to those surfaces, and • that the node spacing is variable in a direction normal to any constant potential (charge) surface in the model. The above requirements will be severely comprised in the regions where the O/W interface meets the particle, especially when the contact angle is very different from 90°. One saving grace is that for contact angles <C 90° the potential varies slowly in the region where the particle and interface meet so the electric field determination is less prone to error. A better approach can be used to calculate the nodes in the various subregions. Currently the nodes are established by calculating equations of straight lines within the 131 Chapter 6. Recommendations subregions in terms of the global co-ordinates and then placing nodes at increments along those hnes. An easier, and more elegant approach involves the concept of isoparametric mapping that was referred to in Chapter 2. Just as one can map a distorted boundary into a regular boundary, it is very easy to map a regular boundary into a distorted one. Recall that a global co-ordinate of a point in a subregion (element) is given in terms of the local co-ordinates by x = Sf Nx k k (6.1) where x axe the global x co-ordinates of the M boundary nodes of the subregion and k Nk are the basis functions (expressed in terms of the local co-ordinates) associated with those nodes (see chapter 2 for details). Mapping a regular shape into an irregular one is performed by: • assigning the boundary nodes in the global domain as the corresponding boundary nodes in the regular local domain, • interpolating between the boundaries in the regular, local domain to get the interior node positions, i.e. evaluate the basis functions at prescribed intervals of the local co-ordinates, and • calculate the global co-ordinates of the interior nodes using equation 6.3. This procedure causes the co-ordinate hnes near a boundary of the subregion to follow the contours of that boundary. A good description of this simple procedure is given by Zienkiewicz and Morgan [60]. The tricky part of using this procedure is defining the subregion boundaries. In the present design the subregion boundaries are always defined in the same way (see Section 3.1), which makes it very simple to use. More sophisticated arrangement of subregions will entail more computing effort but should reap better performance. An interactive Chapter 6. Recommendations 132 graphics package that allows the subregion boundaries to be digitized would obviate a lot of decision-making code that would be otherwise necessary; a good example of this approach is given by Smith et. cd. [61]. 6.2.2 Elliptic G r i d Generation Elliptic grid generation is one of a general family of grid generation schemes that involves solving partial differential equations to obtain the nodal co-ordinates. The domain is usually divided into quadrilateral subregions and the global co-ordinates (dependent variables) are obtained by solving partial differential equations (PDE) which relate the change in the global co-ordinates to the local co-ordinate (independent variables) space of the subregion. When the co-ordinate points (or slopes) are specified on the entire closed surface of the subregion the P D E must be elliptic. A n advantage of generating a grid in this manner is that the grid P D E and its boundary conditions may be manipulated to mirror the physics of the problem being approximated and hence the resulting co-ordinate lines will tend to be similar to the physical field lines. The success of this manipulation requires knowledge of the solution, so the advantage is usually only gained after an iterative process of solving the physics and altering the grid accordingly. Obvious disadvantages of elliptic schemes is that they require more effort (both in software development and execution time) and the analyst has less control over the resulting grid than if an algebraic scheme is used. Thompson et. al. [62] give a good review of the basics of elliptic methods. Applications of Laplace and Poisson-like PDE's are examined. The sink/source term in the Poisson-like equations is called a 'forcing function' which serves to control the co-ordinate line spacings. The forcing function can take on many different forms; it may be attractive in some areas of the grid causing points to be concentrated in that region or repulsive in another region causing the points to be widely spaced. Three separate, published codes Chapter 6. Recommendations 133 are referenced by Thompson [59] which are based on elliptic methods. Examples from the literature are referenced where these codes have been successfully implemented. A major problem with introducing sophisticated forcing functions, which usually result in nonlinear PDE's, is that the grid must be calculated iteratively and thus a good idea of the final grid solution must be available for convergence to occur. It is sometimes recommended that the grid is first formed with a linear set of equations and then the nonlinearities are gradually introduced on subsequent iterations until a satisfactory grid is obtained. This procedure is very time consuming and perhaps not practical for the type of work considered in this thesis. 6.2.3 Adaptive Grids The most active area of grid generation research at present is the development of dynamically adaptive grids. In these systems the grid points move on the grid in response to the physical solution being performed on the grid. A driving force exists that causes the points to congregate in such a fashion that the accuracy of the solution is maximized. A typical driving force used is based on the principle that the solution should vary by the same amount between any pair of points. In a finite difference solution, this causes the second and higher order differences of the field variable with respect to the local co-ordinates to vanish. These adaptive methods are almost all based on variational principles and thus they are not immediately accessible to the average chemical engineer. One interesting (and more palatable) branch of adaptive grid research is that of moving finite elements. In this scheme, the grid points are made additional dependent variables in the Galerkin formulation of the problem [63]. The grid point locations are then solved as part of the finite element solution. These adaptive methods are very complex and unless one is very familiar with grid generation techniques in general, one is hkely to be overwhelmed by the amount of applied Chapter 6. Recommendations 134 mathematics required to develop such schemes. The alternative is to buy the expertise. Reading the following words of J.F. Thompson [59], one wonders whether this type of grid generation software will become a common (black box) tool for engineering applications in the future: "Dynamically adaptive grids coupled with the physical solution evolving thereon can be expected to be the most promising area of (numerical) research in the coming years. It has been noted that when the grid is properly positioned to resolve the solution gradients, most solution algorithms work well. ...the ultimate answer to numerical solution of partial differential equations . . . may well be dynamically adaptive grids, rather than more eleaborate difference representations and solution methods". 6.2.4 PLTMG PLTMG is a software package available at the University of British Columbia Computing Centre [64] that is used to solve boundary value problems. The solution is obtained using the finite element method based on C° linear triangular finite elements. The user must provide a crude triangulation of the solution domain and further triangulation of the domain may be performed explicitly by the user or adaptively by PLTMG. This package was not used in the present application because inconsistencies were reported [65] by several users that could not be resolved. These problems may be corrected in the future so it is advised to keep abreast of the status of this potentially helpful and powerful! package. 135 Chapter 6. Recommendations 6.3 Distortion of the Interface Recall (Section 2.4) that the analyses presented in this thesis assume that the flat interfacial profile is not distorted in any way by the presence of interfacial particles. By conducting a force balance on the interfacial particle, where only the affects of gravity and interfacial tension are considered, it is possible to show [66, 67] that an interfacial particle will not distort the interface providing the Bond number is very small. For spheres and cylinders the Bond number is defined as, B = { p l ~ p 2 ) 9 a (6.2) 2 7 where pi and pi are the densities of the two fluid phases (p\ > p ), g is the acceleration 2 due to gravity, a is the radius of curvature of the particle and 7 is the interfacial tension between the two fluids. Substituting typical numerical values for the experimental conditions of Bowen and Levine [1] (a = .3pm, 7 « .05N/m, Ap 100kg/m ) one obtains 3 Bond numbers on the order of ~ 10~ (low enough to cause negligible distortion of the 9 interface providing the particle density is not many orders of magnitude different than Pi). In Chapter 4 it was revealed that for small na (a = .3pm with an electrolyte concentration of 10 -4 molar) the electrostatic force can be ~ 10 greater than the gravitational 4 force acting on the particle. Just like the force due to interfacial tension acting on the particle, the electrostatic force acts per unit length of the particle perimeter in contact with the 0 / W interface. When the electrostatic force is added to the overall force balance equation for the interfacial particle (see equation (85) in ref. [66]) and the equation is made dimensionless, a factor, L, representing the ratio of the electrostatic force to interfacial tension exists Chapter 6. Recommendations 136 L = F (6.3) h elec where F i e ec is the electrostatic normal force per unit length of the 0/W perimeter. If L is added onto equation (86) in ref. [66] it can be seen that as long as L < 10~ then 2 the distortion of the interface will be very small. If we assume 7 « .05N/m then the upper limit of F t to cause distortion is bxlO~ N/m. 4 e ec Recall (Section 4.2) that for a .Zu-m cyhnder (I/J = 2) in a IO molar electrolyte solution, F i « 2 x 10~ N per unit 12 -4 0 e ec double layer thickness or ~ 6.7 x 10~ N/m. 5 From Fig. D.5 in Appendix D one finds that Feiec ~ 3 x 10~ N/m for the same cyhnder as above but with \\) = 4(« lOOmV). This 4 0 result is very close to the upper limit described above. Similar results are obtainable for spheres. It may be possible, given large surface potentials on the particles, that the electrostatic normal force may cause marginal distortion of the interface. Even marginal distortion may be important due to possible capillary effects so some thought should be given to solving the electrostatic interactions simultaneously with the capillary interactions. This would be an iterative procedure and is certainly beyond the scope of this thesis. Bibliography [1] S.Levine, B.D. Bowen and S.J. Partridge, Report for Agreement 451 with AOSTRA (Alberta Oil Sands and Research Authority, Edmonton, Can.), October, 1986. S.U. Pickering, J. Chem. Soc, 91(1907)2001. J.A. Kitchener and P.R. Musselwhite, "The Theory of Stability of Emulsions" in Emulsion Science (P.Sherman Ed.), Academic Press, London(1968). J.H. Schulman and J. Leja, Trans. Faraday Soc, 50(1954)598. CW. Bowman, Proc 7th World Petroleum Congress, Mexico, 3(1967)188. P. Finkle, H.D. Draper and J.H. Hildebrand, J. Am. Chem. Soc, 45(1923)2780. A.H. Taubman and A.F. Koretskii, in "The Use of Surfactants in the Petroleum Industry", (P.A. Rebinder Ed.), Consultants Bureau, New York (1965), p.188. P. Becher, "Emulsions: Theory and Practice", 2nd. Ed., Reinhold, New York (1965). S. Levine and E. Sanford, Can. J. Chem. Eng., 62(1985)258. P. Pieranski, Phys. Rev. Letters, 45(1980)569. V.B. Menon and D.T. Wasan, Colloids Surfaces, 19(1986)89. K.P. Oza and S.G. Frank, J. Dispersion Sci. Tech., 7(1986)89. S. Levine, B.D. Bowen and S.J. Partridge, Colloids Surfaces, to be published. S. Levine, B.D. Bowen and S.J. Partridge, Colloids Surfaces, to be published. R.D. Void and M.J. Void, "Colloid and Interface Chemistry", Addison-Wesley, Don Mills, Ont.(1983). [16] D.J. Shaw, "Introduction to Colloid and Surface Chemistry", 3rd. Ed., Butterworths, Toronto (1980). [17] A.L. Loeb, J.Th.G. Overbeek and P.H. Wiersema, "The Electrical Double Layer Around a Spherical Colloid Particle", MIT Press, Cambridge, Mass. (1961). 137 Bibliography 138 [18] J.Th.G. Overbeek in "Colloid Science", Vol. I , (H.R. Kruyt Ed.), Elsevier, Amsterdam (1952). [19] E.J.W. Verwey and J.Th.G. Overbeek, "Theory of the Stability of Lyophobic Colloids", Elsevier, Amsterdam (1948). [20] H. Miiller, KoUoidchem. Beih., 26(1928)257. [21] N.E. Hoskin, Trans. Faraday Soc, 49(1953)1471. [22] S.S. Dukhin and B.V. Derjaguin in "Surface and Colloid Science", (E. Matijevic Ed.), Vol. 7, Wiley, New York (1974). [23] J.Y. Parlange, J. Chem. Phys., 57(1972)376. [24] B. Abraham-Schrauner, J. Colloid Interface Sci., 44(1973)79. [25] L.R. White, J. Chem., Soc, Faraday Trans. II, 73(1977)577. [26] S.L. Brenner and R.E. Roberts, J. Chem. Phys., 77(1973)2367. [27] L.R. White, H. Ohsima and T.W. Healy, J. Colloid Interface Sci., 90(1982)17. [28] N.E. Hoskin and S. Levine, Phil. Trans. Roy. Soc. London, A248(1956)449. [29] S. Levine, Proc Phys. Soc, 64A(1951)781. [30] J.M. Crowley,"Fundamentals of Applied Electrostatics", Wiley, Toronto (1986). [31] S. Levine and A. Suddaby, Proc Phys. Soc, 64A(1951)287. [32] S. Levine and A. Suddaby, Proc Phys. Soc, 65A(1952)405. [33] B.V. Derjaguin, Trans. Faraday Soc, 36(1940)203. [34] L.N. McCartney and S. Levine, J. Colloid Interface Sci., 30(1969)345. [35] G.M. Bell, S. Levine and L.N. McCartney, J. Colloid Interface Sci., 33(1970)335. [36] R. Hogg, T.W. Healy and D.W. Fuerstenau, Trans. Faraday Soc, 62(1966)1638. [37] G.R. Wiese and T.W. Healy, Trans. Faraday Soc, 66(1970)490. [38] G.M. Bell and G.C. Peterson, J. Colloid Interface Sci., 41(1972)542. [39] J.Th.G. Overbeek, J. Chem. Soc, Faraday Trans. I,'84(1988)3079. [40] E.P. Honig and P.M. Mul, J. Colloid Interface Sci.,36(1971)258. Bibliography 139 [41] J.R. Philip, J. Chem. Phys., 52(1970)1387. [42] J.E. Ledbetter, T.L. Croxton and D.A. McQuarrie, Canad. J. Chem., 59(1981)1860. [43] B.K.C. Chan and D.Y.C. Chan, J. Colloid Interface Sci., 92(1983)281. [44] H. Sasaki, E. Matijevic and E.J. Barouch, J. Colloid Interface Sci., 76(1980)319. [45] A.L. Nichols and L.R. Pratt, J. Chem. Phys., 77(1982)1070. [46] A.J. Hurd, J. Phys. A., 18(1985)11055. [47] F.H. StiUinger, J. Chem. Phys., 35(1961)1584. [48] B. Jancovici, J. Statistical Phys., 28(1982)43. [49] J. Happel and H. Brenner," Low Reynolds Number Hydrodynamics", Prentice-Hall, Toronto (1965). [50] O.C. Zienkiewicz, "The Finite Element Method", 3rd. Ed., McGraw-Hill, Toronto (1977). [51] FX. Stasa, "Applied Finite Element Analysis for Engineers", Holt, Rinehart and Winston, Toronto (1985). [52] M.H. Protter and C B . Morrey, Jr., "Modern Mathematical Analysis", AddisonWesley, Don Mills (1964). [53] T. Nicol, 'UBC Sparspak', Computing Centre, UBC, Vancouver (1982). [54] A. George and J. W-H. Liu, ' Computer Solution of Large Sparse Symmetric Positive Definite Systems of Linear Equations', Prentice-Hall (1981). [55] F.P. Incropera and D.P. DeWitt, "Fundamentals of Heat and Mass Transfer", 2nd. Ed., J. Wiley and Sons, Toronto (1985). [56] P.R. Eiseman, Ann. Rev. Fluid Mech., 17(1985)487. [57] Z. Turkel, Computers and Fluids, 11(1983)121. [58] J.F. Thompson, Z.U.A. Warsi and CW. Mastin, J. Comp. Phys., 47(1982)1. [59] J.F. Thompson, AIAA Journal, 22(1984)1505. [60] O.C Zienkiewicz and K. Morgan, "Finite Elements and Approximations", J. Wiley and Sons, Toronto (1983). Bibliography 140 [61] R.E. Smith, R.A. Kudlinski and E.L. Everton, "A Grid Spacing Control Technique for Algebraic Grid Generation Methods", AIAA paper 82-0226, Orlando, Fla., Jan. (1982). [62] J.F. Thompson, Z.U.A. Warsi and CW. Mastin, "Numerical Grid Generation", North-Holland, New York (1985). [63] K. Miller and R.N. Miller, SIAM J. Numer. Anal, 18(1981)1019. [64] T. Nicol, 'UBC PLTMG', Computing Centre, UBC, Vancouver (1982). [65] Personal Communication with T. Nicol, Computing Centre, UBC, Vancouver. [66] H.M. Princen in "Surface and Colloid Science", (E. Matijevic Ed.), Vol. 2, Wiley, New York (1969). [67] D.Y.C. Chan, J.D. Henry, Jr. and L.R. White, J. Colloid Interface Sci., 79(1981)410. [68] Personal Communication with S. Levine, Dept. Chemical Engineering, UBC, Vancouver. Appendix A Derivation of Force Equations for Three Phase System The derivation of the expressions for the force components normal and parallel to the oil/water ( 0 / W ) interface in the two cyhnder model are presented. This derivation is due to Levine [68]. Refer to Fig. A . l for the definitions given below. Consider the 'half'-volumes V and Vo of the water and oil separately. The boundaries of V^, and V w are S w a and S respectively, where D S -- Sum -f- Surr -\- Sow H* Syjp w S = 0 (A.l) Som "f Sor -f Sow "i" Sop (A.2) and the particle boundary, S , is p Sp — Swp -t- Sop (A.3) Using the method of Levine and Hoskin [27] we have on S w I JSw {^[(E 47T -n)n-\E n\2 AUn}dS = 0 I ^{(E-n)n-\E n)dS 4TT 2 = 0 2 (A.4) . and on S 0 JSo 2 (A.5) where n is the unit normal directed outwards from V or V and all the above quantities w 0 are as previously defined in the text. If we separate out the integrals over the particle surface we can write — 47T I —[(E-n)n-]-E n}dS+[ Jsw Air 2 / [(E • n)n - \E n]dS 2 JSp 2 2 —[(E-n)n-\E n]dS-l Air 2 + / ATlndS = JSpw ATlndS 2 Js'o 141 JSw (A.6) Appendix A. Derivation of Force Equations for Three Phase System 142 Figure A.l: Two Cylinders at an O/W Interface where S = S — S^ and S = S — Sop. w w 0 0 The left-hand side of equation (A.6) is the force acting on the cylinder. The surface integrals on the right-hand side are carried out over S — S^, and S — S^. The integrals w 0 over Sur and 5^. vanish as these surfaces extend to infinity. The integrals over the median plane surfaces 5 ^ and Sem yield the component of force acting parallel to the line of centres and have already been considered [28] and the expressions for this case are given in Chapter 2. The additional terms arise from the integrals over the 0/W interface and we now proceed tofindthe component normal to the interface and the component normal to thefineof centres and parallel" to the interface. First consider the component normal to the O/W interface starting with the waterside integrals [{E-n)n-\E n\dS 2 47T JSow 2 - / JSow AUndS (A.7) Appendix A. Derivation of Force Equations for Three Phase System 143 Figure A.2: Co-ordinate System Introduce Cartesian co-ordinates x, y, z where z is normal to the interface which lies in the xy plane Since n is outwards from V then n = —k. The k component of equation w (A.7) is ^ 47T f JSow [(E •fc) - \E )dS 2 2 2 + j ARdS JSow where ty is the potential function in the water phase and the (^"-) term is assumed to w be zero since we are ignoring the end effects of the cylinders. The contribution to the normal force from the oil phase is obtained in the same manner except that n = +k and there is no osmotic term where \? is the potential function in the nonpolar oil phase. The total component in n Appendix A. Derivation of Force Equations for Three Phase System the normal direction is given as the sum of equations (A.9) and (A.8), L Since $ w & = Q r - ( ^ ) ' l + £[<£r - ( ^ ) ' l at the interface then it follows that If the interface is uncharged, then (ew^^ = o ^ ) 6 a + = a n ( 144 i.e. Anns (A.10) at the interface as well. l hence equation (A. 10) may be reduced to D ^ i t ^ z ^ ^ ^ 3 ( A n ) which may also be expressed as U ^ 8 - t ? + e t { i £ ) 1 ] + A n } d S ( A 1 2 ) Just as with finding the normal component of force, for the parallel component (xdirection) we consider the contributions from the water and the oil sides to arrive at sLl-(^^) + ^X^)KS (A.13) which is equal to zero when there is either no charge on the interface or a constant potential on the interface. When there is a constant charge on the interface equation (A. 10) becomes where o~ is the charge on the interface. Appendix B Quadratic Order Triangular Elements The element co-ordinate system used for this work , the basis function definitions and their first derivatives are supplied here as a reference. The boundary integral terms for imposing a flux condition are also given. B.l Basis Functions and T h e i r Derivatives Isoparametric, quadratic order triangular elements are defined by six nodes and two independent area co-ordinates; Li and L . A third dependent area co-ordinate is denned 2 by L = 1 — Li — L . Fig. B.l shows the node numbering scheme and the values of the 3 2 area co-ordinates within the element. The basis functions, N(i), associated with the i th JV(1) = Li(2Li - 1) N{2) = L (2L - 1) N(3) = L (2L - 1) 2 3 2 3 nodes are N(4) = AL L X 2 N(b) = AL L 3 N(6) = 4LiL 3 2 (B.l) The basis function derivatives, expressed in terms of the independent area co-ordinates 145 Appendix B. Quadratic Order Triangular Elements 146 Figure B.l: Quadratic Order Triangular Elements Li and L , are 2 8N(i) dLi 8N(i) { 8N(i) dL 2 8N(i) 8Li ) L 2 ' L 3 8L V 8N(i) K ) L l M 3 8N(i) 8L ^ >1J 2 Z 8L v (B.2) Z where the subscripts on the bracketed terms indicate those quantities held constant during the differentiation. The resulting expressions for the derivative terms are then dN^/dL-, = 4Li - 1 dN(2)/dLi = 0 dN{2)/dLi = 1 - 4 L 3 8N(4:)/dLi = 4 L 2 8N(5)/8Li 8N(6)/8Li = -4L 2 = 4(L - L ) 3 x 2 = 0 8N{l)/8L Appendix B. Quadratic Order Triangular Elements 8N{2)/dL = AL - 1 2 2 8N(3)/8L = 1-4L 2 8N{4)/8L 8N(5)/8L 2 3 =4Ii 2 = A{L 3 8N(6)/8L L) 2 = -AU 2 B.2 147 (B.3) Constant Charge B o u n d a r y Condition From Gauss' law, the charge on a surface is given as the net charge flux leaving the surface, which corresponds to the sum of the flux contributions to the load vector, as was discussed in Chapter 2. In the FEM the fluxes are expressed in terms of boundary integrals according to Green's theorem. In Cartesian co-ordinates these integrals have the form fl= J^**N tdC T where cr" = y/^ii^-) — e (B.4) * the constant reduced surface charge given in terms s of the reduced quantities defined in Chapter 2. The integration is performed over the perimeter, C , bounding the element of thickness t. For the axisymmetric case we have e The dC term in the above equations is an infinitesmal arc length of the perimeter e which from basic calculus is yj{dx) + (dy) for the Cartesian system. 2 2 All the parti- cles considered in this thesis have circular cross-sections which are divided into equal segments. Thus the integration interval, dC , is constant for every element around the e perimeter of the particle and it is given by £ 2 and by T E A for models I and III (half particles) for model II (whole particle perimeter) where NEA is the number of radial elements around the particle. Appendix B. Quadratic Order Triangular Elements 148 The grid design is such that only the 1-4-2 leg of the particle boundary elements needs to be considered. Thus the basis functions for the 1,4,2 nodes can be expressed in terms of L (since L is zero and X x 3 = 2 1 ~ -^I) and the above integrals can be evaluated exactly, i.e. one merely integrates the basis function associated with each node with respect to the area co-ordinate L\. The contributions to the load vector elements in the Cartesian co-ordinate system are / ( l ) = <r*dC /6 e e = <T*dC /6 f (2) e e / (4) = 2<r*dC /3 e e (B.6) (B.7) In the axisymmetric case, products of the basis functions must be integrated since r is contained in the integral. The resulting expressions are / ( l ) = <T'dC [2r(l) - l/2r(2) +r(4)]/15 e e / (2) = o-'dC [-l/2r(l) e e + 2f(2) +r'(4)]/15 / (4) = <r*dC [r(l) +r'(2) + 8r(4)]/15 e e (B.8) (B.9) where r(i) is the r co-ordinate of the i th node. The factor of 27T has been left out of the above expressions since it was left out of the calculations for the stiffness matrix. Appendix C FEM Software This appendix contains the source code (written in Fortran 77) used to calculate the potential distribution, and the parallel and normal forces for Model II. The subroutine ENERGY used to calculate the free energy of the double layers and electric fields for Models I and III is also provided. Finally, a typical input file is given. C.l P O W I . C Y L : M o d e l II c C C C C C C C C C C C C C C C POWI.CYL (PARTICLE ON OIL/WATER INTERFACE) CYL: TWO CYLINDRICAL PARTICLES WITH SYMMETRY PLANE. MICHAEL P. LYNE MAY/89 version : 3 Solves f o r the e l e c t r i c p o t e n t i a l d i s t r i b u t i o n around a charged c y l i n d r i c a l p a r t i c l e on a f l a t , charged i n t e r f a c e . The PB equation i s solved i n the aqueous phase while Laplace's i s solved i n the other two phases. Only quadratic order t r i a n g u l a r elements are used. $RUN -LOAD+GRDC.E 5=(input file) ( g r i d software i n GRDC.E) MAX ELEMENTS:2500 MAX NODES:5000 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C SUBROUTINES CALLED cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc GRIDCY c creates the f i n i t e element g r i d . LOAD c i n i t i a l i z e s load v e c t o r (constant charge i n t e r a c t i o n ) . SPSSLV c SPARSPAK i n i t i a l i z a t i o n , system subroutine, INIJIJ c inputs matrix s t r u c t u r e t o SPARSPAK. 0RDRA5 c reorders sparse matrix i n SPARSPAK. SETMAT c f i l l s elemental submatrices. IMAT5 c inputs submatrices t o SPARSPAK. IAIJ5 c adds values t o s t i f f n e s s matrix i n SPARSPAK. IRHS c inputs right-hand side l o a d v e c t o r t o SPARSPAK. S0LV5 c solves the set of equations i n SPARSPAK. BNDINT c performs boundary i n t e g r a l s f o r f o r c e c a l c u l a t i o n s . cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C VARIABLE LIST cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Arrays: C AJAC Jacobian matrix. 149 Appendix C. FEM Software C AJIN C DBS.DBT C DBX.DBZ C EM C IBE C II,JJ C KEY C ND1,ND2,ND3 C NOTE C OLD C RES C SBM C SBOLD C SEXP C SOR C VAL C XAN C ZND C inverse Jacobian matrix. b a s i s func. d e r i v . w.r.t. l o c a l co-ordinates. b a s i s func. d e r i v . w.r.t. g l o b a l co-ordinates. d i e l e c t r i c constants. element i n d i c a t o r s f o r boundary i n t e g r a l s . g l o b a l node numbers f o r elemental s t i f f n e s s matrices. index array f o r elements and nodes. storage key arrays f o r boundary nodes. phase i n d i c a t i n g array. s o l u t i o n vector from previous i t e r a t i o n . right-hand side vector s u p p l i e d t o SPARSPAK. elemental s t i f f n e s s matrices. elemental s t i f f n e s s matrices from previous i t e r a t i o n . c o e f f i c i e n t s f o r expansion of Boltzmann term. new guess values f o r p o t e n t i a l s . lower l e f t t r i a n g l e of elemental s t i f f n e s s matrix. s o l u t i o n vector. contains g l o b a l co-ordinates of nodes and boundary condition indicator. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Input v a r i a b l e s ( i n input o r d e r ) : C NEA elements r a d i a l l y around p a r t i c l e . C NEI elements from p a r t i c l e t o i n f i n i t y on backside. C NED elements from p a r t i c l e t o median plane. C NWA elements from p a r t i c l e t o i n f i n i t y i n H20 phase. C NOB elements from p a r t i c l e t o i n f i n i t y i n o i l phase. C MAXIT maximum number of i t e r a t i o n s . C NALP contact angle i n d i c a t o r . C RA p a r t i c l e r a d i u s . C RB r a d i u s of inner p i e r e g i o n of p a r t i c l e . C ZO depth of o i l phase. C XI extent of i n f i n i t e boundary i n the x - d i r e c t i o n . C ORZ.ORX co-ords. of centre C ZM.CP -,+ valences of the b i n a r y e l e c t r o l y t e . C EH d i e l e c t r i c constants of three phases. C PSI p a r t i c l e surface p o t e n t i a l at i n f i n i t e separation. C ERR convergence c r i t e r i o n ('/,). C ZW depth of water phase. C XD2 distance between p a r t i c l e and median plane. C SKL r e l a x a t i o n parameter (alpha). C FZWA.FZWP.FZOP.FZOB - not used i n t h i s v e r s i o n . C FCZ node spacing f a c t o r i n the z - d i r e c t i o n . C FCXI node spacing f a c t o r i n the x - d i r e c t i o n , backside. C FCXD node spacing f a c t o r i n the x - d i r e c t i o n , f r o n t s i d e . C SOW o/w i n t e r f a c e p o t e n t i a l , (-99. means no charge cond.) C NQ constant charge, constant p o t e n t i a l f l a g f o r p a r t i c l e , C O=constant p o t e n t i a l , l=constant charge. C PINF p a r t i c l e surface p o t l . at i n f i n i t e separation f o r C constant charge case. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Remaining v a r i a b l e s : C ALP immersion angle of p a r t i c l e at o/w i n t e r f a c e . C EMXO l a r g e s t e r r o r fromn previous i t e r a t i o n . C FLG ZND(i,3) value f o r a non-boundary node. C FLUX charge term used i n RHS vector f o r const, charge case. C KSURF surface i n d i c a t o r f o r subroutine BNDINT. C LFG i t e r a t i o n f l a g ; 0 = f i r s t i t e r ' n , l=subsequent i t e r ' n s . C NEL elements i n g r i d . C NEO elements r a d i a l l y around p a r t i c l e i n contact with o i l . C NEW elements r a d i a l l y around p a r t i c l e i n contact with H20. C NIT- i t e r a t i o n counter. C NOD nodes i n the g r i d . C SQL (CP+ZM)/2ZM, lambda**2. 151 Appendix C. FEM Software C C SR ZP scaling factor f o r a l l distances. depth of o/w i n t e r f a c e r e l a t i v e t o p a r t i c l e centre. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc IMPLICIT REAL*8 (A-H.O-Z) DIMENSION XAN(5000),VAL(21),11(21),JJ(21) DIMENSION 0LD(5000),ND1(100),ND2(100),ND3(100) DIMENSION IBE(100,2),SB0LD(2500,21) CHARACTER*! LIST COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), fiSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),NOTE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD 0. N0BI.N0BD COMMON /BLK2/Z0,ZW,XI,RA,RB,0RX,0RZ,ALP,ZP,PI,PSI,CP,S0W COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2),AJIN(2,2) 1, DETJ,EM(3) COMMON /BLK6/XD2,FZWA,FZWP,FZOP,FZOB,FCZ,FCXI,FCXD COMMON /FVEC/RHS(5000) COMMON /SPKERR/ I ERR DATA B0L,EL,AVN,TEMP,C0NC/1.380662D-16,4.80336D-10, ! 6.022045D23.298.15D0.1.D-2/ DATA LIST,BIG/'*',1.D25/ C C Read i n data. C READ(5,1030) NEA,NEI,NED,NWA,NOB,MAXIT,NALP READ(5,1010) RA.RB.ZO.XI READ(5,1010) ORX.ORZ.CP.ZM READ(5,1010) EM.PSI READ(5,1010) ERR,ZW,XD2,SKL READ(5,1010) FZWA,FZWP,FZOP,FZOB READ(5,1010) FCZ,FCXI,FCXD,SOW READ(5,1011) NQ.PINF PI=DATAN(1.D0)*4.D0 1010 F0RMAT(4D8.2) 1011 F0RMAT(1X,I2,D8.2) 1030 F0RMAT(7I5) C C Adjust dimensions f o r v a r i a b l e d i e l e c t r i c constant. C ZW=RA+ZW XI=RA+XI ZO=RA+ZO SR=DSQRT(EM(2)) RA=SR*RA RB=RB*SR Z0=Z0*SR ZW=ZW*SR XI=XI*SR XD2=XD2*SR 0RX=0RX*SR ORZ=ORZ*SR C C Preliminary calculations. C NEA2=NEA/2 NE0=NEA2-NALP+1 NEW=NEA2-NE0 SqL=(CP+ZM)/(2.DO*ZM) RLM=DSQRT(SQL) C1=B0L*TEMP*299.8D0/EL FM=Cl**2*EM(2)/4.DO/PI ALP=(NALP-1)*PI/NEA2 Appendix C. FEM Software ZP=RA*DCOS(ALP)+ORZ IF(DABS(ZP).LT.1.D-14) ZP=O.DO C C C a l c u l a t e constants i n s e r i e s expansion. C FAC=1.DO ZS=1.D0 DO 10 1=1,20 FAC=FAC*ZS SEXP(I)=(ZM**I-(-CP)**I)/FAC ZS=ZS+1.D0 10 CONTINUE C C Create the f i n i t e element g r i d . C FLG=-8888.D0 CALL GRIDCY PRINT*,NEL,NOD IF(N0D.GT.5000.0R.NEL.GT.2500) THEN PRINT *,'**•* EXCEEDED MEMORY LIMIT ****' PRINT *, ' NODES = ',N0D,' ELEMENTS = ' ,NEL STOP END IF C C I f a constant charge - boundary c o n d i t i o n i s chosen f o r the p a r t i c l e C surface then c o l l e c t the boundary elements on the p a r t i c l e and C i n i t i a l i z e the RHS l o a d vector with the proper terms. C IF(NQ.EQ.O) GOTO 29 C C Gather p a r t i c l e boundary elements. C Section 6 NST=3*NEA+NWAI+2*NEI*NE0+NWAD+2*NE0*NED+1 NE=NST LF=3 DO 21 1=1,NEW IBE(I,1)=NE IBE(I,2)=LF NE=NE+1 21 CONTINUE C Section 2 NE=NST+2*NED*NEW+N0BD DO 22 J=1,NEW IBE(I,1)=NE IBE(I,2)=LF NE=NE+1 1=1+1 22 CONTINUE IBE(I,l)=-i DO 23, 1=1,NOD RHS(I)=0.D0 23 CONTINUE CALL LOAD(IBE) C C C a l c u l a t e the reduced charge d e n s i t y u s i n g WHITE'S (81) C method. C BET=DBESK0(RA/SR)/DBESK1(RA/SR) QS=2.D0*DSINH(PSI/2.D0)*DSQRT(1.D0+(1.D0/BET**2-1.D0) $/(C0SH(PSI/4.D0))**2) FLUX=QS*PI*RA/NEA/6.DO/2.D0*EM(2) C C I n i t i a l i z e SPARSPAK Appendix C. FEM Software c 29 M=4*NEL CALL SPSSLV(XAN,VAL,RHS,ICL,II,JJ,IR,N0D,M,1) C C Input matrix s t r u c t u r e . C DO 30 NE=1,NEL K=l DO 25 1=1,6 NI=KEY(NE,I) DO 25 J=1,I JJ(K)=KEY(NE,J) II(K)=NI K=K+1 25 CONTINUE CALL INIJIJ(21,II,JJ) IF(IERR.EQ.O) GOTO 30 WRITE(6,LIST) IERR,' I N I J I J ' PRINT *,' NODES = ',NOD,' ELEMENTS = \NEL STOP 30 CONTINUE C C Zero Submatrices. C DO 35 NE=1,NEL DO 35 1=1,21 SBM(NE,I)=0.D0 35 CONTINUE C C ORDER C CALL 0RDRA5 IF(IERR.EQ.O) GOTO 40 MRITE(6,LIST) IERR,' 0RDRA5' STOP C C Begin convergence procedure. C 40 LFG=0 NIT=0 EMX0=1.D3 C C F i l l submatrices. C 45 CONTINUE CALL SETMAT(LFG) C C Input submatrices. C DO 60 NE=1,NEL DO 50 IN=1,21 VAL(IN)=SBM(NE,IN) IF(N0TE(NE).NE.2.AND.LFG.EQ.1) GOTO 50 IF(N0TE(NE).Eq.2)SBM(NE,IN)=0.D0 50 CONTINUE 51 K=l DO 55 1=1,6 NI=KEY(NE,I) DO 55 J=1,I NJ=KEY(NE,J) II(K)=MAX(NI,NJ) Appendix C. FEM Software 55 JJ(K)=MIN(NI,NJ) K=K+1 CONTINUE CALL IMAT5(21,II,JJ,VAL) IF(IERR.EQ.O) GOTO 60 WRITE(6,LIST) IERR,' IMAT5',NE STOP CONTINUE 60 C C Adjust S t i f f n e s s matrix f o r p r e s c r i b e d values and create RHS. C DO 100 1=1,NOD ZVL=ZND(I,3) IF(ZVL.EQ.5.D0.OR.ZVL.EQ.2.D0) THEN IF(ZVL.EQ.5.D0) PF=O.DO IF(ZVL.EQ.2.D0) THEN IF (NQ.EQ.l) THEN IF(LFG.EQ.O) RHS(I)=FLUX*RHS(I) GOTO 100 ELSE PF=PSI GOTO 103 END IF END I F END IF IF(ZVL.EQ.4.D0) THEN IF(S0W.EQ.-99.D0)G0T0 101 PF=SOW GOTO 103 END IF IF(ZVL.EQ.8.D0) THEN PF=SOW IF(SOW.Eq.-99.DO)PF=0.DO GOTO 103 END IF IF(ZVL.EQ.6.AND.S0W.Eq.-99.D0) THEN PF=0.D0 GOTO 103 END IF GOTO 101 103 RHS(I)=BIG*PF CALL IAIJ5(I,I,BIG) GOTO 100 101 RHS(I)=0.D0 100 CONTINUE C C Insert RHS C CALL IRHS(RHS) C C Solve the s e t of equations. C CALL S0LV5(XAN) IF(IERR.EQ.O) GOTO 190 WRITE(6,LIST) IERR,' S0LV5 , LFG = ' ,LFG STOP 190 CONTINUE 1020 F0RMAT(1X,T10,I4,T20,F9.3,T35,F9.3,T50,E12.4) 1000 F0RMAT(T10,I3,T20,E10.3) IF(LFG.EQ.O) GOTO 220 EMX=0.D0 DO 210 1=1,NOD IF(DABS(OLD(D).LT.l.D-10) GOTO 210 154 Appendix C. FEM Software 210 AB=DABS(XAN(I)-OLD(I)) PCG=AB/DABS(0LD(I))*100.D0 IF(PCG.LT.EMX) GOTO 210 IMK=I EMX=PCG CONTINUE IF(EMX.LE.ERR) GOTO 2000 IF(EMX.GT.(1.D4*EMX0)) GOTO 3000 EMX0=EKX 220 LFG=1 NIT=NIT+1 IF(NIT.GT.MAXIT) GOTO 3000 DO 230 1=1,NOD SOR( I) =OLD ("I) +SKL* (XAN (I) -OLD (I) ) OLD(I)=XAN(I) 230 CONTINUE GOTO 45 2000 CONTINUE WRITE(6,1050) NIT.ERR.EMX.SKL 1050 F0RMAT(T10,' *** SUCCESSFUL CONVERGENCE IN ',12,' ITERATIONS' fi,': ERROR LIMIT ' ,F8.4, '*/.',/, TlO,' EMX = '.F8.4, fi> RELAXATION FACTOR = \F5.3) C C P r i n t out c h a r a c t e r i s t i c v a r i a b l e s . C 1060 F0RMAT(T30,'*** POWI.CYL ***',//,T20,'PARTICLE VARIABLES' *,/,T5,'RADIUS',T15,'SURFACE POTENTIAL',T35,'CONTACT ANGLE',/) 1061 F0RMAT(T30,'*** POWI.CYL ***',//,T20,'PARTICLE VARIABLES' *,/,T5,'RADIUS',T15,'SURFACE CHARGE',T35,'CONTACT ANGLE', *T50,'PSI i n f i n i t y ' , / ) 1065 FORMAT(T5.F8.3,T15,F8.3,T35,F8.3,T50,F8.3//) 1070 FORMAT(T20,'GRID DIMENSIONS',/,T7,'XD2',T17,'OIL DEPTH', *T35,'WATER DEPTH',T50,»X-INFIN',T60,/) 1075 F0RMAT(T5,F8.3,T15,F8.3,T35,F8.3,T50,F8.3,/) 1077 FORMAT(T5,'FZWA',T15,'FZWP',T25,'FZOP',T35,'FZOB',T45,'FCZ',T55 <D,'FCXI',T65,'FCXD') 1078 F0RMAT(T2,F8.5,T12,F8.5,T22,F8.5,T32,F8.5,T42,F8.5,T52,F8.5, fi T62.F8.5,/) 1080 FORMAT(T20,'GRID VARIABLES',/,T5,'NEA',T15,'NEI',T25,'NED', *T35,'NWA',T45,'NOB',T55,'NODES',T65,'NEL'/) 1085 F0RMAT(T2,I5,T12,I5,T22,I5,T32,I5,T42,I5,T52,I7,T62,I5//) IF(NQ.EQ.O) THEN WRITE(6,1060) WRITE(6,1065) RA/SR,PSI,(180.D0-ALP*180.DO/PI) ELSE WRITE(6,1061) WRITE(6,1065) RA/SR,QS,(180.DO-ALP*180.DO/PI),PSI END IF WRITE(6,1070) WRITE(6,1075) XD2/SR,Z0/SR,ZW/SR,XI/SR WRITE(6,1077) WRITE (6,1078 ) FZWA, FZWP, FZOP, FZOB, FCZ, FCXI, FCXD WRITE(6,1080) WRITE(6,1085) NEA,NEI,NED,NWA,NOB,NOD,NEL C C Load r e s u l t s i n t o an array i n a common b l o c k . C DO 245 1=1,NOD SOR(I)=XAN(I) 245 CONTINUE C 155 156 Appendix C. FEM Software C C C C C C Gather elements along boundaries t o be i n t e g r a t e d **** HORIZONTAL FORCE FIRST **** over. I f the 0/W i n t e r f a c e i s charged then the h o r i z o n t a l f o r c e at i n - f i n i t y must be subtracted from the midplane r e s u l t . SIN1=0.D0 SIN2=0.D0 SIN3=0.D0 IF(S0W.EQ.-99.D0) GOTO 250 C C Gather elements along boundaries t o be i n t e g r a t e d C C Plane at i n f i n i t y . C Section 1 NST=NWAI+3*NEA NUM=NWA LF=3 NE=NST DO 251 1=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 251 CONTINUE C Section 3 NST=NST+2*NEI*NE0 NUM=NE0 LF=3 NE=NST DO 252 J=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 252 CONTINUE C Section 2 NST=NST+2*NED*(NEW+NEO)+NWAD+N0BD+2*NEW*NEI NUM=NEW LF=1 NE=NST DO 253 J=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 253 CONTINUE C Section 4 NST=NST+NOBI PRINT *,NST NUM=N0B LF=.l NE=NST DO 254 J=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 254 CONTINUE C KSURF=1 IBE(I,1)=-1 NGF=3 CALL BNDINT(IBE,KSURF,SIN1,SIN2,SIN3,NGF) over. Appendix C. FEM Software 250 CONTINUE C C Median Plane. C Section 5 LF=3 NST=3*NEA+NWAI+2*NE0*NEI+NWAD NUM=NWA NE=NST DO 260 1=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 260 CONTINUE C Section 6 NST=NST+2*NED*(NEO+NEW) NUM=NEW LF=2 NE=NST DO 265 J=1,NUH IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 265 CONTINUE C Section 7 NST=3*NEA+NWAI+2*NE0*(NED+NEI)+NWAD NUM=NE0 NE=NST LF=3 DO 270 J=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 270 CONTINUE C Section 8 NE=NST+2*NEW*NED+N0BD NUM=N0B LF=2 DO 275 J=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE-1 1=1+1 275 CONTINUE NGF=3 KSURF=2 IBE(I,1)=-1 CALL BNDINT(IBE,KSURF,SUM1,SUM2,SUM3,NGF) RSTR=XD2*2.D0/SR FDIM=(CP+ZM)/(PSI**2*EM(2))*SR FMD1=SUM1*FDIM FMD2=SUM2*FDIM FMD3=SUM3*FDIM FMD4=FMD2+FMD3FIN1=SIN1*FDIM FIN2=SIN2*FDIM FIN3=SIN3*FDIM FIN4=FIN2+FIN3 PC0IL=FMD3/FMD4*100 .DO IF(S0W.EQ.-99.D0) THEN WRITE(6,1100) NGF,RSTR,FMD1,FMD2,FMD3,FMD4,PCOIL Appendix C. FEM Software 1100 F0RMAT(/,T10,» ** HORIZONTAL FORCE COMPONENT FROM ', « ' MEDIAN PLANE UNCHARGED INTERFACE** ',//, Q T15,' Number of i n t e g r a t i o n p o i n t s / element : ',12,/, Q T15,' Separation of surfaces : ',F8.4,/, Q T15,' Force - med plane ( l i n e avg) TOTAL : '.E10.4,/, 0 T15,' Force - med plane (element avg) H20 phase: '.E10.4,/, C TIB,' Force - med plane (element avg) o i l phase: ',E10.4,/, 8 T15,' Force - med plane (element avg) TOTAL : '.E10.4,/, Q T15,' '/, c o n t r i b u t i o n from o i l phase t o TOTAL : '.E10.4,/) ELSE FIN1=SIN1*FDIM FIN2=SIN2*FDIM FIN3=SIN3*FDIM FIN4=FIN2+FIN3 FRST=FMD4-FIN4 FRSW=FMD2-FIN2 FRS0=FMD3-FIN3 WRITE(6,1125) S0W,NGF,RSTR,FIN1,FIN2,FIN3,FIN4,FMD1,FMD2, #FMD3,FMD4,FRSW,FRSO,FRST 1125 F0RMAT(/,T10,' ** HORIZONTAL FORCE COMPONENT FROM', «' MEDIAN PLANE : 0/W SURFACE POTENTIAL = '.F8.3,' *** ',//, <D T15,' Number of i n t e g r a t i o n p o i n t s / element : ',12,/, 0 T15,' Separation of surfaces : ',F8.4,/, C T10,' LINE',T25,' H20',T40,' 0IL',T55,' TOTAL',// <0 Tl,'INFIN',T10,E11.5,T25,E11.4,T40,E11.4,T55,E11.4,/ <B Tl,'MEDIAN',T10,E11.5,T25,E11.4,T40,E11.4,T55,E11.4,// C T5,'NET FORCES:',T25,E11.4,T40,E11.4,T55,E11.4,/) END IF C C **** VERTICAL FORCE **** C C The c o n t r i b u t i o n from the separation side and the backside C of the p a r t i c l e are c a l c u l a t e d separately f o r comparison purposes. C As w e l l , the same q u a n t i t i e s are also c a l c u l a t e d from both the o i l C and water s i d e s : these should be equalwhen the 0/W i n t e r f a c e i s C uncharged. C C 0/W i n t e r f a c e , waterside f i r s t . C Front side of c y l i n d e r . IF(S0W.NE.-99.D0) GOTO 3400 C Section 6 LF=3 NST=3*NEA+NWAI+2*NE0*(NED+NEI)+NWAD+2*NEW NUM=NED NE=NST DO 280 I=1,NUM IBE(I,1)=NE IBE(I,2)=LF .NE=NE+2*NEW 280 CONTINUE IBE(I,1)=-1 NGF=3 KSURF=3 CALL BNDINT(IBE,KSURF,SWF1,SWF2,SZFW,NGF) C Backside C Section 2 NST=NST+2*NEW*NED+N0BD NUM=NEI LF=3 NE=NST DO 285 1=1,NUM IBE(I,1)=NE IBE(I,2)=LF Appendix C. FEM Software 285 NE=NE+2*NEW CONTINUE KSURF=3 IBE(I,1)=-1 CALL BNDINT(IBE,KSURF,SWB1,SWB2,SZBW,NGF) C C O i l s i d e : should be the same as waterside r e s u l t . C Frontside C Section 3 NST=3*NEA+NWAI+NE0+1 NUM=NEI NE=NST LF=1 DO 290 1=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE+2*NE0 290 CONTINUE KSURF=4 IBE(I,1)=-1 CALL BNDINT(IBE,KSURF,S0B1,S0B2,SZBO,NGF) C Section 7 NST=NST+2*NE0*NEI+NWAD NE=NST NUM=NED LF=2 DO 295 1=1,NUM IBE(I,1)=NE IBE(I,2)=LF NE=NE+2*NE0 295 CONTINUE KSURF=4 IBE(I,1)=-1 CALL BNDINT(IBE,KSURF,S0F1,S0F2,SZFO,NGF) FDIM=(CP+ZM)/EM(2)/PSI**2*SR FWFT=SWF1*FDIM FWBT=SWB1*FDIM FWF0=SWF2*FDIM FWB0=SWB2*FDIM FWTOT=FWFT+FWBT F0FT=S0F1*FDIM F0BT=S0B1*FDIM F0F0=S0F2*FDIM F0B0=S0B2*FDIM F0T0T=F0FT+F0BT IF (SOW.EQ.-99.DO) THEN PCDIF=(F0T0T-FWT0T)/F0T0T*100. RAT1=SZF0/SZFW RAT2=SZB0/SZBW WRITE(6,1400) FWFT,FWBT,FWFO,FWBO,FWTOT,FOFT,FOBT,FOFO, */. F0B0,F0T0T,PCDIF,RAT1,RAT2 1400 F0RMAT(/,T10,' ** VERTICAL FORCE COMPONENTS FROM', <D' 0/W PLANE ** ',//, CT5,' WATER SIDE RESULTS ',/, 0T2,'FR TOTAL',T15,'BK TOTAL',T30,'FR OSMOT.',T45,'BK OSMOT <D T60,' TOTAL',/, CT2,E10.4,T15,E10.4,T30,E10.4,T45,E10.4,T60,E10.4,//, CT5,' OIL SIDE RESULTS ',/, 0T2,'FR TOTAL',T15,'BK TOTAL',T30, 'FR OSMOT.',T45,'BK OSMOT 6 T60,' TOTAL',/, <DT2,E10.4,T15,E10.4,T30,E10.4,T45,E10.4,T60,E10.4,//, OT10,' **** '/.DIFFERENCE BETWEEN TOTALS ****' .T60.E10.4,// <BT15,' RATIO OF INTEGRATED NORMAL DERIVATIVES OIL/WATER',// Appendix C. FEM Software flT5,' FRONT SIDE = '.E10.4,' BACKSIDE = '.E10.4) ELSE TOTAL=FWTOT+FOTOT WRITE(6,1410) FWFT, FWBT,FWTOT,FOFT,FOBT,FOTOT,TOTAL 1410 F0RMAT(/,T10,' ** VERTICAL FORCE COMPONENTS FROM', fl' 0/W PLANE FOR CONSTANT CHARGE CASE ** ',//, flT5,' WATER SIDE RESULTS ',/, «1T2,'FR TOTAL',T15,'BK TOTAL',T30,'TOTAL',/, flT2,E10.4,T15,E10.4,T30,E10.4,//, flT5,' OIL SIDE RESULTS ',/, CT2,'FR TOTAL',TIS.'BK TOTAL',T30,'TOTAL',/, flT2,E10.4,T15,E10.4,T30,E10.4,//, flTIO,' **** TOTAL VERTICAL FORCE:'.T40.E10.4,//) END IF C P r i n t out s e l e c t values of PSI (when next l i n e i s commented o u t ) . IF(NOD.GT.l) GOTO 3400 NC1=0 NC2=0 EL=l.D-6 DO 300 1=1,NOD C IF(DABS(ZND(I,1)-(0RX-RA-XD2)).LE.EL)THEN IF(ZND(I,1) .LE.(-RA-ORX+EL).AND.DABS(ZND(I,2)-0RZ) $.LE.EL)THEN NC1=NC1+1 ND1(NC1)=I END IF C IF(DABS(ZND(I,1)-(0RX+XI)).LE.EL) THEN IF(ZND(1,2) .GE. (RA+ORZ-EL) . AND.DABS(ZND(1,1)-0RX) $.LE.EL)THEN NC2=NC2+1 ND2(NC2)=I END IF 300 CONTINUE WRITE(6,1105) NCI DO 310 1=1,NCI QR=DABS(ZND(ND1(I),1)/SR) WRITE(6,1110) QR,XAN(ND1(I)) 310 CONTINUE WRITE(6,1105) NC2 DO 320 1=1,NC2 QR=DABS(ZND(ND2(I),2)/SR) WRITE(6,1110) QR,XAN(ND2(D) 320 CONTINUE 1105 F0RMAT(1X,I3) 1110 F0RMAT(2F14.10) 3400 CONTINUE C C Write s o l u t i o n data f o r contour p l o t s (when next l i n e i s commented o u t ) . IF(NOD.GT.l) GOTO 3500 WRITE(8,1300) NEL,NOD,NEA,NOBD,NOBI,NWAD,NWAI,NEW,NEO,NED,NEI WRITE(8,1305) PSI,ZW/SR,XI/SR,XD2/SR,RA/SR,ZP/SR WRITE(8,1310) ((KEY(NE,J),J=1,6),NE=1,NEL) WRITE(8,1320) (ZND(J,1)/SR,ZND(J,2)/SR,XAN(J),J=1,N0D) 1300 F0RMAT(1X,11I6) 1305 F0RMAT(1X,6(1X,F10.4)) 1310 F0RMAT(1X,6(1X,I6)) 1320 F0RMAT(1X,3(1X,E15.5)) GOTO 3500 3000 CONTINUE WRITE(6,1055) NIT,ERR,EMX,SKL 1055 F0RMAT(T10,' ??? UNSUCCESSFUL CONVERGENCE IN ',12,' ITERATIONS' Appendix C. FEM Software Q ,':ERR0R LIMIT = ' ,F8.4,"/.' ,/,T10,' « ' Relaxation Factor = '.F8.4) 3500 CONTINUE END EMX = '.F8.4.T30, C SUBROUTINE SETMAT(LFG) C C C C C C C C C C Subroutine t o f i l l the lower t r i a n g l e of the elemental matrices. Upon f i r s t entry a l l the submatrices a r e c a l c u l a t e d . During subsequent c a l l s only the element submatrices corresponding t o the water phase domain are r e c a l c u l a t e d . INPUT : LFG - f l a g f o r non-linear convergence. LFG=0 ; f i r s t pass of f i t , use l i n e a r approx. LFG=1 ; subsequent passes, use p r e v i o u s nodal p o t e n t i a l s i n s e r i e s expansion of s i n k term. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Variables: C GWT C L1.L2.L3 C NG C B i n t e g r a t i o n wieghts. integration points. i n t e g r a t i o n order. b a s i s f u n c t i o n values at i n t e g r a t i o n points. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c C C C C C IMPLICIT REAL*8 (A-H.O-Z) CHARACTER* 1 LIST REAL*8 L1(4),L2(4),L3(4) DIMENSION GWT(4),B(6) DATA NG,LIST/4,'*'/ DATA L1/0.5D0,0.5D0,0.D0,0.D0/ DATA L2/O.DO,0.5DO,0.5DO,O.DO/ DATA L3/0.5D0,O.DO,0.5D0,O.D0/ DATA GWT/0.3333333333333333D0.0.3333333333333333D0, 00.3333333333333333D0,0.3333333333333333D0/ DATA L1/0.3333333333333333D0.0.6DO,0.2DO,0.2DO/ DATA L2/0.3333333333333333D0,0.2D0,0.6D0,0.2D0/ DATA L3/0.3333333333333333D0,0.2D0,0.2D0,0.6D0/ C COMMON /BLK1/ZND(5000,3),SBM(2500,21),SOR(5000),SEXP(20), (DSqL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD, ffl.NOBI.NOBD COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2) ,AJIN(2,2) i,DETJ,EM(3) C GWT(1)=-27.DO/48.DO GWT(2)=25.D0/48.D0 GWT(3)=GWT(2) GWT(4)=GWT(2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NUMERICAL INTEGRATION SECTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c C Perform quadrature over t r i a n g u l a r elements. C C C EVALUATE CONTRIBUTION OF EACH POINT. C DO 100 11=1,NG B(1)=L1(II)*(2.D0*L1(II)-1.D0) B(2)=L2(II)*(2.D0*L2(II)-1.D0) Appendix C. FEM Software B(3)=L3(II)*(2.D0*L3(II)-1.D0) B(4)=4.D0*L1(II)*L2(II) B(5)=4.D0*L3(II)*L2(II) B(6)=4.D0*L1(II)*L3(II) DBS(1)=4.D0*L1(II)-1.D0 DBS(2)=0.D0 DBS(3)=1.D0-4.D0*L3(II) DBS(4)=4.D0*L2(II) DBS(5)=-4.D0*L2(II) DBS(6)=4.D0*(L3(II)-L1(II)) DBT(1)=0.D0 DBT(2)=4.D0*L2(II)-1.D0 DBT(3)=1.D0-4.D0*L3(II) DBT(4)=4.D0*L1(II) DBT(5)=4.D0*(L3(II)-L2(II)) DBT(6)=-4.D0*L1(II) C C CALCULATE CONTRIBUTION OF EACH ELEMENT FOR THE QUADRATURE POINT. C DO 50 NE=1,NEL IF(NOTE(NE).NE.2.AND.LFG.EQ.1) GOTO 50 CALL JACOB(NE,II,LFG) IF(N0TE(NE).Eq.2) THEN IF(LFG.EQ.l) THEN YE=0.D0 DO 45 1=1,6 YE=YE+B(I)*SOR(KEY(NE,I)) 45 CONTINUE SERIES=SEXP(1) DO 46 1=2,20 IM=I-1 SERIES=SERIES+SEXP(I)*YE**IM 46 CONTINUE SERIES=SERIES/(2.DO+ZM) ELSE SERIES=SQL END IF END IF C C FILL MATRIX C K=l DO 40 1=1,6 DO 40 J=1,I IF(N0TE(NE).Eq.2) THEN FST=((DBX(I)*DBX(J)+DBZ(I)*DBZ(J))*EM(NOTE(NE))+SERIES # * B ( I ) * B ( J ) ) *GWT(II)*DETJ ELSE FST=(DBX(I)*DBX(J)+DBZ(I)*DBZ(J))*DETJ*GWT(II)*EM(NOTE(NE)) END IF SBM ( NE, K) =FST+SBM ( NE, K ) C WRITE(6,LIST) FST K=K+1 40 CONTINUE 50 CONTINUE 100 CONTINUE RETURN END C SUBROUTINE JACOB(NE,II,LFG) C C Program t o c a l c u l a t e the Jacobian plus x & z d e r i v a t i v e s of the C Basis f u n c t i o n s . Appendix C. FEM Software c C Input NE - Current element being considered. C I I - Quadrature p o i n t . C LFG - i t e r a t i o n f l a g (only check signs on f i r s t pass) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Variables : C SJAC array f o r checking a s i g n change of the determinate. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c IMPLICIT REAL*8 (A-H.O-Z) DIMENSION SJAC(2500,5) COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), BSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD, C.NOBI.NOBD COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2),AJIN(2,2) 1,DETJ,EM(3) C C C a l c u l a t e the Jacobian,AJAC, i t s determinate.DETJ, and i t s i n v e r s e , C AJIN. C DO 30 1=1,2 DO 30 J=l,2 AJAC(I,J)=0.D0 30 CONTINUE C DO 40 JK=1,6 AJAC(1,1)=AJAC(1,1)+DBS(JK)*ZND(KEY(NE,JK),1) AJAC(1,2)=AJAC(1,2)+DBS(JK)*ZND(KEY(NE,JK),2) AJAC(2,1)=AJAC(2,1)+DBT(JK)*ZND(KEY(NE,JK),1) AJAC(2,2)=AJAC(2,2)+DBT(JK)*ZND(KEY(NE,JK),2) 40 CONTINUE DETJ=AJAC(1,1)*AJAC(2,2)-AJAC(1,2)*AJAC(2,1) IF(LFG.NE.O) GOTO 45 C C Check f o r a s i g n change i n the determinate. C NS=0 IF(DETJ.GT.O.DO) NS=1 IF(II.GT.l) THEN IM=II-1 IF(NS.NE.SJAC(NE,IM)) THEN WRITE(6,1010) NE WRITE(6,1005) (ZND(KEY(NE,J),1),J=1,4) WRITE(6,1005) (ZND(KEY(NE,J),1),J=5,6) WRITE(6,1005) (ZND(KEY(NE,J),2),J=1,4) WRITE(6,1005) (ZND(KEY(NE,J),2),J=5,6) END IF END IF SJAC(NE,II)=NS 1010 FORMAT(IX,' SIGN CHANGE FOR ELEMENT ',13) IF(DETJ.NE.O.DO) GOTO 45 WRITE(6,1000) NE.DETJ WRITE(6,1005) AJAC(1,1),AJAC(1,2),AJAC(2,1),AJAC(2,2) WRITE(6,1005) (DBS(J),J=1,4) WRITE(6,1005) (DBS(J),J=5,6) WRITE(6,1005) (ZND(KEY(NE,J),1),J=l,4) WRITE(6,1005) (ZND(KEY(NE,J),1),J=5,6) WRITE(6,1005) (ZND(KEY(NE,J),2),J=1,4) WRITE(6,1005) (ZND(KEY(NE,J),2),J=5,6) 1000 FORMAT(T10,' ELEMENT #',14,' HAS DETERMINATE = '.E10.3) 1005 FORMAT(T5.4E13.4) Appendix: C. FEM Software c 45 AJIN(1,1)=AJAC(2,2)/DETJ AJIN(1,2)=-AJAC(1,2)/DETJ AJIN(2,1)=-AJAC(2,1)/DETJ AJIN(2,2)=AJAC(1,1)/DETJ C C C a l c u l a t e DB/DX.ftDB/DZ f o r each node i n the element. C DO 50 1=1,6 DBX(I)=DBS(I)*AJIN(1,1)+DBT(I)*AJIN(1,2) DBZ(I)=DBS(I)*AJIN(2,1)+DBT(I)*AJIN(2,2) 50 CONTINUE RETURN END C SUBROUTINE BNDINT(IBE,KS,SUM1,SUM2,SUM3,NG) C C Program -that performs the boundary i n t e g r a l s r e q u i r e d f o r the C the f o r c e c a l c u l a t i o n s . C C INPUTS: IBE(N.l) - l i s t of elements containing boundary. C IBE(N,2) - f l a g i n d i c a t i n g whether L2 or L3 C i s zero on boundary. C KS - surface (1=N/A , 2=median p l a n e ) . C (3 = waterside, 4 = o i l s i d e ) C OUTPUT: SUM(x) - value of the t o t a l i n t e g r a l and the C c o n t r i b u t i o n s t o the i n t e g r a l (see code). C NG - i n t e g r a t i o n order. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Variables: C GW i n t e g r a t i o n weights. C GP integration points. C IN key array f o r boundary i n t e g r a t i o n s . C B b a s i s f u n c t i o n values at i n t e g r a t i o n p t s . C DBL basis f u n c t i o n d e r i v a t i v e s . C PE p o t e n t i a l value at i n t e g r a t i o n p t . C XE x co-ordinate a t i n t e g r a t i o n p t . C ZE z co-ordinate a t i n t e g r a t i o n p t . C DPDX p o t e n t i a l d e r i v a t i v e w.r.t. x. C DPDZ p o t e n t i a l d e r i v a t i v e w.r.t. z. C DRDZ s e r e n d i p i t y co-ord. d e r i v . w.r.t. z. C DC i n f i n i t e s m a l a r c length along a boundary. C OP osmotic term at the i n t e g r a t i o n p o i n t . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IMPLICIT REAL*8(A-H,0-Z) DIMENSION IBE(100,2),GP(5),GW(5),B(3),IN(6),DBL(3) REAL*8 L1.L2.L3 COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), CSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD, 0. NOBI.NOBD COMMON /BLK2/Z0,ZW,XI,RA,RB,ORX,ORZ,ALP,ZP,PI,PSI,CP,SOW COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2),AJIN(2,2) 1, DETJ,EM(3) C GOTO (1,2,3,4,5),NG 1 GW(1)=2.D0 GP(1)=0.D0 GOTO 9 2 GW(1)=1.D0 GW(2)=1.D0 Appendix C. FEM Software GP(1 =-.577350269189626 GP(2 =0.577350269189626 GOTO 9 GW(1 =5.DO/9.DO GW(3 =GW(1) GW(2 =8.DO/9.DO GP(1 =-.774596669241483 GP(3 =0.774596669241483 GP(2 =0.D0 GOTO 9 G«(l =0.347854845137454 GW(4 =0.347854845137454 GW(2 =0.652145154862546 GW(3 =0.652145154862546 GP(1 =-.861136311594053 GP(4 =0.861136311594053 GP(2 =-.339981043584856 GP(3 =0.339981043584856 GOTO 9 GW(1 =0.236926885056189 GW(5 =0.236926885056189 GW(2 =0.478628670499366 GW(4 =0.478628670499366 GW(3 =0.568888888888889 GP(1 =-.906179845938664 GP(5 =0.906179845938664 GP(2 =-.538469310105683 GP(4 =0.538469310105683 GP(3 =0.D0 SUM1=0.D0 SUM2=0.D0 SUM3=0.D0 SUM4=0.DO SUM5=0.D0 SR=DSQRT(EM(2)) EDF=(EM(2)-EM(3))/2.D0 LFG=1 C Perform i n t e g r a t i o n . C DO 100 1=1,50 NE=IBE(I,1) EFC=EM(NOTE(NE)) IF(NE.LT.O) GOTO 110 DO 50 11=1,NG R=GP(II) IF(IBE(I,2).EQ.3) THEN L1=(1.D0-GP(II))/2.D0 L2=l.DO-LI L3=0.D0 DO 10 J=l,6 IN(J)=J 10 CONTINUE GOTO 18 END IF IF(IBE(I,2).EQ.2) THEN IN(1)=3 IN(2)=1 IN(3)=2 IN(4)=6 IN(5)=4 IN(6)=5 Appendix C. FEM Software 18 19 20 L1=(1.D0+GP(II))/2.D0 L2=0.D0 L3=1.D0-L1 GOTO 18 END I F IF(IBE(I,2).EQ.l) THEN IN(1)=2 IN(2)=3 IN(3)=1 IN(4)=5 IN(5)=6 IN(6)=4 L2=(1.D0-GP(II))/2.D0 L1=0.D0 L3=1.D0-L2 END IF B(1)=-.5D0*R*(1.D0-R) B(2)=.5D0*R*(1.D0+R) B(3)=(1.D0+R)*(1.D0-R) DBL(1)=R-.5D0 DBL(2)=R+.5D0 DBL(3)=-2.D0*R DBS((1))=4.D0*L1-1.D0 DBS((2))=0.D0 DBS((3))=1.D0-4.D0*L3 DBS((4))=4.D0*L2 DBS((5))=-4.D0*L2 DBS((6))=4.D0*(L3-L1) DBT((1))=0.D0 DBT((2))=(4.D0*L2-1.D0) DBT((3))=(1.D0-4.D0*L3) DBT((4))=(4.D0*L1) DBT((5))=(4.DO*(L3-L2)) DBT((6))=(-4.D0*L1) PE=O.DO XE=O.DO ZE=O.DO DPDXL=0.D0 DPDX=0.D0 DPDZ=O.DO DRDZ=0.D0 DX=0.D0 DZ=0.D0 JJ=0 IF(KS.EQ.2.0R.KS.EQ.1)THEN DO 19 J=l,4 IF(J.Eq.3)G0T0 19 JJ=JJ+1 DRDZ=DRDZ+DBL(JJ)*ZND(KEY(NE,IN(J)),2) CONTINUE DRDZ=1.D0/DRDZ END IF JJ=0 DO 20 J=l,4 IF(J.EQ.3)G0T0 20 JJ=JJ+1 PE=PE+B(JJ)*SOR(KEY(NE,IN(J))) XE=XE+B(JJ)*ZND(KEY(NE,IN(J)),1) ZE=ZE+B(JJ)*ZND(KEY(NE,IN(J)),2) DX=DX+DBL(J J)*ZND(KEY(NE,IN(J)),1) DZ=DZ+DBL(JJ)*ZND(KEY(NE,IN(J)),2) DPDXL=DPDXL+DBL(J J)*SOR(KEY(NE,IN(J))) CONTINUE Appendix C. FEM Software DPDXL=DPDXL*DRDZ CALL JACOB(NE,II,LFG) DO 25 J=l,6 DPDX=DPDX+DBX((J))*SOR(KEY(NE,(J))) DPDZ=DPDZ+DBZ((J))*SOR(KEY(NE,(J))) 25 CONTINUE DC=DSQRT(DX**2+DZ**2) OP=O.DO IF(KS.EQ.2.0R.KS.EQ.l) THEN IF(ZE.GT.ZP) THEN OP=(ZM*(DEXP(-CP*PE)-1.DO)+CP*(DEXP(ZM+PE)-1))/(ZM**2*2.DO '/. *CP) C SUM2=SUM2+(EFC*(DPDX**2+DPDZ**2)/2.DO+OP)*DZ*GW(II) SUM2=SUM2+0P*DZ*GW(II) C SUM4=SUM4+(EFC*DPDXL**2/2.DO+OP)*DZ*GW(II) SUM4=SUM4+(EFC*DPDXL**2/2.DO)*DZ*GW(II) ELSE SUM3=SUM3+(EFC*(DPDX**2+DPDZ**2)/2.DO+OP)*DZ*GW(II) SUM5=SUM5+(EFC*DPDXL**2/2.DO+OP)*DZ*GW(II) END IF SUM1=SUM1+(EFC*DPDXL**2/2.D0+0P)*DZ*GW(II) GOTO 50 END IF 0P=(ZM*(DEXP(-CP*PE)-1.DO)+CP*(DEXP(ZM*PE)-1))/(ZM**2*2.DO '/. *CP) IF(KS.EQ.3) THEN IF(S0W.NE.-99.D0) THEN SUM1=SUM1-(EM(2)/2-D0*DPDZ**2-0P)*DX*GW(II) SUM2=SUM2+(OP)*DX*GW(II) SUM3=SUM3-(EM(2)/2.D0*DPDZ**2)*DX*GW(II) ELSE SUM1=SUM1+(EDF*(DPDX**2+EM(2)/EM(3)*DPDZ**2)+0P)*DX*GW(II) SUM2=SUM2+(OP)*DX*GW(II) SUM3=SUM3+EDF*DPDX**2*DX*GW(II) SUM4=SUM4+EDF*DPDZ**2*DX*GW(II)*EM(2)/EM(3) SUM5=SUM5+DPDZ*DX*GW(II) END IF GOTO 50 END IF IF(KS.EQ.4) THEN IF(S0W.NE.-99.D0) THEN SUM1=SUM1+(EM(3)/2.D0*DPDZ**2)*DX*GW(II) ELSE SUM1=SUM1+(EDF*(DPDX**2+EM(3)/EM(2)*DPDZ**2)+0P)*DX*GW(II) SUM2=SUM2+(OP)*DX*GW(II) SUM3=SUM3+EDF*DPDX* * 2*DX*GW(11) SUM4=SUM4+EDF*DPDZ**2*DX*GW(II)*EM(3)/EM(2) SUM5=SUM5+DPDZ*DX*GW(II) END IF END IF 50 CONTINUE 100 CONTINUE 110 CONTINUE IF(KS.EQ.2) THEN PC2=SUM2/SUM1*100 PC4=SUM4/SUM1*100 WRITE(6,1046) SUM1,SUM2,SUM4,PC2,PC4 1046 FORMAT(/,T5,'TOTAL '.E10.3.T25,' 0S= '.E10.3.T43, '/.' DPDZ= '.E10.3,/, '/.T5,' '/. CONTRIBUTIONS: ' ,T30,E10.3,T50,E10.3,//) END IF IF(KS.GT.2.AND.S0W.EQ.-99.D0) THEN PC2=SUM2/SUM1*100 21 Appendix C. FEM Software PC3=SUM3/SUM1*100 PC4=SUM4/SUM1*100 WRITE(6,1044) SUM1,SUM2,SUM3,SUM4,PC2,PC3,PC4 1044 FORMAT(/,T5,'T0TAL= '.E10.3.T25,' 0S= '.E10.3.T43,' DPDX= '/•E10.3.T62,' DPDZ= '.E10.3,/, '/.TS,' '/. CONTRIBUTIONS: ',T30,E10.3,T50,E10.3,T69,E10.3,//) END IF SUM3=SUM5 RETURN END C SUBROUTINE LOAD(IBE) C C Program that i n i t i a l i z e s the RBS load v e c t o r with the r e s u l t o f C the f l u x boundary i n t e g r a l s . C C INPUT IBE - l i s t of elements containing boundary. C IMPLICIT REAL*8 (A-H.O-Z) DIMENSION IBE(100,2) COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), OSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP", NWA,NOB,NEL,NOD,NWAI,NWAD (D.NOBI.NOBD COMMON /BLK2/Z0,ZW,XI,RA,RB,ORX,ORZ,ALP,ZP,PI,PSI,CP,SOW COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2),AJIN(2,2) 1,DETJ,EM(3) COMMON /BLK6/XD2,FZWA,FZWP,FZOP,FZOB,FCZ,FCXI,FCXD COMMON /FVEC/RHS(5000) C DO 20 1=1,100 NE=IBE(I,1) IF (NE.LT.O) GOTO 30 N1=KEY(NE,1) N2=KEY(NE,2) N4=KEY(NE,4) RHS(N1)=1.D0 RHS(N2)=1.D0 RHS(N4)=4.D0 20 CONTINUE 30 CONTINUE RETURN C SUBROUTINE GRIDCY C C VERSION 2: MICHAEL P. LYNE MAR/89 C C Program t o generate a g r i d over the c r o s s - s e c t i o n o f a C c y l i n d e r on a f l a t i n t e r f a c e . One side of the g r i d corresponds C t o a symmetry plane. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINES CALLED CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SPACE spaces nodes a t increments on a s t r a i g h t l i n e . C FILSEC f i l l s each s e c t i o n ' s nodal co-ordinates and numbers. C FILKEY f i l l s KEY array f o r the g r i d . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C VARIABLE LIST CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DT angular increments of p i e s e c t i o n s i n p a r t i c l e . C NDOB # of nodes i n z - d i r ' n below p a r t i c l e i n o i l phase. Appendix C. FEM Software C NDOP # of nodes i n z - d i r ' n on p a r t i c l e p r o f i l e i n o i l phase. C NDWA # of nodes i n z - d i r ' n on p a r t i c l e p r o f i l e i n H20 phase. C NDWP # of nodes i n z - d i r ' n above p a r t i c l e i n H2D phase. C NDT # of nodes i n z - d i r ' n on p a r t i c l e p r o f i l e . C SEC contains nodal co-ords. and numbers f o r each s e c t i o n . C SIX,SIZ,S2X,S2Z contain end p t s . of l i n e s t o be i n t e r p o l a t e d . C YSPC contains p o i n t s f o r a l i n e . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT REAL*8 (A-H.O-Z) CHARACTER*1 LIST DIMENSION YSPC(51) DATA LIST/'*'/ COMMON /BLK5/SEC(2601,8,3),S1X(51),S1Z(51),S2X(51),S2Z(S1),IDZ(8) COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), CSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD <B,N0BI,N0BD COMMON /BLK2/Z0,ZW,XI,RA,RB,ORX,ORZ,ALP,ZP,PI,PSI,CP COMMON /BLK6/XD2,FZWA,FZWP,FZOP,FZOB,FCZ,FCXI,FCXD C C Preliminary c a l c u l a t i o n s . C NEA2=NEA/2 NDT=NEA+1 NDWP=2*NEW+1 ND0P=2*NE0 NDWA=2*NWA ND0B=2*N0B C C This code ensures that the node spacing i n z - d i r ' n i s equal f o r nodes e i t h e r C side of the l i n e d e l i m i t i n g the middle s e c t i o n s . C DZWP=FZWP*(RA-ZP)/(FZWP+l.D0+(FZWP+l.D0)/2.D0*(NDWP-3)) XFC=1.D0+(NDWA-2.DO)/2.DO FZWA=((ZW-RA)/DZWP-XFC)/XFC DZ0P=FZ0P*(RA+ZP)/(FZ0P+l.D0+(FZ0P+l.D0)/2.D0*(ND0P-2)) XFC=1.D0+(ND0B-2.D0)/2.D0 FZOB=((ZO-RA)/DZOP-XFC)/XFC Fl=(FZWP+l.D0+(FZWP+l.D0)/2.D0*(NDWP-3)) F2=FZWP*(FZWA+1.D0+(FZWA+1.DO)/2.DO*(NDWA-2)) F3=-(RA-ZP)/(ZW-RA) FDW=F1/F2+F3 IF(FDW.LT.O.DO) FDW=O.DO FIW=FDW Fl=(FZ0P+l.D0+(FZ0P+l.D0)/2.D0*(ND0P-2)) F2=FZ0P*(FZOB+1.D0+(FZ0B+i.D0)/2.D0*(NDOB-2)) F3=-(RA+ZP)/(ZO-RA) FD0=F1/F2+F3 IF(FDO.LT.O.DO) FD0=0.D0 FI0=FD0 C C C a l c u l a t e the number of elements i n the domain. C NWAD=2*NED*NWA NWAI=2*NEI*NWA N0BD=2*NED*N0B N0BI=2*NEI*N0B C NEL=3*NEA+NEA*(NEI+NED)+NWAI+NWAD+N0BI+N0BD C C F i l l phase i n d i c a t i n g array: NOTE. 169 Appendix C. FEM Software 20 C DO 20 NE=1,3*NEA N0TE(NE)=1 CONTINUE NE=NE-1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C CALCULATE NODAL CO-ORDINATES C C C a l c u l a t e the nodal co-ordinates of each of the C 8 sections surrounding the sphere. Each of the sections c o n s i s t s C of groups of nodes which l i e on l i n e s with end points PIftP2. C To f i l l the 8 sections one must supply FILSEC with the arrays C of end p o i n t s (S1X,S1Z,S2X,S2Z) d e s c r i b i n g the l i n e s i n that ~C s e c t i o n . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Section 1 C YH=ZW+0RZ YL=RA+0RZ NDL=NDHA+1 CALL SPACE(YH,YL,NDL,YSPC,FZWA) DO 50 1=1,NDL S1Z(I)=YSPC(I) S1X(I)=0RX 50 CONTINUE YL=FIW*(ZW-RA)+YL YH=YL+(ZH-RA) CALL SPACE(YH,YL,NDL,YSPC,FZWA) DO 55 1=1,NDL S2Z(I)=YSPC(I) S2X(I)=0RX+XI 55 CONTINUE NDR=2*NEI+1 CALL FILSEC(NDL,NDR.l) C C SECTION 2 C NDL=NDWP DT=PI/NEA DO 60 1=1,NDL TT=(I-1)*DT S1Z(I)=RA*DC0S(TT)+0RZ S1X(I)=RA*DSIN(TT)+0RX 60 CONTINUE YH=YL YL=ZP+0RZ CALL SPACE(YE,YL,NDL,YSPC,FZWP) DO 65 1=1,NDL S2Z(I)=YSPC(I) S2X(I)=0RX+XI 65 CONTINUE CALL FILSEC(NDL,NDR,2) C C S e c t i o n 3. C J=l DO 70 I=NDWP,NDT TT=(I-1)*DT S1Z(J)=RA*DC0S(TT)+0RZ S1X(J)=RA*DSIN(TT)+0RX J=J+1 Appendix C. FEM Software 70 75 CONTINUE YH=YL YL=-RA-FI0*(Z0-RA)+0RZ NDL=ND0P+1 CALL SPACE(YL,YH,NDL,YSPC,FZOP) DO 75 I=NDL,1,-1 S2Z(I)=YSPC(NDL-I+1) S2X(I)=0RX+XI CONTINUE CALL FILSEC(NDL,NDR,3) C C Section 4. C YH=YL YL=YH-(ZO-RA) NDL=ND0B+1 CALL SPACE(YL,YH,NDL,YSPC,FZOB) DO 80 I=NDL,1,-1 S2Z(I)=YSPC(NDL-I+1) S2X(I)=0RX+XI 80 CONTINUE YH=-RA+ORZ YL=ORZ-ZO CALL SPACE(YL,YH,NDL,YSPC,FZOB) DD 85 I=NDL,1,-1 S1Z(I)=YSPC(NDL-I+1) S1X(I)=0RX 85 CONTINUE CALL FILSEC(NDL,NDR,4) C C SECTION 5 C NDL=NDWA+1 NDR=2*NED+1 DO 90 I=1,NDL SiX(I)=SEC(I,l,l) S1Z(I)=SEC(I,1,2) 90 CONTINUE YL=FDW*(ZW-RA)+RA+ORZ YH=YL+(ZW-RA) CALL SPACE(YH,YL,NDL,YSPC,FZWA) DO 95 1=1,NDL S2Z(I)=YSPC(I) S2X(I)=0RX-RA-XD2 95 CONTINUE CALL FILSEC(NDL,NDR,5) C C Section 6 C NDL=NDWP DO 100 1=1,NDL S1Z(I)=SEC(I,2,2) SIX(I)=2.DO*ORX-SEC(I,2,1) 100 CONTINUE YH=YL - YL=ZP+ORZ CALL SPACE(YH,YL,NDL,YSPC,FZWP) DO 105 1=1,NDL S2Z(I)=YSPC(I) S2X(I)=0RX-RA-XD2 105 CONTINUE CALL FILSEC(NDL,NDR,6) Appendix C. FEM Software C Section 7 C NDL=ND0F+1 DO 110 1=1,NDL S1Z(I)=SEC(I,3,2) S1X(I)=2.D0*0RX-SEC(I,3,1) 110 CONTINUE YH=YL YL=0RZ-RA-FD0*(ZO-RA) CALL SPACE(YL,YH,NDL,YSPC,FZOP) DO 115 I=NDL,1,-1 S2Z(I)=YSPC(NDL-I+1) S2X(I)=0RX-RA-XD2 115 CONTINUE CALL FILSEC(NDL,NDR,7) C C Section 8 C NDL=ND0B+1 DO 120 1=1,NDL S1Z(I)=SEC(I,4,2) S1X(I)=SEC(I,4,1) 120 CONTINUE YH=YL YL=YH-(ZO-RA) CALL SPACE(YL,YH,NDL,YSPC,FZOB) DO 125 I=NDL,1,-1 S2Z(I)=YSPC(NDL-I+1) S2X(I)=0RX-RA-XD2 125 CONTINUE CALL FILSEC(NDL,NDR,8) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c C Assign co-ordinates t o nodes C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ZND(1,1)=0RX ZND(1,2)=0RZ ZND(1,3)=FLG NC=2 INC=3 DR=RB/2.D0 DT=PI/NEA NDA=2*NEA L0IL=NDA-NDWP+2 C C Triangular p i e region. C DO 130 N=l,2 RD=N*DR INC=INC-1 DO 130 I=1,NDA,INC TT=(I-1)*DT ZND(NC,3)=FLG ZND(NC,1)=RD*DSIN(TT)+ORX ZND(NC,2)=RD*DC0S(TT)+0RZ NC=NC+1 130 CONTINUE C C C i r c u l a r region. C DR=(RA-RB)/2.D0 DO 135 N=l,2 Appendix C. FEM Software RD=N*DR+RB DO 135 1=1,NDA TT=(I-1)*DT IF(N.EQ.2.AND.(I.LE.NDWP.0R.I.GE.L0IL))THEH ZND(NC,3)=2.D0 ELSE ZND(NC,3)=FLG END IF ZND(NC,1)=RD*DSIN(TT)+0RX ZND(NC,2)=RD*DC0S(TT)+ORZ NC=NC+1 CONTINUE 135 C C 'Z-axis' nodal co-ordinates. C ZND(NC,3)=5.D0 ZND(NC,1)=SEC(1,1,1) ZND(NC,2)=SEC(1,1,2) NC=NC+1 DO 140 I=2,NDWA ZND(NC,3)=FLG ZND(NC,1)=SEC(I,1,1) ZND(NC,2)=SEC(I,1,2) NC=NC+1 140 CONTINUE C DO 150 I=2,ND0B ZND(NC,3)=FLG ZND(NC,1)=SEC(I,4,1) ZND(NC,2)=SEC(I,4,2) NC=NC+1 150 CONTINUE ZND(NC,3)=5.D0 ZND(NC,1)=SEC(NDOB+1,4,1) ZND(NC,2)=SEC(ND0B+1,4,2) NC=NC+1 IDZ(l)=NDWA+l IDZ(5)=IDZ(1) IDZ(2)=NDVfP IDZ(6)=IDZ(2) IDZ(3)=ND0P+1 IDZ(7)=IDZ(3) IDZ(4)=ND0B+1 IDZ(8)=IDZ(4) C IF(NED.NE.NEI) THEN IF(NED.GT.NEI) THEN NLIM=2*NEI+1 NFN=2*NED MF=0 ELSE NLIM=2*NED+1 NFN=2*NEI MF=1 END IF ELSE NFN=2*NED NLIM=NFN+1 MF=2 END IF INC=0 IST=1 Appendix C. FEM Software IFN=8 C C Assign nodal co-ordinates f o r each of 8 s e c t i o n s . C DO 200 N=1,NFN INC=INC+1 IF(N.EQ.NLIM) THEN IF(MF.EQ.O) THEN IST=5 IFN=8 ELSE IST=1 IFN=4 END I F END I F DO 200 I=IST,IFN JST=N*IDZ(I)+2 JFN=(N+1)*IDZ(I) IF(I.Eq.l.0R.I.Eq.5)THEN JST=N*IDZ(I)+1 ELSE SEC(N*IDZ(I)+1,I,3)=DBLE(NC-1) END IF IF(JST.GT.JFN) GOTO 200 DO 200 J=JST,JFN ZND(NC,1)=SEC(J,I,1) ZND(NC,2)=SEC(J,I,2) ZND(NC,3)=FLG IF(I.Eq.2.0R.I.Eq.6) THEN IF(J.Eq.JFN)THEN ZND(NC,3)=4.D0 GOTO 190 END IF END IF IF(N.Eq.NLIM) THEN IF(MF.Eq.O) THEN IF(I.LT.5) ZND(NC,3)=6.DO ELSE IFQ.GT.4) ZND(NC,3)=1.D0 END IF GOTO 190 END IF IF(N.Eq.NFN)THEN IFQ.LT.5) ZND(NC,3)=6.D0 IF(I.GT.4) ZND(NC,3)=1.D0 END I F IF(I.EQ.1.0R.I.Eq.5)THEN IF(J.Eq.JST) THEN ZND(NC,3)=5.D0 GOTO 190 END IF END IF IF(I.Eq.4.0R.I.Eq.8)THEN IF(J.Eq.JFN) THEN ZND(NC,3)=8.D0 GOTO 190 END IF END IF 190 SEC(J,I,3)=DBLE(NC) NC=NC+1 200 CONTINUE N0D=NC-1 C 174 Appendix C. FEM Software C F i l l key code. C NE=1 C C Triangular portion f i r s t . C DO 250 11=1,HEA KEY(NE,1)=1 KEY(NE,3)=2+NEA+2*(N-l) KEY(NE,2)=KEY(NE,3)+2 IF(N.EQ.NEA) KEY(NE,2)=2+NEA KEY(NE,5)=KEY(NE,3)+1 KEY(NE,4)=N+2 IF(N.EQ.NEA) KEY(NE,4)=2 KEY(NE,6)=1+N NE=NE+1 250 CONTINUE C C Outer t r i a n g l e p o r t i o n . C KST=NEA+1 NFEL=3*NEA KNT=1 DO 260 NE=KST,NFEL-1,2 KEY(NE,1)=2+KNT*2+NEA KEY(NE,2)=KEY(NE,1)+4*NEA KEY (NE, 3) =KEY(NE, 1)-2 KEY(NE,6)=KEY(NE,1)-1 KEY(NE,4)=KEY(NE,1)+2*NEA KEY(NE,5)=KEY(NE,4)-1 IF(NE.EQ.NFEL-l) THEN KEY (NE, 1) =NEA+2 KEY(NE,4)=3*NEA+2 KEY(NE,2)=5*NEA+2 KEY(NE,5)=KEY(NE,2)-1 KEY(NE,6)=KEY(NE,4)-1 KEY(NE,3)=KEY(NE,6)-1 END IF KNT=KNT+1 260 CONTINUE KNT=1 DO 265 NE=KST+1,NFEL,2 KEY (NE, 3 ) =NEA+2*KNT KEY (NE, 2 ) =5*NEA+2*KNT KEY(NE,5)=3*NEA+2*KNT KEY(NE,4)=KEY(NE,2)+1 KEY(NE,1)=KEY(NE,4)+1 IF(NE.EQ.NFEL) KEY(NE,l)=5*NEA+2 KEY(NE,6)=KEY(NE,5)+1 KNT=KNT+1 265 CONTINUE C C F i l l key code f o r each of 8 sections; NOTE f i l l e d here too. C NE=NFEL+1 CALL FILKEY(NE,NDWA,NDOB) C C P r i n t out g r i d numbering scheme. C C500 FORMAT(//,TlO,'ELEMENT',T20,'ND1',T25,'ND2',T30,'ND3', C 1T35,'ND4',T40,'ND5',T45,'ND6',T50,//) 510 F0RMAT(T12,I3,T19,6I5) C WRITE(6,500) Appendix C. FEM Software C DO 300,I=1,NEL C WRITE(6,510) I,(KEY(I,J),J=1,6) C300 CONTINUE C WRITE(6,500) C DO 49,I=1,NEL C WRITE(6,520)I (ZND(KEY(I,J),1),J=1,6) C WRITE(6,520)I,(ZND(KEY(I,J),2),J=1,6) C WRITE(6,LIST) ' ' C49 CONTINUE 520 F0RMAT(T12,I3,T19,6(1X,E10.3)) C WRITE(6,LIST) 'NODES BEFORE EXITING GRID :',N0D RETURN END C SUBROUTINE FILSEC(NDZ,NDR,NSEC) C C Program t o c a l c u l a t e node p o i n t s along l i n e s described by the C end p o i n t s s u p p l i e d i n S1X,S1Z,S2X,S2Z. The nodal co-ordinates C are stored i n SEC. C C INPUTS: NDZ - # of nodes i n Z - d i r e c t i o n . C NDR ' along l i n e P1P2 C NSEC - s e c t i o n number being c a l c u l a t e d . C IMPLICIT REAL*8 (A-H.O-Z) DIMENSION RV(51) " COMMON /BLK5/SEC(2601,8,3),S1X(51),S1Z(51),S2X(51),S2Z(51),IDZ(8) COMMON /BLK6/XD2,FZWA,FZWP,FZOP,FZOB,FCZ,FCXI,FCXD C FCX=FCXD IF(NSEC.LT.5) FCX=FCXI DO 100 J=1,NDZ SLP=(S2Z(J)-S1Z(J))/(S2X(J)-S1X(J)) B=S1Z(J)-SLP*S1X(J) RL=DSQRT((S2Z(J)-S1Z(J))**2+(S2X(J)-S1X(J))**2) CALL SPACE(RL,0.D0,NDR,RV,FCX) NS=0 IF(S1X(J).GT.S2X(J))NS=1 DO 100 1=2,NDR-1 K=1+NDR-I IM=I-1 XAD=DSQRT(RV(K)**2/(1.D0+SLP**2)) IF(NS.EQ.O) THEN X=S1X(J)+XAD ELSE X=S1X(J)-XAD END IF Z=SLP*X+B NUM=J+IM*NDZ SEC(NUM,NSEC,1)=X SEC(NUM,NSEC,2)=Z 100 CONTINUE DO 110 J=1,NDZ SEC(J,NSEC,1)=S1X(J) SEC(J,NSEC,2)=S1Z(J) 110 CONTINUE JS=NDZ*(NDR-1)+1 JF=NDZ*NDR 1=1 DO 120 J=JS,JF SEC(J,NSEC,1)=S2X(I) SEC(J,NSEC,2)=S2Z(I) 1=1+1 > Appendix C. FEM Software 120 CONTINUE RETURN END C SUBROUTINE SPACE(YH.YL,NOD,YSPC,FC) C C Program t o space nodes i n decreasing order along a l i n e C bounded by YHftYL. The f i r s t spacing w i l l be FC times as great C as f i n a l spacing. C C INPUTS: YH - upper bound of l i n e . C YL - lower bound of l i n e . C NOD - # of nodes on the l i n e ( i n c l u d i n g ends) C (ihax. i s 50) C FC - spacing f a c t o r . C OUTPUT: YSPC- nodal co-ordinates. C IMPLICIT REAL*8 (A-H.O-Z) DIMENSION YSPC(51) C FAC=FC SEP=YH-YL SP=(FAC-1.DO)/(NOD-2) TNL=FAC+1.D0+(FAC+1.D0)/2.D0*(N0D-3) SCL=SEP/TNL YSPC(1)=YH YSPC(NOD)=YL DO 10 1=2,NOD-1 IM=I-1 YSPC(I)=YSPC(IM)-FAC+SCL FAC=FAC-SP 10 CONTINUE RETURN END C SUBROUTINE FILKEY(NE,NDWA,NDOB) C C Program t o f i l l the KEY array f o r the g r i d arrangement C created i n GRIDCY. C IMPLICIT REAL*8(A-H,0-Z) CHARACTER*1 LIST DATA LIST/'*'/ COMMON /BLK5/SEC(2601,8,3),S1X(51),S1Z(51),S2X(51),S2Z(51),IDZ(8) COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), CSQL,FLG,ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD (D.NOBI.NOBD COMMON /BLK2/Z0,ZW,XI,RA,RB,ORX,ORZ,ALP,ZP,PI,PSI,CP COMMON /BLK4/DBS(6),DBT(6),DBX(6),DBZ(6),AJAC(2,2),AJIN(2,2) 1,DETJ,EM(3) C C SECTIONS 1ft3 C MX=2*NEI DO 100 KK=1,2 IF (KK.EQ. 1) THEN MZ=NWA NS=1 KN=7*NEA+2 N0TE(NE)=2 Appendix C. FEM Software 10 15 ELSE MZ=NEO NS=3 KN=2+5*NEA+2*NEW N0TE(NE)=3 END IF KEY(NE,l)=KN+2 KEY(NE,6)=KN+1 KEY(NE,3)=KN KEY(NE,5)=INT(SEC(IDZ(NS)+2,NS,3)) KEY(NE,4)=KEY(NE,5)+1 KEY(NE,2)=INT(SEC(2*IDZ(NS)+3,NS,3)) NE=NE+1 DO 1 5 1=2,MZ DO 1 0 J = l , 8 KEY(NE,J)=KEY(NE-l,J)+2 CONTINUE N0TE(NE)=2 IF(KK.EQ.2) N0TE(NE)=3 , NE=NE+1 CONTINUE IF(KK.EQ.1)KEY(NE-1, l)=5*NEA+2 C C Fill across the s e c t i o n . C C JFN=MZ DO 4 0 1=2,MX IF(KK.EQ.1)JFN=JFN-1+M0D(I,2) DO 2 0 J = 1 , J F N I F ( M 0 D ( I , 2 ) . E Q . O ) THEN K E Y ( N E , 3 ) =KEY ( N E - J F N , 3 ) KEY(NE,6)=KEY(NE-JFN,5) KEY(NE,1)=KEY(NE-JFN,2) KEY(NE,5)=KEY(NE,6)-1 KEY(NE,4)=KEY(NE,1)-1 KEY(NE,2)=KEY(NE,4)-1 ELSE K E Y ( NE, 1 ) =KEY ( N E - J F N , 1 ) KEY(NE,6)=KEY(NE-JFN,4) KEY(NE,3)=KEY(NE-JFN,2) KEY(NE,5)=INT(SEC((IDZ(NS)*I+2*(MZ-JFN+J)),NS,3)) KEY(NE 4 ) K E Y ( N E 5)+1 KEY(NE|2)=INT(SEC((IDZ(NS)*(I+1)+2*(MZ-JFN+J)+1),NS,3)) END I F N0TE(NE)=2 IF(KK.EQ.2) N0TE(NE)=3 NE=NE+1 CONTINUE CONTINUE CONTINUE = 20 40 100 C C SECTIONS 5 & 7 C MX=2*NED DO 2 0 0 KK=1,2 I F ( K K . E Q . l ) THEN MZ=NWA NS=5 NPM=1 KN=7*NEA+2 N0TE(NE)=2 ELSE MZ=NEO Appendix C. FEM Software 110 115 C C Fill C NS=7 NPM=-1 KN=2+5*NEA+2*(NEW+2*NE0) N0TE(NE)=3 END IF KEY(NE,2)=KN+2*NPM KEY(NE,5)=KN+1*NPM KEY(NE,3)=KN KEY(NE,6)=INT(SEC(IDZ(NS)+2,NS,3)) KEY(NE,4)=KEY(NE,6)+1 KEY(NE,1)=INT(SEC(2*IDZ(NS)+3,NS,3)) NE=NE+1 DO 115 1=2,MZ DO 110 J=l,8 KEY(NE,J)=KEY(NE-l,J)+2 CONTINUE N0TE(NE)=2 IF(KK.EQ.2) THEN N0TE(NE)=3 KEY(NE,3)=KEY(NE-1,3)-2 KEY(NE,5)=KEY(NE-l,5)-2 KEY(NE,2)=KEY(NE-l,2)-2 END IF NE=NE+1 CONTINUE IF(KK.EQ.l) KEY(NE-l,2)=5*NEA+2 across the s e c t i o n . JFN=MZ DO 140 1=2,MX DO 120 J=1,JFN IF(M0D(I,2).EQ.O) THEN KEY(NE,3)=KEY(NE-JFN,3) KEY (NE, 5) =KEY (NE-JFN, 6 ) KEY(NE,2)=KEY(NE-JFN,1) KEY(NE,6)=KEY(NE,5)-1 KEY(NE,4)=KEY(NE,2)-1 KEY(NE,1)=KEY(NE,4)-1 ELSE KEY(NE,2)=KEY(NE-JFN,2) KEY(NE,5)=KEY(NE-JFN,4) KEY(NE,3)=KEY(NE-JFN,1) KEY(NE,6)=INT(SEC((IDZ(NS)*I+2*(MZ-JFN+J)),NS,3)) KEY(NE,4)=KEY(NE,6)+1 KEY(NE,1)=INT(SEC((IDZ(NS)*(I+1)+2*(MZ-JFN+J)+1),NS,3)) END IF N0TE(NE)=2 IF(KK.Eq.2) N0TE(NE)=3 NE=NE+1 CONTINUE CONTINUE CONTINUE 120 140 200 C C Sections 6&8 C MX=2*NED DO 300 KK=1,2 IF (KK.Eq.l) THEN MZ=NEW NS=6 KN=5*NEA+2 KEY(NE,4)=7*NEA+1 Appendix C. FEM Software 210 215 C C Fill C KEY(NE,1)=7*NEA K0TE(NE)=2 ELSE MZ=NOB NS=8 KN=2+6*NEA KEY(NE,4)=7*NEA+2+NDWA KEY(NE,1)=KEY(NE,4)+1 N0TE(NE)=3 END IF KEY(NE,2)=KN KEY(NE,5)=INT(SEC(IDZ(NS)+1,NS,3)) KEY(NE,6)=KEY(NE,5)+1 KEY(NE,3)=INT(SEC(2*IDZ(NS)+1,NS,3)) NE=NE+1 DO 215 1=2,MZ DO 210 J=l,8 KEY(NE,J)=KEY(NE-l,J)+2 CONTINUE N0TE(NE)=3 IF (KK.EQ.l) THEN KEY (NE, 4) =KEY (NE-1,4)-2 KEY(NE,1)=KEY(NE-1,1)-2 KEY(NE,2)=KEY(NE-l,2)-2 N0TE(NE)=2 IF(I.EQ.2) KEY(NE,2)=7*NEA ELSE IF(I.EQ.2) KEY(NE,2)=KEY(NE-1,1) END I F NE=NE+1 CONTINUE across the s e c t i o n . JFN=MZ NOF=JFN DO 240 1=2,MX DO 220 J=1,JFN IF(M0D(I,2).EQ.O) THEN KEY(NE,2)=KEY(NE-NOF,1) KEY(NE,5)=KEY(NE-NOF,6) KEY(NE,3)=KEY(NE-NOF,3) KEY(NE,6)=KEY(NE,3)+1 KEY(NE,l)=KEY(NE,3)+2 KEY(NE,4)=KEY(NE,5)+1 ELSE KEY (NE, 2 ) =KEY (NE-NOF ,3) KEY ( NE, 4 ) =KEY ( NE-NOF, 6 ) KEY (NE, 1) =KEY (NE-NOF, 1) KEY(NE,5)=INT(SEC((IDZ(NS)*I+2*(J-1)+1),NS,3)) KEY(NE,6)=KEY(NE,5)+1 KEY(NE,3)=INT(SEC((IDZ(NS)*(I+1)+2*(J-1)+1),NS,3)) END I F N0TE(NE)=2 IF(KK.EQ.2) N0TE(NE)=3 NE=NE+1 CONTINUE CONTINUE CONTINUE 220 240 300 C C Sections 2&4 C MX=2*NEI Appendix C. FEM Software 310 315 C C Fill C 320 340 400 DO 400 KK=1,2 IF (KK.Eq.l) THEN MZ=NEW NS=2 KN=5*NEA+2 KEY(NE,4)=KN+1 N0TE(NE)=2 ELSE MZ=N0B NS=4 KEY(NE,4)=2+7*NEA+NDWA KN=2+6*NEA N0TE(NE)=3 END IF KEY(NE,1)=KN KEY(NE,2)=KEY(NE,4)+1 KEY(NE,6)=INT(SEC(IDZ(NS)+1,NS,3)) KEY(NE,5)=KEY(NE,6)+1 KEY(NE,3)=INT(SEC(2*IDZ(NS)+1,NS,3)) NE=NE+1 DO 315 1=2,MZ DO 310 J=l,8 KEY(NE,J)=KEY(NE-l,J)+2 CONTINUE N0TE(NE)=2 IF (KK.EQ.2) THEN IF(I.Eq.2) KEY(NE,1)=KEY(NE-1,2) N0TE(NE)=3 END IF NE=NE+1 CONTINUE across the s e c t i o n . JFN=MZ N0F=JFN DO 340 1=2,MX DO 320 J=1,JFN IF(MOD(I,2).Eq.O) THEN KEY(NE,1)=KEY(NE-NOF,2) KEY(NE,6)=KEY(NE-NOF,5) KEY(NE,3)=KEY(NE-N0F,3) KEY(NE,5)=KEY(NE,3)+1 KEY(NE,2)=KEY(NE,3)+2 KEY(NE,4)=KEY(NE,6)+1 ELSE KEY(NE,2)=KEY(NE-NOF,2) KEY(NE,4)=KEY(NE-N0F,5) KEY(NE,1)=KEY(NE-NOF,3) KEY(NE,6)=INT(SEC((IDZ(NS)*I+2*(J-1)+1),NS,3)) KEY(NE,5)=KEY(NE,6)+1 KEY(NE,3)=INT(SEC((IDZ(NS)*(I+1)+2*(J-1)+1),NS,3)) END IF N0TE(NE)=3 IF(KK.Eq.l) N0TE(NE)=2 NE=NE+1 CONTINUE CONTINUE CONTINUE RETURN END Appendix C. FEM Software C.2 182 E N E R G Y Subroutine The subroutine ENERGY, which calculates the free energy of the system by integrating the square of the electric field and the osmotic pressure over the system, is presented here. This subroutine was used for Models I and III. The common block variables have the same definitions as found at the beginning of POWI.CYL. This integration procedure in ENERGY is identical to the procedure found in the subroutine SETMAT, which was presented above. c SUBROUTINE C C C C C C C C C C C C C ENERGY(ELEC.OSM) Subroutine t o c a l c u l a t e the f r e e energy of the double l a y e r based on the volume i n t e g r a l of the square of the e l e c t r i c f i e l d and t h e osmotic pressure term. OUTPUT: ELEC ; c o n t r i b u t i o n OSM ; c o n t r i b u t i o n from e l e c t r i c f i e l d from osmotic term. term. IMPLICIT REAL*8 (A-H.O-Z) REAL*8 L1(4),L2(4),L3(4) DIMENSION GWT(4),B(6) DATA NG/3/ DATA L1/0.5D0,0.5D0,0.D0,0.D0/ DATA L2/0.D0,0.5D0,0.5D0,0.D0/ DATA L3/0.5D0,0.D0,0.5D0,0.D0/ DATA GWT/0.3333333333333333D0.0.3333333333333333D0, CO.3333333333333333D0,0.3333333333333333D0/ DATA Ll/0.3333333333333333D0,0.6D0,0.2D0,0.2D0/ DATA L2/0.3333333333333333D0,0.2D0,0.6D0,0.2D0/ DATA L3/0.3333333333333333D0,0.2D0,0.2D0,0.6D0/ C COMMON /BLK1/ZND(5000,3),SBM(2500,21),S0R(5000),SEXP(20), (DSQL.FLG.ZM COMMON /BLK3/KEY(2500,6),N0TE(2500) COMMON /BLK7/NEA,NEI,NED,NEW,NEO,NALP,NWA,NOB,NEL,NOD,NWAI,NWAD,NOBI,NOBD COMMON /BLK4/DBS(6),DBT(6),DBR(6),DBA(6),AJAC(2,2),AJIN(2,2) l.DETJ,XE,ZE,RE,EM(3),DBX(6),DBZ(6) C GWT(1)=-27.D0/48.D0 GWT(2)=25.D0/48.D0 GWT(3)=GWT(2) GWT(4)=GWT(2) ELEC=0.D0 0SM=0.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NUMERICAL INTEGRATION SECTION CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Perform quadrature over t r i a n g u l a r elements. C C C EVALUATE CONTRIBUTION OF EACH POINT. C Appendix C. FEM Software VOL=O.DO DO 100 11=1,NG B(l)=Ll(II)*(2.D0*Ll(II)-l.D0) B(2)=L2(II)*(2.D0*L2(II)-1.D0) B(3)=L3(II)*(2.D0*L3(II)-1.DO) B(4)=4.D0*L1(II)*L2(II) B(5)=4.D0*L3(II)*L2(II) B(6)=4.D0*L1(II)*L3(II) DBS(1)=4.'DO*L1(II)-1.DO DBS(2)=0.D0 DBS(3)=1.DO-4.D0*L3(II) DBS(4)=4.D0*L2(II) DBS(5)=-4.D0*L2(II) DBS(6)=4.D0*(L3(II)-L1(II)) DBT(1)=0.D0 DBT(2)=4.D0*L2(II)-1.D0 DBT(3)=1.D0-4.D0*L3(II) DBT(4)=4.D0*L1(II) DBT(5)=4.D0*(L3(II)-L2(II)) DBT(6)=-4.D0*L1(II) C C CALCULATE CONTRIBUTION OF EACH ELEMENT FOR THE QUADRATURE POINT. C LFG=1 DO 50 NE=1,NEL CALL JACOB(NE,II,LFG) PE=0.D0 0S=O.D0 NSG=1 IF(N0TE(NE).EQ.2) THEN C C CALCULATE OSMOTIC TERM C NSG=-1 DO 45 1=1,6 PE=PE+B(I)*SOR(KEY(NE,I)) 45 CONTINUE 0S=(ZM*(DEXP(-CP*PE)-1.DO)+CP*(DEXP(ZM*PE)-1.DO))/(ZM**2*2.DO $ *CP) END IF C C CALCULATE FIELD TERMS C DPDR=O.DO DPDZ=O.DO DO 40 1=1,6 DPDR=DPDR+DBR(I)*SOR(KEY(NE,I)) DPDZ=DPDZ+DBZ(I)*SOR(KEY(NE,I)) 40 CONTINUE C C INTEGRATE C DO 60 1=1,6 OSM=OSM+OS*GWT(II)*DABS(DETJ)*B(I)*ZND(KEY(NE,I),2) ELEC=ELEC+(DPDR**2+DPDZ**2)*GWT(II)*DABS(DETJ)*NSG $ *B(I)*ZND(KEY(NE,I),2)*EM(N0TE(NE))/2.D0 V0L=V0L+GWT(11)*DABS(DETJ)*B(I)* ZND(KEY(NE,I),2) 60 CONTINUE 50 CONTINUE 100 CONTINUE PRINT*,'VOLUME = '.VOL RETURN END Appendix C. FEM Software C.3 184 T y p i c a l Input File A typical input file that may be used for each of model programs is shown below. The numerical data is supplied as is with a template description following. 32 10 10 5.DO 2.5D0 0.D0 0.D0 4.DO 80. DO l.D-2 15.DO 1.D0 1.D0 1.D0 12.DO 0 2.DO INPUT FILE FOR 10 10 30 80.DO 20. DO 1.D0 1.D0 2.DO 2.DO 1.D0 1.D0 1.D0 1.D0 12.DO -99.DO POWI.CYL, POWI.SP, PH20.SRC NEA,NEI,NED,NWA,NOB,MAIXT,NALP RA.RB.ZO.XI ORX.ORZ.CP.ZH EM(3),PSI ERR('/.),ZW,XD2,SKL FZWA,FZWP,FZOP,FZOB FCZ,FCXI,FCXD,SOW NQ.PINF (715) (4D8.2) (4D8.2) (4D8.2) (4D8.2) (4D8.2) (4D8.2) (1X,I2,D8.2) S0W=-99., uncharged 0/W i n t e r l a c e . S0W=? , p o t e n t i a l on i n t e r f a c e . NQ =0 , constant p o t e n t i a l p a r t i c l e . NQ =1 , constant chargel p a r t i c l e . Appendix D Supplemental Results This appendix contains results that are referenced i n chapter 4. T h e y are presented here for convenience sake. 185 Appendix D. Supplemental Results Table D . l : Force of Repulsion (/**) Between Two Cylindrical Particles (/ca Embedded in an Uncharged O / W Interface 2a] K[R- 45° 1 2 3 4 5 6 7 8 9 10 11 12 13 15 17.4 19 25 30 35 40 50 60 70 80 water phase / * 90° 5.347E00 2.203E00 8.602E-1 3.248E-1 1.20213-1 4.392E-2 1.599E-2 5.823E-3 2.125E-3 7.787E-4 2.877E-4 1.078E-4 4.143E-5 7.259E-6 1.696E-6 9.953E-7 3.813E-7 2.215E-7 1.381E-7 9.027E-8 4.259E-8 2.205E-8 1.218E-8 7.047E-9 2.601E-9 2.673E00 1.104E00 4.329E-1 1.646E-1 6.149E-2 2.282E-2 8.485E-3 3.182E-3 1.211E-3 4.722E-4 1.907E-4 8.117E-5 3.716E-5 1.048E-5 3.954E-6 2.593E-6 9.722E-7 5.455E-7 3.332E-7 2.150E-7 1.002E-7 5.189E-8 2.885E-8 1.687E-8 6.404E-9 100 Degrees re: er to contact angle 135° 45° 1.913E-3 7.431E-4 3.022E-4 1.307E-4 6.123E-5 3.153E-5 1.794E-5 1.129E-5 7.730E-6 5.665E-6 4.367E-6 3.493E-6 2.870E-6 2.045E-6 1.445E-6 1.174E-6 5.947E-7 3.662E-7 2.373E-7 1.600E-7 7.924E-8 4.281E-8 2.458E-8 1.478E-8 5.887E-9 1.607E1.510E1.366E1.206E1.050E9.085E7.852E6.798E5.906E5.156E4.524E3.991E3.538E2.823E2.203E1.889E1.134E7.861E5.658E4.188E2.438E1.499E9.547E6.226E2.763E- oil phase f* 90° 618E-4 228E-3 092E-3 232E-3 991E-3 617E-3 236E-3 896E-3 608E-3 372E-3 178E-3 019E-3 ,885E-4 896E-4 .251E-4 .449E-4 .597E-4 780E-4 .274E-4 .419E-5 .502E-5 .421E-5 .214E-5 472E-5 833E-6 10, V'o 135 c 2.646E-3 2.205E-3 1.857E-3 1.580E-3 1.358E-3 1.177E-3 1.030E-3 9.064E-4 8.028E-4 7.151E-4 6.402E-4 5.758E-4 5.201E-4 4.289E-4 3.461E-4 3.026E-4 1.920E-4 1.376E-4 1.017E-4 7.706E-5 4.666E-5 2.984E-5 1.979E-5 1.346E-5 6.543E-6 Appendix D. Supplemental Results 187 o Figure D.l: Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (V'o = 2) Embedded in an Uncharged 0/W Interface 188 Appendix D. Supplemental Results o K[R-2Q] Figure D.2: Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (V'o = 2) Embedded in an Uncharged 0 / W Interface Figure D.3: Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (ip = 2) Embedded in an Uncharged 0/W Interface Q Appendix D. Supplemental Results 190 o o Figure D.4: Parallel Interaction Force Acting Between Two Identical Cylinders with a Constant Surface Potential (xp = 4) Embedded in an Uncharged 0/W Interface 0 Figure D.5: Normal Force Acting on Constant Potential Cylinders (V'o = 4) Embedded in an Uncharged 0 / W Interface Appendix D. Supplemental Results Figure D.6: Parallel Force Acting between Two Constant Potential Cylinders (-00 Embedded in a Constant Potential 0 / W Interface Appendix D. Supplemental Results Figure D.7: Parallel Force Acting between Two Constant Potential Cylinders (xp Embedded in a Constant Potential 0/W Interface Appendix D. Supplementad Results Figure D.8: Parallel Force Acting between Two Constant Potential Cylinders (tp Embedded in a Constant Potential 0/W Interface Appendix D. Supplemental Results tn o o 4R-2a] Figure D.9: Parallel Force Acting between Two Constant Potential Cylinders (V'o Embedded in a Constant Potential 0 / W Interface Appendix D. Supplemental Results Figure D.10: Parallel Force Acting between Two Constant Potential Cylinders (V'o Embedded in a Constant Potential 0 / W Interface Appendix D. Supplemental Results m o 6 - 135. ICQ = 15. interface o = - 1 0 += 1 A = x = 2 O * "o D Q_ g o o o o 0.0 1.0 2.0 3.0 4.0 5.0 Ac[R-2a] Figure D.ll: Parallel Force Acting between Two Constant Potential Cylinders (V>i Embedded in a Constant Potential 0/W Interface Appendix D. Supplemental Results Figure D.12: Parallel Force Acting between Two Constant Potential Cylinders (ip Embedded in a Constant Potential 0 / W Interface Appendix D. Supplemental Results Figure D.13: Parallel Force Acting between Two Constant Potential Cylinders (V'o Embedded in a Constant Potential 0/W Interface 200 Appendix D. Supplemental Results KQ =2. tea = 5. Contact Angle a = 45 o = 90 A = 135 Contact Angle o = 45 o = 90 A = 135 o <r 0.12S 0.600 1.075 1.550 2.025 2.500 0.25 US 1.60 2.05 2.50 /cd *ca = 10. o E Contact Angle o = 45 o = 90 o = 135 2 ° 0.5 0.9 1.3 1.7 2.1 2.5 Kd Figure D.14: Normal Force Acting on a Constant Potential Sphere (Vo ' = 4) Embedded in an Uncharged 0/W Interface Appendix E Nomenclature a radius of curvature of particle (m) B Bond number; equation (6.2) C electrolyte concentration (molar) d separation from particle surface to median plane (m) e elementary charge (C) E electric field (V/m) r elemental load vector r reduced force; equations (2.33) and (2.34) r scaled reduced force; equations (4.1) and (4.2) Eelec electrostatic normal force (N/m) 9 acceleration due to gravity (m/s ) G electrostatic free energy (J) G* reduced free energy; equations (2.36) and (2.37) G** scaled reduced free energy; equation (4.36) h height of point on particle surface from line joining particle centres (m) H separation of particle surfaces at height h (m) k Boltzmann's constant (J/K) K 2 e elemental stiffness matrix Kl contribution to K from Laplacian operator K contribution to K from charge term e e B 201 Appendix E. Nomenclature 202 L equation (6.3) M number of nodes in element n unit normal vector N Avagadro's number (molecules/mole) a N basis function of i R separation of centres of particles (m) Rn residual T absolute temperature (K) v + th node reduced parallel plate energy/unit area ; equation (5.15) in [38] ( J / m ) 2 W weighting function x,z,r displacement co-ordinates (m) a:,etc. reduced displacement co-ordinates; equation (2.7) x,etc. alernative reduced displacement co-ordinates; equation (2.31) z± ion valences L Greek Symbols a under/over relaxation parameter P immersion angle of particle (6 — 180°) 7 interfacial tension (N/m) r Maxwell stress tensor 6 tolerance parameter; equation (2.62) e dielectric constant ; subscripts: w-water, o-oil, p-particle e three phase contact angle (°) double layer thickness (m); equation (2.8) A ratio of ion valences; equation (2.9) Appendix E. Nomenclature A factor from series expansion of charge term; equation (2.61) v± local ion concentrations (number/cm ) n osmotic pressure ( N / m ) p charge density ( C / m ) 3 2 3 constant surface charge (C) a* reduced surface charge T reduced osmotic pressure; equation (2.34) $ e $ elemental parameter function vector of field variables value of field variable at i th node electrostatic potential (V) reduced electrostatic potential; equation (2.6) V'oo surface potential of an isolated particle 4>o constant surface potential of a particle constant surface potential of oil/water interface potential function in the particle phase potential function in the oil phase potential function in the water phase ft elemental domain
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Electrostatic interactions between interfacial particles
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Electrostatic interactions between interfacial particles Lyne, Michael Peter 1989
pdf
Page Metadata
Item Metadata
Title | Electrostatic interactions between interfacial particles |
Creator |
Lyne, Michael Peter |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | The Finite Element Method (FEM) is employed, using Galerkin weighted residuals, to simulate the electrostatic interactions between interfacial particles in two separate models. The three phase problem of two identical, parallel, infinitely long, solid cylinders at a given separation and each embedded to the same extent in an oil/water interface is considered. It is revealed that due to the asymmetric charge distribution surrounding an interfacial particle, a force (normal to the interface) acts on the particle in the direction of the electrolyte phase. The electrostatic forces acting on the particles normal and parallel to the oil/water interface are calculated and it is found that under certain conditions they may be of comparable magnitude. The Debye-Hiickel form of the Derjaguin method is also applied to this model; this method showed good agreement with the low potential numerical results for the parallel force when the oil/water interface is assumed to be uncharged. Interactions between the cylinders on a constant potential interface are also considered and it is revealed that for some circumstances the magnitude of these interactions will be significantly different from the corresponding interactions on an uncharged interface. For practical situations, however, increases in the interaction force were not large enough to account for a 10⁶kT change in the free energy associated with the interaction. The more complex interaction between particles in a dense, two-dimensional structure on an oil/water interface is simulated using a cell model approach. The normal forces acting on the spheres and the electrostatic free energy change associated with forming the structure are calculated. The normal force per unit of particle perimeter in contact with the oil/water interface is shown to be comparable to that acting on cylinders while the free energy calculations are only accurate enough for order of magnitude estimations. Additionally, to test the effectiveness of the numerical approach, the familiar case of two identical interacting spheres completely immersed in a single electrolyte is modeled for both constant potential and constant charge boundary conditions on the sphere surfaces. The implicit assumption inherent in the Derjaguin method when it is applied to constant charge interactions (i.e. zero potential gradient inside the particle) is examined and found to cause errors of > 10% in cases of low ka and large particle dielectric constants. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0058823 |
URI | http://hdl.handle.net/2429/27899 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1989_A7 L96.pdf [ 9.25MB ]
- Metadata
- JSON: 831-1.0058823.json
- JSON-LD: 831-1.0058823-ld.json
- RDF/XML (Pretty): 831-1.0058823-rdf.xml
- RDF/JSON: 831-1.0058823-rdf.json
- Turtle: 831-1.0058823-turtle.txt
- N-Triples: 831-1.0058823-rdf-ntriples.txt
- Original Record: 831-1.0058823-source.json
- Full Text
- 831-1.0058823-fulltext.txt
- Citation
- 831-1.0058823.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0058823/manifest