UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
SURFACE TENSION AND COMPUTER SIMULATION OF POLYATOMIC FLUIDS By JAMES MITCHELL HAILE A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1976 ACKNOWLEDGEMENTS It is a pleasure to express my gratitude to those who have freely contributed to this work through their instruction, guidance, and advice. Keith Gubbins initiated the research reported herein and enthu siastically stimulated and supported its development, as well as my own professional growth. Dr. Gubbins consistently provided an extraordinary environment for learning, provided a strong example of the scientific method, and exposed me to numerous knowledgeable scientists and engineers on both sides of the Atlantic. In addition, he expended considerable effort in obtaining the financial support, travel and computer funds which made this work possible. Bill Streett generously allowed me to spend several months in his laboratory at the U.S. Military Academy and patiently taught me molecular dynamics. He obtained the copious amounts of computer time used in the molecular dynamics work reported here and kept the program running in my absence. Many of the ideas for presenting the molecular dynamics results came to light in discussions with Colonel Streett. Further, I am grateful to Colonel and Mrs. Street for the hospitality extended to me during my visits to West Point. John O'Connell, University of Florida, continually inspired me through openended questioning concerning classical and statistical thermodynamics, science, engineering, and, most importantly, the character of life. Chris Gray, University of Guelph, instructed me in spherical trigonometry, spherical harmonic expansions, Racah algebra, etc., thereby developing in me a healthy respect for the physicist's view of applied science. I have benefited greatly from countless discussions with my colleague ChorngHorng Twu on various aspects of thermodynamics, statistical mechanics, numerical methods, and Chinese cooking. Skren Toxvaerd, University of Copenhagen, contributed much valuable advice on the theory and associated calculations for fluid interfaces. Thanks are also due Dr. Toxvaerd for providing a copy of his computer program for calculating the vaporliquid interfacial density profile for LennardJones fluids. Peter Egelstaff allowed me to spend several months in the stimulating atmosphere of the Physics Department at the University of Guelph. I am grateful to the faculty and staff for their hospi tality and for the large amount of NOVA 2 computer time made available to me. I am especially thankful to Dan Litchinsky for useful advice on NOVA 2 software and to Ross McPherson for timely hardware support on the NOVA. I am also indebted to ShienShion Wang for spending many hours in teaching me the Monte Carlo method. Dick Dale and Ron Franklin of the Engineering Information Office, University of Florida, gave timely and enthusiastic photographic tech nical assistance in producing the filmed animation of molecular dynamics simulations. Larry Mixon in the Northeast Regional Data Center, Univer sity of Florida, provided valuable software support in developing the filmed animation technique. I am grateful to Dr. J. W. Dufty for serving on the supervisory committee. I would also like to remember Dr. T. M. Reed who was an original member of the committee and who strongly encouraged me in the initial phases of the research reported here. P.S.Y. Cheung, T. Keyes, C. G. Gray and R. L. Henderson, W. B. Street, and S. Toxvaerd kindly provided manuscripts of their work prior to publica tion. Mrs. J. Ojeda, University of Florida, performed the remarkably excellent typing of the manuscript. Finally, I am grateful to Tricia who, in addition to all the usual annoyances with which wives of Ph.D. students are plagued, quietly endured our being separated for the greater part of the last year and a half of this work. The three dimensional drawings presented in Chapter 7 were done on a Gould 5100 electrostatic plotter driven by the IBM 370/165 at the Northeast Regional Data Center, University of Florida. The associated software was the SYMVU Computer Graphics Program, Version 1.0, of the Laboratory for Computer Graphics and Spatial Analysis, Harvard University. I thank the Petroleum Research Fund (administered by the American Chemical Society) and the National Science Foundation for financial support of this study. TABLE OF CONTENTS ACKNOWLEDGEMENTS................................................ LIST OF TABLES.................................................... LIST OF FIGURES................................................ .. KEY TO SYMBOLS................................................... ABSTRACT......................................................... CHAPTERS: 1 INTRODUCTION............................................ 1.1 Theory of Surface Properties..................... 1.2 Computer Simulation Methods ...................... 1.3 Outline of Dissertation.......................... 2 THEORY OF SURFACE TENSION.............................. 2.1 General Expressions for Surface Tension of Polyatomic Fluids................................ 2.2 General First Order Perturbation Theory for Surface Tension.................................. 2.3 Perturbation Theory for Surface Tension using a Pople Reference................................ 2.4 Fowler Model Expressions for Perturbation Terms Y2A' Y2B' 3A' 3B in Pople Expansion ............ 2.5 Superficial Excess Internal Energy from the Pade Perturbation Theory for Surface Tension..... 3 NUMERICAL CALCULATIONS OF SURFACE TENSION.............. F F F F 3.1 Evaluation of 2A' 2B' 3A' and 3B ............ 3.2 Surface Tension Calculations for Model Fluids.... Page ii ix xv xxii xxix 1 6 9 10 13 13 18 23 29 32 34 35 40 TABLE OF CONTENTS (Continued) CHAPTERS: Page 3.3 Calculation of the Superficial Excess Internal Energy for Model Fluids........................... 46 3.4 Surface Tension Calculations for Real Fluids...... 48 3.5 Correlation of Surface Tension for Pure Poly atomic Liquids.................. ... ............. 58 4 VAPORLIQUID DENSITYORIENTATION PROFILES................ 72 4.1 First Order Perturbation Theory for p(z1wl) ....... 72 4.2 Calculations of p(zl w) for Overlap and Dis persion........................................... 78 5 MONTE CARLO SIMULATION OF MOLECULAR FLUIDS ON A MINICOMPUTER.......................................... 94 5.1 Introduction...................................... 94 5.2 Monte Carlo Method for Nonspherical Molecules...... 96 5.3 Description of the Minicomputer System............ 101 5.4 Monte Carlo Program for the NOVA.................. 102 5.5 Comparison of NOVA Results with FullSize Computer Results.................................. 107 5.6 Conclusions............................. ......... 112 6 MOLECULAR DYNAMICS METHOD FOR AXIALLY SYMMETRIC MOLECULES................................................ 114 6.1 Introduction ................... .................. 114 6.2 Expressions for the Force and Torque for Axially Symmetric Molecules............................. 117 6.3 Method of Solution of the Equations of Motion and the Molecular Dynamics Algorithm.............. 129 6.4 Evaluation of Pair Correlation Functions.......... 136 6.5 Equilibrium Properties from the g Z 2m(rl2)....... 146 1l2m ).. TABLE OF CONTENTS (Continued) CHAPTERS: Page 7 MOLECULAR DYNAMICS RESULTS.............................. 158 7.1 Potential Models.................................. 158 7.2 Equilibrium Properties.................. ......... 171 7.3 Spherical Harmonic Coefficients,g 1 2m(rl2)....... 191 7.4 Angular Pair Correlation Function.................. 215 7.5 SiteSite Pair Correlation Functions.............. 238 7.6 Filmed Animation of Molecular Motions............. 250 8 CONCLUSIONS........................................... 259 8.1 Theory for Surface Tension of Polyatomic Fluids... 259 8.2 Theory for the Interfacial DensityOrientation Profile of Polyatomic Fluids....................... 261 8.3 Computer Simulation of Polyatomic Fluids........... 262 APPENDICES: A EXPRESSIONS FOR THE ANGLE AVERAGES IN EQUATIONS (34) TO (37)................................................ 267 B COORDINATE TRANSFORMATION AND INTEGRATION OVER EULER ANGLES TO OBTAIN EQUATIONS (289) AND (290)............ 270 B.1 Choice of Euler Angles............................ 270 B.2 Evaluation of Integral I ......................... 273 z C MODELS FOR ANISOTROPIC POTENTIALS OF LINEAR MOLECULES... 277 F F F F D EXPRESSIONS FOR y2A' Y2B' 3A AND y3B FOR VARIOUS ANISOTROPIC POTENTIALS FOR AXIALLY SYMMETRIC MOLECULES.. 283 E THE INTEGRALS KY(k'V";nn'n") AND LY (;nn') ............. 288 F EXPRESSIONS FOR THE SPHERICAL HARMONIC COEFFICIENTS g 1 2m(ri 2) IN EQUATION (677).......................... 292 G THE INTEGRAL I USED TO CALCULATE THE ANGULAR COR RELATION PARAMETR G2 FOR QUADRUPOLES................... 294 vii TABLE OF CONTENTS (Continued) APPENDICES: Page H VALUES FOR THE g. m(rl2) COEFFICIENTS................. 297 12 I VALUES FOR THE J INTEGRALS ......................... 328 n J VALUES OF THE SITESITE CORRELATION FUNCTIONS........... 349 K VALUES OF THE INTEGRAL H126.......................... 355 11 LITERATURE CITED........ ............................ ........... 357 BIBLIOGRAPHY ........................................................ 365 BIOGRAPHICAL SKETCH....... ....................................... 376 viii LIST OF TABLES Table Page 1 Examples of Macroscopic and Microscopic Interfacial Properties Related by Equations of the Form (11)......... 5 2 Test of the GibbsHelmholtz Equation in the Fowler Model Perturbation Theory for LennardJones plus Quadrupole Fluids ......................................... 47 3 Potential Parameter Values used in Calculating Surface Tension... ................................................. 57 4 Values for the Parameters al and a2 in the Surface Tension Correlation of Equation (361).................... 62 5 Equilibrium Properties in the Form of Ensemble Averages... 100 6 Approximate Number of Monte Carlo Configurations Generated per Hour on the NOVA 2........................... 106 7 Comparison of NOVA and CDC Results for Property Values of LennardJones + Quadrupole Model Fluid. kT/E = 0.719, pC3 = 0.80, Q/(ca5)1/2 = 1................................ 109 8 Expressions for in Terms of g 1 m(r2) a 1lj2 1g 2m l2 for Various Model Potentials.............................. 150 9 Expressions for the Configurational Energy in Terms of g 2m(r12) for Various Model Potentials................... 151 10 Expressions for the Pressure in Terms of g. (r12) for Various Model Potentials...................1.2............ 152 11 Expressions for the Fowler Model Surface Tension in Terms of g I 2m(rl2) for Various Model Potentials........ 153 12 Expressions for the Fowler Model Surface Excess Internal Energy in Terms of g. (r1) for Various Model Potentials.......... 1.2 1................................. 154 13 Expressions for the Mean Squared Force in Terms of g 1 2m(r12) for Various Model Potentials................. 155 RR 12 Table Page 14 Expressions for the Mean Squared Torque in Terms of S1 2m(r12) for Various Model Potentials................. 156 15 Expressions for the Angular Correlation Functions GL in Terms of g 1 2m(r2 ) for Various Model Potentials......... 157 16 Primary Orientations for Pairs of Linear Molecules....... 162 17 Property Values of a LennardJones plus Quadrupole Fluid Obtained in this Work and Compared with those given by Berne and Harp................................... 177 18 Equilibrium Properties for LennardJones plus Quadru pole Fluid at pa3 = 0.85, Q/(EC5)1/2 = 1/2............... 181 19 Equilibrium Properties for LennardJones plus Quadru pole Fluid at po3 = .931, Q/(Ec5)1/2 = 0.707............. 182 20 Equilibrium Properties for LennardJones plus Quadru pole Fluid at po3 = 0.85, Q/(Eo5)1/2 = 1.0............... 183 21 Anisotropic Contributions to Equilibrium Properties for LennardJones plus Quadrupole Fluid at po3 = 0.85, Q/i(E5)1/2 = 1/2 ........................................ 188 22 Anisotropic Contributions to Equilibrium Properties for LennardJones plus Quadrupole Fluid at po3 = .931, Q/(E~ 5)1/2 = 0.707....................................... 189 23 Anisotropic Contributions to Equilibrium Properties for LennardJones plus Quadrupole Fluid at po3 = 0.85, Q/(Ea 5)1/2 = 1.0......................................... 190 24 Equilibrium Properties for LennardJones plus Overlap Fluid at po3 = 0.85, 6 = 0.10............................ 192 25 Equilibrium Properties for LennardJones plus Overlap Fluid at po3 = 0.85, 6 = 0.30............................ 193 26 Effect of Potential Model and State Condition on the First Peak Height of the g 12m Coefficients............. 201 27 Comparison of Molecular Dynamics Results for Equilibrium Properties with Values Obtained from S12m J Integrals for LennardJones plus Quadrupole Fluid with po3 = 0.85, kT/c = 1.277, Q/(Eo5)]/2 = 0.5.........210 Table Page 28 Comparison of Molecular Dynamics Results for Equilib rium Properties with Values Obtained from 2^ m 1i2m J Integrals for LennardJones plus Quadrupole n 3 Fluid with p3 = 0.85, kT/E = 1.294, Q/(Ea5) /2 = 1.0... 211 29 Comparison of Molecular Dynamics Results for Equilib rium Properties with Values Obtained from 1 2m J Integrals for LennardJones plus Overlap Fluid with po3 = 0.85, kT/E = 1.291, 6 = 0.10.................. 212 30 Comparison of Molecular Dynamics Results for Equilib rium Properties with Values Obtained from k1 2m J Integrals for LennardJones plus Overlap Fluid with po3 = 0.85, kT/E = 1.287, 6 = 0.30................. 213 31 Range of Values for Orientational Contributions to Property Integrands for Quadrupole and Overlap Fluids... 243 Cl Expressions for the Expansion Coefficients E for Various Interaction Potentials for Linear Molecules..... 279 C2 Expressions for Anisotropic Potential Models in the Intermolecular Frame of Figure 32....................... 280 C3 Expressions for Anisotropic Potential Models in the Intermolecular Frame, using Y rather than j ........... 281 C4 Derivatives of Various Anisotropic Potentials for Evaluating the Force and Torque from Equations (625) and (634) .............................................. 282 D1 Expressions for yF for Various Anisotropic Potentials for Linear Molecues ..................................... 284 F D2 Expressions for y for Multipole Potentials for 3A Linear Molecules ........................................ 285 F D3 Expressions for y for Various Anisotropic Potentials for Axially Symmetric Molecules......................... 286 D4 Expressions for y3B for Multipole Potentials for Linear Molecules......................................... 287 El The Integrals KY(9'k";nn'n") for Pure Fluids........... 289 E2 The Integrals LY(;nn') for Pure Fluids ................. 290 Table Page E3 The Constants in Equation El............................ 291 G1 The Integral I554 for Pure Fluids....................... 296 HI Values of g000(rl2) g400(r1l) for LennardJones plus Quadrupole Fluid at kT/E = 1.277, po3 = 0.85, Q/(Eo5)1/2 = 0.5........................................ 298 H2 Values of g420(r12) g442(rl2) for the Fluid of Table HI1................................................ 300 H3 Values of g443(r12) g660(rl2) for the Fluid of Table HI................................................ 302 H4 Values of g000(rl2) g840(rl2) for LennardJones plus Quadrupole Fluid with kT/c = 0.765, po3 = 0.931, and Q/(Eo5 1/2 = 0.707...................................... 304 H5 Values of g420(rl2) g442(r12) for the Fluid of Table H4 ................................................ 306 H6 Values of g443(r12) g660(r12) for the Fluid of Table H4................................................ 308 H7 Values of g000(rl2) g400(rl2) for LennardJones plus Quadrupole Fluid with kT/E = 1.294, po3 = 0.85, and Q/(Eo5)1/2 = 1.0........................................ 310 H8 Values of g420(rl2) 8442(r12) for the Fluid of Table H7 ................................................ 312 H9 Values of g443(r12) g640(rl2) for the Fluid of Table H7................................................ 314 H10 Values of g000(rl2) g400(rl2) for LennardJones plus Anisotropic Overlap Fluid with kT/E = 1.291, po3 = 0.85, and 6 = 0.10............................................ 316 H11 Values of g420(rl2) g442(r12) for the Fluid of Table H10 ............................................... 318 H12 Values of g443(r12) g660(r12) for the Fluid of Table H10 ............................................... 320 H13 Values of g000(rl2) g400(rl2) for LennardJones plus Anisotropic Overlap Fluid with kT/E = 1.287, po3 = 0.85, and 6 = 0.30............................................. 322 xii Table H14 Values of g420(r12) g442(r12) for the Fluid of Table H13................................................ H15 Values of g443(r12) g660(r12) for the Fluid of Table H13............................................... 000 222 II The Integrals J J for a LennardJones plus Quadrupole Fluid. pa3"= 0.85, kT/E = 1.277, Q/(Eo5)1/2 = 0.5........................................ The Integrals J4 n The Integrals J44 n The Integrals J620 n The Integrals J000 n Quadrupole Fluid. Q(EG5)1/2 = 0.707. The Integrals J400 n4 The Integrals J4 n The Integrals J6 n The Integrals J0 Quadrupole Fluia. Q/(EO5)1/2 = 1.0.. 400 The Integrals J400 n4 The Integrals J4 n The Integrals J620 n 440 J for the Fluid of Table II..... n 600 for the Fluid of Table II..... n 660 J660 for the Fluid of Table II..... n J2 for a LennardJones plus po3 = 0.931, kT/c = 0.765, .. .... ..440......... .... J for the Fluid of Table 15..... n S600 for the Fluid of Table 15..... n  j660 for the Fluid of Table 15..... n  J for a LennardJones plus po3n= 0.85, kT/E = 1.294, 440 J for the Fluid of Table 19..... n  j600 for the Fluid of Table 19..... n 640  640 for the Fluid of Table 19..... n 000 222 113 The Integrals J J for a LennardJones plus Anisotropic Overlap Fluid. pa3 = 0.85, kT/E = 1.291, 6 = 0.10................................................. 400 440 114 The Integrals J J440 for the Fluid of Table 113.... n n 115 The Integrals J J600 for the Fluid of Table 113.... n n 620 660 116 The Integrals J J660 for the Fluid of Table 113.... n n 000 222 117 The Integrals J 0 J2 for a LennardJones plus n n Anisotropic Overlap Fluid. po3 = 0.85, kT/E = 1.287, 6 = 0.30 .............................................. .. xiii Page 324 326 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 Table Page 400 440 118 The Integrals J 4 J for the Fluid of Table 117.... 346 n n 441 600 119 The Integrals J 600 for the Fluid of Table 117.... 347 n n 620 660 120 The Integrals J J6 for the Fluid of Table 117.... 348 n n Jl SiteSite Correlation Function for LennardJones plus Quadrupole Fluid with /c = 0.3292, kT/c = 1.277, po3 = 0.85, Q/(co5)1/2 = 0.5.............................. 350 J2 SiteSite Correlation Function for LennardJones plus Quadrupole Fluid with /C = 0.2955, kT/E = 0.765, PO3 = 0.931, Q/(ao5)1/2 = 0.707......................... 351 J3 SiteSite Correlation Function for LennardJones plus Quadrupole Fluid with i/a = 0.3292, kT/E = 1.294, po3 = 0.85, Q/(Ec 5)1/2 = 1.0....................... ..... 352 J4 SiteSite Correlation Function for LennardJones plus Anisotropic Overlap Fluid with /a = 0.3292, kT/E = 1.287, po3 = 0.85, 6 = 0.10 ..................................... 353 J5 SiteSite Correlation Function for LennardJones plus Anisotropic Overlap Fluid with /a = 0.3292, kT/E = 1.287, po3 = 0.85, 6 = 0.30..................................... 354 K1 The Integral H 1261 for Pure Fluids..................... 356 11 xiv LIST OF FIGURES Figure Page 1 Relation of Theory, Experiment, and Computer Simula tion in the Study of Liquids............................ 2 2 Variation of Fluid Density with Position through a Planar VaporLiquid Interface........................... 4 3 Two Possible Values for the Maximum zi Value for a Pair of Molecules in the Fowler Model Interface............... 22 4 Fowler Model Surface Tension for a Fluid of Axially Symmetric Molecules Interacting with LennardJones plus Dipole Model Potential. po3 = 0.85, kT/E = 1.273 ............................................ 42 5 Fowler Model Surface Tension for a Fluid of Axially Symmetric Molecules Interacting with LennardJones plus Dipole Model Potential. pa3 = 0.45, kT/E = 2.934 ............................................ 43 6 Fowler Model Surface Tension for a Fluid of Axially Symmetric Molecules Interacting with LennardJones plus Quadrupole Model Potential. po3 = 0.85, kT/E = 1.273............................................ 44 7 Fowler Model Surface Tension for Fluids of Axially Symmetric Molecules Interacting with LennardJones plus Various Anisotropic Potentials. pa3 = 0.85, kT/E = 1.273 ............................................ 45 8 Corresponding States Plot for Surface Tension of Simple Liquids.......................................... 52 9 Surface Tension for CO2 Comparing Perturbation Theory Calculations with Experimental Values.................... 55 10 Surface Tensions for C2H2 and HBr Comparing Perturba tion Theory Calculations with Experimental Values....... 56 11 Test of Surface Tension Correlation for CO2............. 64 12 Test of Surface Tension Correlation for Acetic Acid..... 65 13 Test of Surface Tension Correlation for Methanol......... 66 Figure Page 14 Comparison of Surface Tensions Calculated from the Correlation with Experimental Values for Several Polyatomic Liquids...................................... 68 15 Comparison of Surface Tensions Calculated from the Correlation with Experimental Values for Several Polyatomic Liquids...................................... 69 16 Test of Surface Tension Correlation for nHexane and nOctane................................................ 70 17 Comparison of Surface Tensions Calculated from the Correlation with Experimental Values for Several Hydrocarbons... ......................................... 71 18 Interfacial Density Profile for LennardJones Fluid..... 80 19 Interfacial DensityOrientation Surface for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Dispersion Model Potential. kT/E = 0.85, K = 0.25.................................................. 83 20 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Dispersion Model Potential. kT/E = 0.85, K = 0.25................................................ 84 21 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Dispersion Model Potential. kT/e = 0.85, K = 0.25................................................ 85 22 Difference in Normal and Tangential Components of Stress Tensor for LennardJones Fluid........................... 87 23 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Anisotropic Overlap Model Potential. kT/E = 0.85, 6 = 0.10................................... 89 24 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Anisotropic Overlap Model Potential. kT/E = 0.85, 6 = 0.10................................... 90 25 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Anisotropic Overlap Model Potential. kT/c = 0.85, 6 = 0.10.................................. 91 Figure Page 26 Interfacial DensityOrientation Profiles for a Fluid of Axially Symmetric Molecules Interacting with Lennard Jones plus Anisotropic Overlap Model Potential. kT/e = 0.85, 6 = 0.10.................................. 92 27 Simplified Schematic Flow Diagram of Fortran Monte Carlo Program Developed for NOVA 2...................... 105 28 Comparison of CDC and NOVA 2 Monte Carlo Results for the CenterCenter Pair Correlation Function for a LennardJones plus Quadrupole Fluid....................... 110 29 Comparison of CDC and NOVA 2 Monte Carlo Results for the Angular Pair Correlation Function for a Lennard Jones plus Quadrupole Fluid for Molecular Pairs in the Tee Orientation..................................... 111 30 Methods of Specifying the Orientation of an Axially Symmetric Molecule........................................ 119 31 Orientation Angles for Axially Symmetric Molecules in an Arbitrary Space Fixed Frame........................... 121 32 Orientation Angles for Axially Symmetric Molecules in the Intermolecular Frame.................................. 122 33 Geometry of a Pair of Diatomic Molecules.................. 142 34 Pair Potential for LennardJones plus Dipole Model Fluid at Primary Pair Orientations....................... 161 35 Pair Potential for LennardJones plus Quadrupole Model Fluid at Primary Pair Orientations....................... 165 36 Surface of the LennardJones plus Quadrupole Pair Potential for the Tee Orientation as a Function of the Quadrupole Strength .................................. 166 37 Pair Potential for LennardJones plus Dipole, Dipole Quadrupole, and Quadrupole Model Fluid at Primary Pair Orientations. p/(co3)1/2 = 1.0, Q/(Eo5)1/2 = 1.75....... 167 38 Pair Potential for LennardJones plus Dipole, Dipole Quadrupole, and Quadrupole Model Fluid at Primary Pair Orientations. p/(Eo3) /2 = 1.75, Q/(co5)1/2 = 1.0....... 168 39 Pair Potential for LennardJones plus Anisotropic Overlap Model Fluid at Primary Pair Orientations ........ 170 xvii Figure Page 40 MeanSquared Displacement of Molecular Centers of Mass for LennardJones plus Quadrupole Fluid............... 173 41 Fluctuation in Temperature for LennardJones plus Quadrupole Fluid........................................... 175 42 Fluctuation in the Ratio of Translational to Rotational Kinetic Energy for LennardJones plus Quadrupole Fluid..... 176 43 Effect of Quadrupole Moment on the CenterCenter Pair Correlation Function at pa3 = 0.85.......................... 194 44 Spherical Harmonic Coefficients g2Z for LennardJones 3 2m2 plus Quadrupole Fluid at po = 0.85, kT/E = 1.294, Q/(Ea 5)l/2 = 1.0.......................................... 196 45 Spherical Harmonic Coefficients g4m for the Fluid of Figure 44.......................... .. .................... 197 46 Spherical Harmonic Coefficients gm for the Fluid of 44m Figure 44.......................................... ....... 198 47 Spherical Harmonic Coefficients g6Z 0 for the Fluid of Figure 44 ........................... ? ................... 199 48 Effect of Anisotropic Overlap Parameter on the Center Center Pair Correlation Function at po3 = 0.85.............. 203 49 Spherical Harmonic Coefficients g m for LennardJones 22m plus Anisotropic Overlap Fluid at po3 = 0.85, kT/E = 1.287, 6 = 0.30................................................... 204 50 Spherical Harmonic Coefficients g4m for the Fluid of Figure 49........................... ....................... 205 51 Spherical Harmonic Coefficients g44m for the Fluid of Figure 49.................................................. 206 52 Spherical Harmonic Coefficients g6, 0 for the Fluid of Figure 49.......................... ?....................... 207 53 Integrands [g220(r) 2g221(r) + 2g222(r)]r2 and [g220(r) + 4/3g221(r) + 1/3g222(r)]r3 for G2 and Ua' respectively, for LennardJones plus Quadrupole Fluid...... 214 xviii Figure Page 54 Angular Pair Correlation Function for the LennardJones plus Quadrupole Fluid of Figure 44 for the Tee Orienta tion (01 = 90, 82 = 0, ( undefined)...................... 217 55 Angular Pair Correlation Function for the LennardJones plus Quadrupole Fluid of Figure 44 for the Cross and Parallel Orientations (01 = 02 = 900)..................... 218 56 Angular Pair Correlation Function for the LennardJones plus Quadrupole Fluid of Figure 44 for a Skewed Orienta tion (1 = 82 = = 45 )........... .. ........................ 220 57 Angular Pair Correlation Function for the LennardJones plus Anisotropic Overlap Fluid of Figure 49 for the Tee Orientation (81 = 900, 82 = 0, < undefined)............... 222 58 Angular Pair Correlation Function for the LennardJones plus Anisotropic Overlap Fluid of Figure 49 for the Endon Orientation (81 = 82 = 1 = 0)............................. 223 59 Surface of the Angular Pair Correlation Function for LennardJones plus Quadrupole Fluid for 8 = 90, 0 = 0, with Q* = 0.5, T* = 1.277, p* = .85 ...................... 226 60 Surface of the Angular Pair Correlation Function for LennardJones plus Quadrupole Fluid for 0 = 90, 0 = 0, with Q* = 1.0, T* = 1.294, p* = .85....................... 227 61 Surface of the Angular Pair Correlation Function for LennardJones.plus Quadrupole Fluid for 8 = 90, = 90, with Q = 0.5, T = 1.277, p* = .85....................... 228 62 Surface of the Angular Pair Correlation Function for LennardJones plus Quadrupole Fluid for 0 = 90, 0 = 90, with Q = 1.0, T* = 1.294, p* = .85...................... 229 63 Surface of the Angular Pair Correlation Function for LennardJones plus Quadrupole Fluid for 1 = 45, 0 = 45, with Q = 0.5, T* = 1.277, p* = .85...................... 231 64 Surface of the Angular Pair Correlation Function for LennardJones plus Quadrupole Fluid for 1 = 45, 92 = 45, with Q = 1.0, T* = 1.294, p* = .85....................... 232 65 Comparison of Peak Heights in the Angular Pair Correla tion Function with Well Depths in the Pair Potential for the LennardJones plus Quadrupole Fluid with Q/(Ca5)1/2 = 1.0 ......................................... 234 xix Figure Page 66 Surface of the Angular Pair Correlation Function for LennardJones plus Anisotropic Overlap Fluid for 6 = 90, 0 = 0, with 6 = 0.10, T* = 1.291, p* = 0.85...... 235 67 Surface of the Angular Pair Correlation Function for LennardJones plus Anisotropic Overlap Fluid for e1 = 90, ( = 0, with 6 = 0.30, T* = 1.287, p* = 0.85...... 236 68 Integrand for Internal Energy, r2u(12)g(12), for LennardJones plus Quadrupole Fluids for the Tee Orientation............................................... 239 69 Integrand for Pressure, r3 du(12 g(12), for Lennard dr Jones plus Quadrupole Fluids for the Tee Orientation...... 240 2 70 Integrand for Internal Energy, r u(12)g(12), for Lennard Jones plus Anisotropic Overlap Fluids for Parallel Orientation............................................... 241 3 du(12) 71 Integrand for Pressure, r d g(12), for Lennard dr Jones plus Anisotropic Overlap Fluids for Parallel Orientation............................................... 242 72 SiteSite Pair Correlation Function for LennardJones plus Quadrupole Fluids................................... 245 73 Possible Square Packing of LennardJones plus Quadru pole Molecules for Interpreting g (r) .................... 247 74 SiteSite Pair Correlation Function for LennardJones plus Anisotropic Overlap Fluids............................ 249 75 Box Representing the Molecular Dynamics System with the Volume Element Sampled for the Filmed Animation Indicated. 253 76 Initial FCC Lattice Configuration of LennardJones Molecules in the Volume Element Sampled in the Filmed Animation.................................................. 256 77 Frame from the Filmed Animation of LennardJones Mole cules Corresponding to the Sixth TimeStep in the Molecular Dynamics Calculation............................ 257 78 Frame from the Filmed Animation of LennardJones Mole cules Corresponding to the 101st TimeStep in the Molecular Dynamics Calculation............................ 258 xx Figure Page Bl Rotations Defining the Euler Angles {(OX) ................ 271 B2 Rotations in the Triangle 123 in the Fowler Model Interface to Define Values for z ...................... 274 max KEY TO SYMBOLS Roman Upper Case A = Helmholtz free energy Ac = Configurational Helmholtz free energy th A. = The i term in the perturbation expansion for i Helmholtz free energy A(zl,zl2) = Function defined in Equation (426) B(zl,Zl2) = Function defined in Equation (427) R CV = Residual contribution to constant volume heat capacity 2 ;mlm2m) = ClebschGordan coefficient D = Diffusion coefficient D = Representation coefficient mn D = Complex conjugate of representation coefficient mn F = Force on molecule 1 1 FI = Function defined in Equation (232) F2A,F2B = Functions defined in Equation (262) and (263), respectively F3A,F3B = Functions defined in Equations (277) and (278), respectively GL = Angular correlation parameter H = Integral defined in Equation (334) xxii I = Moment of inertia I(r) = Functions defined in Equations (665) and (669) I = Integral defined in Equations (240) and (B7) z Inn = Integral defined in Equation (Gl) J = Integral defined in Equation (316) 11 2 J = Integral defined in Equation (689) n K = Coefficient of ellipticity for plane polarized light K = Integral defined in Equation (325) K = Integral defined in Equation (318) KL,KV = Functions defined in Equations (430) and (431), respectively L = Angular momentum operator L = Integral defined in Equation (324) L = Integral defined in Equation (317) N = Number of molecules NA = Avogadro's number P = Pressure th P. = The i component of the local polarization vector 1 Pi(x) = Legendre polynomial of order i P i = Probability of the ith state occurring, defined in Equation (58) Q = Quadrupole moment Q = Reduced quadrupole moment = Q/( )1/2 QZ = General multiple moment T = Temperature T = Critical temperature c xxiii TR T U s V V c Ym Y^m Z = Reduced temperature, T/T c = Reduced temperature, kT/E = Superficial excess internal energy = Volume = Critical volume = Spherical harmonic = Complex conjugate of Yzm = Configurational integral Roman Lower Case a = Semiempirical parameters defined in Equations (350) to (353) c = Constant defined in Equations (421) and (422) c = Cosine of angle .ij c. = Cosine of angle 8. 1 1 c(y) = Cosine of angle y c(zlz2r12) = Interfacial direct correlation function d = Hard sphere diameter defined in Equation (339) d,dd = Maximum allowable step lengths for translational and rotational motion, respectively, of molecules in one Monte Carlo generated configuration f(z rl2) = Interfacial pair distribution function for spherical molecules f(zlI12rl3) = Interfacial triplet distribution function for spherical g(rl2) molecules = Radial distribution function xxiv g (r12) ga8(rl2) g(zlz2rl2) g(r12 W12) g(12) Z l 2m(rl2) fi k 1 k m n n n n s PN PT r12 r12 1 rl2 r m s. 5i 1 t u(rl2) = Centercenter pair correlation function = Sitesite pair correlation function = Interfacial pair correlation function = Angular pair correlation function = g(r12W12) = Coefficients in the expansion of g(12) in spherical harmonics of the molecular orientations = Unit vector aligned along the axis of linear molecule i = Boltzmann's constant = Molecular bond length between atoms = Molecular mass = Index of refraction = Constant, values given in Equations (421) and (422) = Exponent of r in repulsive part of Mie potential, Equation (332) = Power of r in various model potentials = Normal component of the stress tensor = Tangential component of the stress tensor = Vector location of center of mass of molecule i = Vector separation of centers of mass of molecules 1 and 2 = Magnitude of r12 = Reduced distance, rl2/o = Value of r where pair potential is a minimum = Sine of angle 0. 1 = Time = Spherically symmetric intermolecular pair potential xxV u(r12w1 2) u(12) Uo (r12) u (12) u (r12) u s v x12 Y12 y(rl2) z1 z12 z max z1 = Orientation dependent intermolecular pair potential = u(rl2wIw2) = Isotropic, reference contribution to u(12) = Anisotropic contribution to u(12) = Reduced potential, u(rl2)/E = Superficial excess internal energy per unit of surface area = Molecular translational velocity = xcomponent of r2 12 = ycomponent of r2 12 = Function defined in Equation (338) = zcoordinate location of molecule 1 = zcomponent of r12 = Maximum value of zl, defined in Equations (242) and (B6) = Reduced coordinate, zl/o Roman Script = Interfacial area = Total intermolecular potential = Coefficients in the spherical harmonic expansion of the anisotropic pair potential Greek Upper Case = Surface adsorption of component i = Triplet of indices, 2 k for example = Volume in angle space xxvi = Angular velocity th i = The I component of 2 Greek Lower Case a = Adjustable parameter in Equations (426) and (427) ai = Interior angle in the triangle formed by rl2, r13, r23, at molecule i ai,8i = Constants in the predictorcorrector algorithm; values given in Equations (645) and (646), respectively ai,i = Azimuthal and polar angles, respectively, for molecular orientation in the space fixed frame S= /kT y = Angle between the axes of a pair of linear molecules Y = Surface tension F S= Fowler model surface tension th Yi = The i term in the perturbation expansion for surface tension 2 y = Surface tension reduced by potential parameters, yo /E YR = Surface tension reduced by critical constants, 2/3 l Y(Vc/NA) (k ) 6 = Dimensionless anisotropic overlap parameter 6.. = Kronecker delta c = Intermolecular potential energy parameter 0(x) = Unit step function i. = Polar aigle for molecular orientation in the inter molecular frame xxvii 8 = Angular displacement vector 812 = Polar angle for rl2 in spherical coordinates K = Dimensionless anisotropic polarizability A = Wavelength of incident light in light scattering experiments S= Perturbation parameter p = Dipole moment = Reduced dipole moment, 1/(c3 )1/2 S= Random numbers p = Fluid number density p(zl) = Interfacial number density profile p(z 1w) = Interfacial number densityorientation profile pl(zl l) = First order term in perturbation expansion for p(z1w1) 0 = Intermolecular potential distance parameter ai = Standard deviation of property i T1 = Torque on molecule 1 i = Azimuthal angle for molecular orientation in the inter molecular frame ^12 = Azimuthal angle for r2 in spherical coordinates = Function defined in Equation (319) i = Set of variables specifying the orientation of molecule i xxviii Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SURFACE TENSION AND COMPUTER SIMULATION OF POLYATOMIC FLUIDS By James Mitchell Haile December, 1976 Chairman: Keith E. Gubbins Major Department: Chemical Engineering A third order thermodynamic perturbation theory for the surface tension of polyatomic liquids has been developed using a Pople reference fluid. The expansion includes terms containing the unknown pair and triplet interfacial correlation functions, go(zrl12) and go(z l1223 ) for the reference fluid. These are removed by using the Fowler approxima tion for the interface. The resulting theory has been tested against Monte Carlo results for the Fowler model surface tension of a fluid whose molecules interact with a LennardJones plus dipole potential. The behavior of the theory compared with similation parallels that of bulk fluid properties; namely, the second order theory agrees with the Monte Carlo results for small values of the dipole moment up to P/(o3)1/2 = 0.6. For larger dipole strengths neither the second nor third order theories agree with the computer simulation results. How ever, when the third order expansion is recast in the form of a simple [1,2] Pade approximant, the theory agrees with the Monte Carlo results up to dipole moments as large as w/(Ca3)1/2 = 1.75. This Pade theory has been used to calculate the surface tension of pure dipolar and xxix quadrupolar liquids. The theory agrees with experimental values of surface tension in the neighborhood of the triple point; however, the Fowler model does not give the correct temperature dependence of the surface tension. The theory has also been used as a basis for develop ing useful correlations of surface tensions of pure polyatomic liquids. A first order perturbation theory has been developed for the interfacial densityorientation profile p(z 1W) for polyatomic fluids. Upon introduction of a Pople reference, the first order term p1(z 11) vanishes for multipolar anisotropies, but does not vanish for anisotropic overlap or dispersion potential models. Calculations of p(z 1W) for axially symmetric molecules interacting with each of the latter potentials have been performed, using a LennardJones reference fluid and the inter facial pair correlation function model used by Toxvaerd. The calculations indicate that the axially symmetric molecules have preferred orientations in the interfacial region. A method has been devised for using a NOVA 2 minicomputer with 32K words of core and external disc storage to perform Monte Carlo simulations of 128 nonspherical molecules. The simulations generally require several days of continuous calculation; however, several equilibrium property values and values for the angular pair correlation function g(rl2 W12) at five to seven specific orientations may be obtained at a fraction of the cost of doing the calculations on a full size machine. Results from the minicomputer compare within the statistical precision with results previously obtained on CDC and IBM machines. The method of molecular dynamics has been used to study systems of 256 molecules interacting with LennardJones plus quadrupole and xxx LennardJones plus anisotropic overlap potentials. The equilibrium properties determined include configurational internal energy, pressure, Fowler model surface tension, Fowler model surface excess internal energy, mean squared force, and mean squared torque. In addition, the coefficients g 2 m(r12) in an expansion for the angular pair correlation function g(r12W1W2) in terms of products of spherical harmonics of the molecular orientations are determined. Sitesite pair correlation functions g a(r) are also found. Relations are developed between the g l2m(r12) coefficients and the above listed equilibrium properties for several anisotropic potential models. Study is made of orientational structure in quadrupolar and overlap fluids via g(r 12W12) as obtained from the recombined spherical harmonic expansion. A method of producing filmed animations of molecular motions from molecular dynamics data is described. xxxi CHAPTER 1 INTRODUCTION The development of rigorous methods for prediction of properties of liquids must come from an understanding of how molecules are dis tributed and interact with one another. Such fundamental understanding of liquids may be pursued from three directions: theory, experiment, and computer simulation. Molecular theory is the domain of statistical mechanics; molecular level experimentation usually involves light, xray, or neutron scattering, while computer simulation embodies the Monte Carlo and molecular dynamics techniques. Computer simulation adds a powerful new dimension to the study of matter which in no way supplants the other two methods of study. As indicated in Figure 1, simulation complements both theory and experiment. Thus, comparison of simulation with laboratory experiment provides information about the molecular interactions in the liquid. On the other hand, com parison of simulation with theory provides a stringent test of the theory. In spite of current gains being made in the study of liquids, fundamental understanding of fluidfluid interfacial phenomena has been slow to develop. The difficulty with study of interfacial properties arises from the inhomogeneity, nonuniformity and anisotropy of the interfacial region between homogeneous fluid phases. In general, microscopic fluid properties vary through the interface from their values in one bulk phase to their values in the other bulk phase, as THEORY Potential Theory + Potential EXPERIMENT Theory Figure 1. Relation of Theory, Experiment, and Computer Simulation in the Study of Liquids COMPUTER SIMULATION in Figure 2. These properties include density, refractive index, dielectric constant, etc. The oft measured macroscopicc) properties of interfaces are, in many cases, related to integrals over the microscopic properties [1]: X = k [MB M(z)]dz (11) 00 where X represents a macroscopic property, M is a microscopic property, subscript B indicates, usually, a bulk phase value and k is a pro portionality constant. Examples of specific properties having the general form of (11) are given in Table 1. Equation (11) indicates that prediction of macroscopic interfacial properties and comparison with experiment is not a completely satisfactory test of theory. Of more value is understanding of the microscopic properties and their relations to macroscopic properties of interest. In the case of theoretical study, such an approach must lie in statistical mechanics. The prediction of macroscopic property values in a variety of situations, e.g., over a large portion of the phase diagram and in multicomponent systems, remains a vital engineering concern. In the case of interfacial problems, the macroscopic property of interest is the interfacial tension. The ability to predict interfacial tensions is required in the equipment and process design and operation of numerous fluidfluid contacting operations, such as distillation, solvent extrac tion, liquid membrane separation techniques, and tertiary oil recovery. Interfacial tension is also important in understanding various biological processes such as blood oxygenation and eye lubrication. Of particular P Figure 2. Variation of Fluid Density with Position through a Planar VaporLiquid Interface TABLE 1 Examples of Macroscopic and Microscopic Interfacial Properties Related by Equations of the Form (11) 00 X = k [MB M(z)]dz 0 CO S(z) L ]dz + [i(z) p ]dz (12) f 0 000 Y = j [PN pT(z)]dz (13) 0T 7T 2 P (z) P (z) z x K = (n + dz (14) A P (z) P (X) Reference [2] References [3,4]" F. = surface adsorption of component i pi = density of component i V = bulk vapor phase density i L = bulk liquid phase density Y = surface tension PN = normal component of stress tensor pT(z) = tangential component of stress tensor K = coefficient of ellipticity for light plane polarized at 450 to the plane of incidence, when incident at Brewster's angle th P. = i component of local polarization vector n = bulk phase index of refraction A = wavelength of incident beam importance is the determination of how molecular characteristics of the fluid contribute to the interfacial tension. An important goal in design considerations is the improved efficiency of fluidfluid contacting operations by modifying the system (e.g., by introduction of additives) in order to lower the interfacial tension. This, in turn, leads to consideration of how molecules absorb and orient themselves in the interfacial region. 1.1 Theory of Surface Properties Readily used methods currently available for predicting inter facial tension have been reviewed [5,6]. These methods include corresponding states approaches, techniques based on regular solution theory and scaled particle theory, and more empirical methods. These methods are generally applicable to spherical and nearly spherical molecules, nonpolar polyatomics and their mixtures. These methods fail to accurately predict interfacial tensions for strongly polar, quadrupolar or associating liquids. Almost all rigorous theoretical work on interfacial phenomena has been for fluids with spherically symmetric molecular interactions. Two rigorous relations have been developed for surface tension. The KirkwoodBuff equation relates the surface tension y to the inter facial density profile p(zl) and the interfacial pair correlation function g(zlz2r12) [2,7] S2 2 Y = 2 dz d12 dr P(Z l)p(z2) g(zlz2rl2 r2 x12 12 (15) 00 The KirkwoodBuff equation is valid for spherical molecules and assumes the intermolecular potential to be a sum of pair potentials. This relation is intractable as it stands because of the unknown function g(zlz2r12). Consequently, the KirkwoodBuff relation has been studied using various simplifying models for the interfacial region. The simplest such model is due to Fowler [8] and assumes an abrupt transition from liquid to vapor phase. Further, the vapor phase density is assumed to be negligible compared to the liquid density. Thus, the model may be expressed as: p(zl)p(z2) g(zlz2rl2) = 6(zl)e(z2) Pg (rl2) (16) where 6 is the unit step function, 1 if x > 0 6(x) = (17) 0 if x < 0 subscript L indicates a bulk liquid property and the negative z direction is into the liquid. The resulting FowlerKirkwoodBuff expression for surface tension is: 2 m P u 4 du(r2) Y =8 dr12 12 dr12 gL(r2) (18) 0 As could be expected, the FowlerKirkwoodBuff theory works well near the fluid triple point but gives increasing errors in surface tension as the temperature is raised towards the critical point. Recent evidence indicates that the good agreement at the triple point is due to cancellation of errors [6]. The second rigorous relation for surface tension is a generalized van der Waals equation which gives the surface tension in terms of the density gradient dp(zl)/dzl and the interfacial direct correlation function c(zlz2r12) [9]: kTf f d f dp(z ) dp(z2 2 2 Y = dzI dz2 dxl2 dyl2 dz dz c(zlZ2r2)(xl2 y12) (19) 00 00 CO This relation is more general than the KirkwoodBuff expression since no assumption of pairwise additive potential is made in its derivation. It has the further advantage that the direct correlation function c(rl2) is generally of shorter range than the pair correlation function g(r12). The generalized van der Waals equation (19) has not been as thoroughly studied as the older KirkwoodBuff formula. The KirkwoodBuff equation has recently been generalized to nonspherical molecules [10]. From both a practical standpoint and a desire to gain understanding of interfacial phenomena, there is strong need for using this new relation as a basis for developing predictive methods for surface tension of polar and quadrupolar systems. Evaluation of the interfacial density profile p(zl) is required in order to determine the surface tension from the rigorous expressions (15) and (19). Study of the interfacial density profile is also of interest per se. Nearly all studies thus far have been for spherical molecules. A variety of approaches have been used: a) van der Waals theory [11], b) constant chemical potential through the interface [12,13], c) first BornGreenYvon equation [14], d) constant normal pressure through the interface [15,16], e) minimiza tion of system free energy obtained by perturbation theory [17,18], and f) computer simulation [19,20,21]. Determination of surface tension for nonspherical molecules will require knowledge of the interfacial densityorientation profile p(zl 1), where wl is a set of Euler angles specifying the orientation of molecule 1. Further, there is considerable interest in determining how modification of molecular orientation affects "the surface tension. 1.2 Computer Simulation Methods In the computer simulation approach to the study of liquids either the Monte Carlo or molecular dynamics technique may be used. These methods provide detailed information on equilibrium and time dependent properties and on liquid structure for molecules interacting with known force laws. In addition to use of simulation results as a standard for evaluating theories of liquids, there is now serious interest in using simulation as a vehicle for evolving and evaluating realistic potential models for liquids. Much of the simulation work has been performed for simple, spherical molecules. Only in the past few years has simulation of nonspherical molecules been undertaken. There is need of exploiting these simulation methods as fully as possible since they provide a wealth of detailed information for study. Much of the effort in the history of computer simulation has been expended in developing the methods, demonstrating their usefulness and potential applicability and, in general, making simulation a respectable research tool. There is currently a need to explore methods for improving the efficiency of these calculations with a view towards conserving computer resources and making the simulation techniques accessible to more users. One possible approach is the use of a minicomputer for performing simulations of liquids. 1.3 Outline of Dissertation This work is divided into two parts. Part I includes Chapters 24 and is concerned with the theory of vaporliquid interfaces for nonspherical molecules. In Chapter 2 a statistical mechanical perturbation theory is developed for the surface tension of polyatomic fluids. The general first order perturbation term is obtained and the second and third order terms are found when a Pople reference is used. Specific relations for the perturbation terms are given for several anisotropic potential models: dipole, quadrupole, overlap and dispersion. The corresponding perturbation terms are also derived when the Fowler'approximation for the interface is introduced. In Chapter 3 numerical calculations are reported for surface tension using the perturbation theory in the Fowler model. The results are compared with experiment and computer simulation. Semiempirical methods based on perturbation theory are explored for correlating surface tensions of polyatomic and polar substances. Chapter 4 develops a perturbation theory for determining the vaporliquid interfacial densityorientation profile of fluids with anisotropic potentials. Calculations are presented for overlap and dispersion model interactions. The calculations predict that these nonspherical molecules exhibit preferred orientations in the inter facial region. Part II of this work describes computer simulation studies of linear molecules. Chapter 5 reports development of a method for performing Monte Carlo calculations for linear molecules on a NOVA 2 minicomputer. Chapter 6 describes the molecular dynamics method for linear molecules, including: derivation of expressions for efficient evalua tion of the force and torque exerted on a molecule due to various anisotropic potential interactions, the method used to solve Newton's equations of motion and the molecular dynamics algorithm for linear molecules, evaluation of the coefficients in a spherical harmonic expansion of the angular pair correlation function, and development of relations between these coefficients and various fluid equilibrium properties. Chapter 7 reports the results of the molecular dynamics calculations for equilibrium property values obtained for Lennard Jones plus quadrupole and LennardJones plus anisotropic overlap fluids. Comparisons of the equilibrium property values are made with perturbation theory predictions in the case of the quadrupole fluid calculations. Study is made of the local structure in these fluids using the molecular dynamics determined angular pair cor relation functions. A method for producing filmed animations of molecular motions from molecular dynamics calculations is presented. Chapter 8 draws conclusions from this study. PART I THEORETICAL STUDY OF FLUID INTERFACES CHAPTER 2 THEORY OF SURFACE TENSION 2.1 General Expressions for Surface Tension of Polyatomic Fluids There are two rigorous expressions for the surface tension of polyatomic fluids. One is the generalized van der Waals equation (19). The second is a generalization of the KirkwoodBuff expression (15) which has been previously derived [10,22]. In this section the deriva tion for the general KirkwoodBuff equation is summarized and its simplified form obtained when the Fowler approximation is made. 2.1.1 Generalized KirkwoodBuff Formula Consider a planar interface between vapor and liquid phases with a coordinate system oriented such that the xy plane is in the plane of the interface and the +z direction points into the vapor. The two phase system has N molecules.in volume V at temperature T and interfacial area S. Thermodynamically, the surface tension y is related to the Helmholtz free energy A of the two phase system by [23]: Y= (21) Y NVT Since only the configurational part of the free energy Ac depends on the interfacial area, (21) may be written as: y = kT nn Z T (22) NVT where Ac = kT nn Z (23) is the statistical mechanical definition of the configurational Helmholtz free energy, k = Boltzmann's constant, and Z is the con figurational integral. For nonspherical molecules: N N U(r Nw ) (24) Z= J*** dr dw e U W where g = 1/kT and U is the full intermolecular potential. For non spherical molecules U depends on the orientations w of the molecules, as well as the positions of their centers of mass r. The orientations w are usually specified by a set of Euler angles w = O6X between a bodyfixed reference frame on the molecule and a spacefixed frame located external to the system. Substituting (24) into (22): kT f rr N N U(rNWN)) Y = Z dr dw e JNVT (25) The differentiation in (25) cannot be done immediately since the N integration limits on the integrals over r depend on S. We, therefore, follow Green [24] and change r variables using the transformation: 1/2 , x = S x y = S12 y' (26) V , z =z Performing the differentiation in (25) and changing back to the old variables gives: N N N NN Y 1 f f e8(r NN) S= dr dr eN Z J J LauN NN aU(r w ) as INVT Assuming the potential to be pairwise additive, U = Z u(r.i W .) i (27) becomes wd f(((Du(12)( S i f2...f drl dr2 ld 2 f l m2 aS NVT YNVT where u(12) E u(r12 1W2) and the definition of the angular pair dis tribution function has been used: f/ ( N(N1) f( l21 l2) = Z f N N U(r w ) dr 3* dr dw3. dwN e  3 N 3 N Integrating (29) over xl, Yl and transforming r2 to r12 gives: Y = 2 u dzl dr12 dwldw2 () Su(i12) N Y = 2 l 2 2 f(zl 2 2) S S Su(12) S can be evaluated for nonspherical molecules by considering: as au(12) 12 au(12) BS 9S Brr 12 1 Su(12) 32 u(12)1 2 12 r z12 az 212z12) Hence, (211) becomes: (27) (28) (29) (210) (211) (212) (213) Y =I dzl drl2 dwldw2 f(zlrl2l2) rl au(12) r12 312 Du(12) 12 8Z 12 (214) Defining the angle average by: 1Wm2 (''')dw dw2 Equation (214) can be written: I dr12 au(12) Dr12 au(12)  3zl2 (z l2 ]> z12 z 12J 1)2 (217) Equation (217) is one form of the general KirkwoodBuff equation. Another form may be obtained by transforming (217) to spherical coordinates. Sau(12) z 12'Y12 Since [25] = cos 612 (218) then _3z2 u(12)  az u12) 12 12 au(12) = 2P(cos )r2 12 2 12 12 Dr12 au(12) + 3 sin 012 cos 612 u12) S12 a 12 1 2 where P2(x) is the second order Legendre polynomial, P2(x) = (3x 1). 2 2)2 Substituting (219) into (217) gives: 12  2 j where S= f dw (215) (216) co f2 f , S= 4 jdz1 0o fau(12) I r12 sin 012 . rl2 (au(12)] a "12 (12 8u(12) ar12 (219) Y = YA + YB (220) YA 2 dz drl2 P2(cos 12) r2 2 ur12 12 o00 Y = 4 dz I dr12 sin 812 cos 012 B 4 21 1212 H2 I 2 00 Equations (220,21,22) give the general KirkwoodBuff equation. The equation applies to general shaped molecules, the only assumption being a pairwise additive potential. In the case of spherical molecules, the potential u(12) goes to u(r12) so that the derivative in the yB term vanishes and yA reduces to the KirkwoodBuff equation (15). 2.1.2 FowlerKirkwoodBuff Equation for Nonspherical Molecules The general expression for the surface tension (220) is in tractable as it stands due to the unknown distribution function f(zl r2wl2). Various simplifying assumptions can be made for f to enable calculations to be performed. The simplest assumption is that due to Fowler [8]: 2 pL f(z1 12l2) = 0(z1) e(z2) 2 gL(12) (223) where 0 is the unit step function as in (17), pL is the bulk liquid density, and gL is the bulk liquid angular pair correlation function. Introducing the Fowler model (223) into (221) allows the integration over z and w12 to be performed giving: 2 00 F L T 4 u(12) YA = 8 dr2 r12 A 8 112 1lm2 0 where the superscript F indicates Fowler model. Similarly, (222) reduces to: Y2 dr r F 3 2 2 3 <() u(12) 2 B = 32 PL dr12 r12 e W F Further, the y term can be shown to vanish by symmetry arguments. One B such argument is that 3u(12)/6 12 is proportional to the 812 component of the force on molecule 2 due to molecule 1. This should vanish by symmetry when averaged over l1 and w2. Thus, the FowlerKirkwoodBuff equation for nonspherical molecules is: 2 20 F L 4 au(12) y f dr12 12 8 1212 L r12 1Wl2 0 For the special case of spherical molecules, (226) reduces to the usual FKB expression (18). 2.2 General First Order Perturbation Theory for Surface Tension 2.2.1 Rigorous First Order Term The difficulty with using the general KirkwoodBuff expression (220) lies with the, in general, unknown interfacial distribution function f(zlr2wlw2). This problem may be avoided by use of per turbation theory, if the interfacial pair distribution function is available (or may be approximated) for some reference substance. Perturbation theory for fluid properties may be developed by con sidering the anisotropic pair potential u(12) to be the sum of a (known) reference term u and a perturbing term ul: u(12;A) = uo(12) + Au (12) (227) where X is a perturbation parameter such that when A = 0, (227) gives the reference potential and when A = 1, (227) gives the full potential. Expanding the twophase system Helmholtz free energy (23) in powers of A and setting X = 1 gives: A = A + A1 + *" (228) where A1 = drldr2 1 f1O W 1 W2 and f is the interfacial angular pair distribution function for the 0 reference system. The corresponding expansion for surface tension is obtained by applying (21) to (228): Y = Y + Y1 + *** (230) where Y= I ( dr dr F(r2) (231) as 1 2 1 112) INVT and F1(zll2) E ( <1f ( 2 u(12;A) > (232) l A=0 wlw2 Applying Green's method as in section (2.1), integrating over xl and Y,' and transforming rl to r12, (231) becomes: 20 1 f f 12 flzs I = 2 dl r S[f 2NVT S 3Fl(zrl2)/@S may be evaluated by considering: 12 ,F (zlrl2) 3 I (zl 1 2 z1 =S + S as r12 a1z as 12 3Fl(z 12) 2 r 12 12 3 z Fl(z l12) 2 12 z12 aFl(Zlr12) Z F (zIl  zz Thus, (233) becomes: CO F1 f(z 1r 2) l2 3 Fl(zil12) Y 2 dZ 2 1 zz 2 "12 az12 00 _12 Fl(z rl2) 2 Dr 12 The second two terms in (236) can be integrated by parts and shown to cancel, leaving: 1 fA f2fu(l2;X) Yl = 2 12 dzl A O 00 af z (z) D 21 1 2 Equation (237) is valid for spherical or nonspherical reference fluids; the only assumption has been a pairwise additive potential. Note that if the reference fluid is taken to be a Pople reference, defined by (248) below, then Y1 vanishes. (233) (234) (235) (236) (237) SFl(Z142) S 2.2.2 First Order Term in the Fowler Model To make (237) amenable to calculation, the Fowler approximation (223) may be introduced: 2 PL 2Q2 Sau(12;X)( dr2 dzl Z 1, 6(zl) e(z2)< fu(12; 0g (12)> 12 21 1 z L 8)2 (238) Substituting (215) and changing the order of integration: 2 L 2Q4 dr dw dw 12 1 2 fau(12;A)J ( ax X=O (239) where Iz I dzl Z1 (240) The integral I may be evaluated by parts to give: z I = go,(12) zmax z oL max z = max r 12 cos e12 (241) (242) if 6 >2 if 12 if e < 7T 12 2 and 012 is the spherical coordinate polar angle Figure 3. Equation (239) becomes, therefore: for r12' as shown in Sdwdw2 u(12;A) 1 2 ( @A ]= goL(12) zmax (243) where 2 PL 2Q4 S6(zl) 6(z2) goL(12) 12 Two Possible Values for the Maximum z Value (zmx) for a Pair of Molecules in the Fowler Model Interface Figure 3. Reintroducing the angle average (215) and noting that g(12) in the bulk liquid is independent of 12, (243) can be written as: 2 o F _L 2 S 2 dr2 12 1 A X= oL W2 12 max 0 Performing the integration over w12 gives: F L 3 [u(12;X) Y fdr r < (12)> (245) 1 22 12 12 X =0 oL (242 0 Equation (245) is the general result for the first order perturbation term in the Fowler approximation. 2.3 Perturbation Theory for Surface Tension using a Pople Reference When the general expansion for surface tension (230) is in creased to higher order, the second order term is found to include a term containing the reference fourbody distribution function. Higher order terms in the expansion contain even higher order multibody terms. These complicated terms can be made to vanish up to at least the second order term in the expansion by using the isotropic reference potential first suggested by Pople [26]: uo(rl2) = 2 (246) With this choice of reference, (227) becomes: u(12;X) = u (rl2) + Aua(12) (247) (246) and (247) together give the simplification: a3u(12; A) > X I=0 loW2 = 1 2)0 = 0 (248) Then (229) and (237) have: (249) Al = Y1 = 0 and the second and third order terms (A2,A3,Y2,y3) simplify. Thus, in the Pople expansion, (228) and (230) become, to third order: AC = AC + A + A (250) o 2 3 Y = YO + Y2 + Y3 (251) 2.3.1 Derivation of Second Order Terms, Y2A and Y2B The second order term in (250) for the bulk phase liquid is [27]: 2 A= dr dr A2 4 2 a 12 (252) where g1(12) is the first order term in a perturbation theory expansion of the angular pair correlation function. When the expansion is about a Pople reference, g1(12) is given by [281: g (12) = Bu (12) g (r ) dr [ + ] x 1 a o 12 3 a W3 a go(r12r13r23) (253) where g (r12r13r23) is the triplet correlation function for the reference. When the anisotropic potential contains only spherical harmonics of order R # 0, such as multipoles, then the angle averages in the second term in (253) vanish. Such potential models as anisotropic overlap and dispersion contain X = 0 spherical harmonics, in which case the second term must be included. Equation (252) may be written, therefore, as A2 = A2A + A2B (254) where 2A dr dr go(r ) 2A 4 2 (r12) a l> 2 (255) p3 A2B = dr dr2d go(r12r13r23) r~, r2 and r3 are each integrated indices and (256) may, hence, be pN3 dr dd 2 drd2dr3 o(r12r13r23) l2m3 (256) over in (256), the indices are written as: (257) For a fluid nonuniform in the zdirection, (255) and (257) generalize to: A = dr dr p (zl o(z2) o(12 2A 4 =l 2 o 1 o 2 0'zlrE12) Since dummy A2B 2B I dr drdr3 (z )p(z2)P (z3) go(z12 13)< 12)Ua(13 > 23 (259) The corresponding terms in the expansion for surface tension (251) may be found by applying (21) to (258) and (259): Y2A = 4 f 1 drl2 F2A(zr2) NVT 2B = dr dr dr F2B(Z l ) 2B 2 f 1 2 3 2B 1 12 NVT where F2A (zl)P (z2) 8o(Zr2) 2B Po(zPP(z2)Po(z3) go(z12r13) (260) (261) (262) (263) Equations (260) and (261) have the same form as (231); thus they may be evaluated in a similar manner to give, analogous to (237): Y = 2A 4 Y = f 2B 2 dr12 dz z 2A2 L12 1 1z 1 dr12d13 dz 12 13 1j aF2B (Zl112) z 1 zI Substituting (262) and (263) into (264) and (265), respectively, and using the relations: (264) (265) fo(z12) = o(zl)Po(z2) go(zlr2) (266) fo(Zlr12r13) = Po(zl)po(z2)po(z3) go(Zlr12r13) (267) gives, S d 2 rfo(zlr12) Y2A dr dzl z (268) 2A f4 12 'a l2 z 1 zI 2B = dr dr13 J dz1 z f (z 2113) 2B 2 12 13 a a 23 1 (269) The derivation of (268) and (269) only assumes a pairwise additive potential and use of a Pople reference fluid. If the aniso tropic potential contains only k # 0 spherical harmonics, then y2B vanishes because of (253). 2.3.2 Derivation of Third Order Terms, Y3A and Y3B The third order term in (250) for the bulk phase liquid is: A dr dr (270) 3 =6 1 2 a g1 2 where g2 is the second order term in the perturbation theory expansion for g(12) [28]. If only anisotropic potentials containing Z # 0 spherical harmonics are considered, g2(12) simplifies. Using arguments anologous to those for simplifying the A2B term in (256), and generalizing to a fluid which is nonuniform in the zdirection, (270) becomes (for A # 0 harmonics): A3 = A3A + A3B 2 A3A 12 A 2 f 3B 6 dr dr2 o(z )P (z2) go(rl2) i2 d12 o o(2 o(z3) o(Z1223) x ldr 2d Lr3 P.(zl)p.(z2)po(z3) g.(1z1213) x Applying (21) to (271) gives: Y3 = Y3A + Y3B dr dr F3A( 12 NVT NVT g2 Y3B =6 7 dr rd2d3r F3B(zlrl2l3) NVT F3A(zlrl2 Po(zl)Po(2) go(z 12) (275) (276) (277) 1 o3(12)u (13)u (23)> F3B(zll2rl3) Po(zl)Po(z2)Po(z3) go(zl23) (278) Again, (275) and (276) have the same form as (231) and, consequently, they may be evaluated in the same manner to give (271) (272) (273) where (274) B2 Y = 12 3A 12 analogous to (237) (the Jacobian of the transformation dr dr dr = J dr dr dr is unity): 1 12 13 S2 f (z r12) 3A 12 dzl z1 zl (279) 3A 12 12 < ul2 1 1 azl 00 R2 f (z 1l2l 13) Y = dr dr3 dz z 0 3B 6 12 13 a a a w2m3 J 1 z (280) where fo(z112) and fo(zl2E1 13) are given by (266) and (267), respec tively. The derivation of (279) and (280) assumes: a) pairwise additive potential, b) use of a Pople reference fluid, and c) the anisotropic potential contains only terms with spherical harmonics of order R # 0. 2.4 Fowler Model Expressions for Perturbation Terms Y2A 2B_ 3A' Y3B in Pople Expansion The Fowler approximation for the spherically symmetric Pople reference may be expressed as: fo(zll2 = (zl)6(z2) ~ goL(r12) (281) fo(z!121 = (z)(z2) (z3 goL(r2) (282) Putting (281) into (268), (269), and (279) for y2A' Y2B and Y3A' respectively, and using (282) in (280) for Y3B' all give terms analo gous to (239) for y1 in the Fowler model. In each case, the integration over zI may be done by parts giving, analogous to (244): 2  F PL r 2 2 f Y1A dr r2 < (12) g(rL ) d z (283) 2A 4 12 12 a ( 12) goL r j 12 max 0 o,3 F L 2 2 r Y dr r dr r 2B 2 12 12 13 13 u a a () u 3) 0 0 x f dw12 dw13 go(rl2rl3r23) Zmax (284) 2 2 Y dr 2 (r f d 12 285) 3A = 12 2 12 a 2 oL 12 12 max 0 2 3 F L 2 2 Y 1 dr 1 2 dr r32 3B 6 12 12 13 13 a a a W 12m3 0 0 x dWl2 d13 goL(r12rl3r23) Zmax (286) In the case of the twobody terms, Y2A and Y3A' the values for z maxare given by (242); see Figure 3. Hence, the integration over S12 can be performed giving: 2 co 2A = drl2 r2 goL(rl2) 0 F Lp 3 3 YA = dr r2 g oL(rl2) (288) 3A 12 12 12 goL 12 a W1 2(288) 0 The integration over 12 and w13 in the threebody terms, y2B and 3B' cannot be performed immediately since zmax and g (r ,2'r13 r ) 3B max goL 12 13 23 depend on these angles. A portion of this angle dependence may be integrated out, however, if the angles 12012 13 13 are transformed to a set of angles which include Euler angles specifying the orientation of the triangle whose vertices are the molecules 1, 2 and 3. Details of the coordinate transformation and integration over the Euler angles are reserved for Appendix B. The resulting expressions are: 2 3 r12+ r13 F 7 L fd f B 2 fdr2 r2 dr13 r13 J dr23 r23 goL(r12r13r23) 0 0 r12 r131 x (r + + r23) (289) r + r 2a2 2 3 0 m 12 r13 F L f f 3B 6 dr2 r12 dr3 r13 dr23 r23 gL(r12r3r23 0 0 ir12 r131 x l 3 (rl2 + r3 + 23) (290) a1a 12w3 12 F F F F Computationally convenient forms for 2A' 3A 2B and y3B are derived from (287), (288), (289), and (290), respectively, in the F F next chapter. Equations for each of these perturbation terms Y2A' Y2B' F F Y3A' Y3B for specific anisotropic potentials are tabulated in Appendix D. In the case of bulk fluid properties, Stell et al. [29] have obtained improved results over the perturbation theory by resumming the series (250) in the form of a Pade approximant: A Ac = A + (291) o A 3 1 A2 The analogous form for the surface tension expansion (251) is: 32 1 Y2 Calculations based on (292) in the Fowler model are given in the next chapter. 2.5 Superficial Excess Internal Energy from the Pade Perturbation Theory for Surface Tension The surface excess internal energy U is related to the temperature s dependence of the surface tension by the classical GibbsHelmholtz equa tion: Us = Y T l(I (293) NV where u is the surface excess internal energy, U per unit of surface area S: U u = (294) s S The Fowler approximation does not predict the correct temperature depen dence of the surface tension for fluids of spherical molecules. Further, the Fowler model may be used to obtain an equation for u in terms of S the pair distribution function for the bulk liquid [2], and this second equation is inconsistent with (293) [30,31]; Freeman and McDonald show that these two expressions give quite different results for use in the case of LennardJones liquids [31]. To determine the validity of (293) for the Pade perturbation theroy for surface tension presented above, we combine (292) and (293) to obtain: u = u + u (295) s so sa where u is the isotropic reference contribution, and so DY3 2 3 Y T T 12 2 u = (296) sa 2 21 I 2^ CHAPTER 3 NUMERICAL CALCULATIONS OF SURFACE TENSION The Fowler model perturbation theory developed in Chapter 2 has been used to calculate surface tensions for pure polyatomic fluids. The intermolecular potential used for the fluids is of the form: u(rl2 1W2) = uo(rl2) + ua(r l2WW2) (31) where u is the LennardJones potential, Uo(rl2) = 4E[(/rl2)12 (a/r12)6] (32) and u is the dipole, quadrupole, anisotropic overlap or anisotropic dispersion potential. Equations for these potentials are given in Appendix C. In these calculations the superposition approximation is made for go(r12r13r23): go(r12r13r2) = go(r12) g(r23) go(r3) (33) and Verlet's molecular dynamics results are used for go(r) [32]. In Section 3.1 forms amenable to calculation are presented F F F F for y2A' Y2B' Y3A, and y3B In Sections 3.2 and 3.3 calculations of surface tension and surface excess internal energy are presented for various model fluids which obey (31); comparisons with Monte Carlo calculations are made for surface tension of polar liquids. In Section 3.4 results for real fluids are given and compared with experimental measurements. In Section 3.5 the perturbation theory is used as a basis for developing a correlation of surface tension for polyatomic liquids. F F F F 3.1 Evaluation of Y2A F 2B and 3B Rewriting (287,88,89,90) in dimensionless form: *2 p* F* F 02 p L ( *3 *2 Y2A 2A 4 7 dr12 r12 g o(r 1) (34) 0 *2 oo A 3A 12 T*2 12 12 12 0L L *3 *3* 3A C 12 0 *. S*3 co CO 12 r13 F 2 2 T L r '2B f dr r dr r 23 23 2B E *2 12 12 f 13 1 23 23 T A 0 0 r12 r13 * SgoL(r2r3 r23) (r+ r + r 23) S12 13 23 12 13 23 a a w1w2m3 * r + r *3 o o 12+ 13 F 02 2T * =T L dr2 r12 dr13 r1 dr r 3B 6 *2 12 12 13 13 23 23 0 0 Jr12 r 13 g r r13 ( * x goL(r1r3 r2) (r 2+ r13+ r23) (37) 3 * where p = po T = kT/E, u = u /E, and r = r/a. a a F* _ Y2B = Y F* _ Y3B = Y The angle averages in (34,5,6,7) have been evaluated [33] and are tabulated in Appendix A. Substituting those expressions into (34,5,6,7) gives: F* r F* Y* = Y2A(A;ss') (38) A ss F* F* Y2B Y2B(AA';ss') (39) AA' ss yA* = YA(AA'A";ss'S") (310) AA'A" ss's" YB y3B (FAA'A"; ss's") (311) AA'A" ss's" In each of these equations A = 1 2 ; in (310) A' = 1'', A"= 2 while in (39) and (311) A' = A'A'', A" = 2 In these equations, *2 F* PL (2 + 1) * Y (A;ss') = J (p IT 2A (21+1)(22+1) n +n 1 16T s s n1n2 x Es(A;nn2 ) E,(A;nn2)/E2 (312) *3 F* PL y YB (AA';ss') = 6 6 6 6 6 L ( ;n n 2B 8 12 1 1 20 3'0 1ss S n1 F*A 73A(AA'A";ss's") = *2 PL 96 1/2T*2 96ff T (2+l)(22'+1)(2Z"+l) * x J (P ,T ) ns+n ,+ns ,,1 ' S S' s x 0 x 1 n2n1 2 n n '2 " n,,n, Ti 0' ORY S1 Pt" 2 2 0 0 [ Zv n all 1 22 a' I 2' _a1j Ln F* s") Y3B (AA'A"; ss's") *3 2 p*3 PL+Z,+. 72 pL 2+c'+ 6 T*2 6 k 11 2 2 [(2k+l)(2k'+1)]l/2 3A" (2a1+1)(2k2+1)(2k'+1) r3 SKY(Z'a";n n'n") kIl sss nl+n2+ n ()1 x Es,,(A";n )/E In these equations, is a 3j symbol, M 9 6. 1j is the Kronecker delta, n n, m' m"i is a 6j symbol [34,35], and 22 2 2 2 (314) 0 0 (315) a' a x Es(A;n n2) Es'(A';n'n2) Es,,(A";nl"n)/E3 x 2n nIn2n3 Es(A;n1n2) E s,(A';nn') 1 1 1 ." R2 22 is a 9j symbol [36]. The E are coefficients in a 2 2 2 s spherical harmonic expansion of the anisotropic potential u The a superscript on Es indicates a complex conjugate. Details of the expansion and equations for E for specific interactions are given in Appendix C. In (312) and (314) J (p*,T ) is the single integral: 00 Jn(* ,T ) E dr 2 r2 g(r2) (316) 0 Values of the J integrals are tabulated elsewhere and have been fitted * to an empirical equation in p and T [37]. In (313) LY(1 ;n n ,) is the triple integral: * 0 r12+r13 LYrV/o *(nl) d *(n'l) L (E;nn') = ~ dr r dr3 r , f 12 12 13 13 0 0 r 12 r13 x dr23 r23(r12+ r13+ r23) goL(r12r13r23) P(cos al1) (317) Here P is the th order Legendre polynomial and a1 is the interior angle at molecule 1 in the triangle formed by r12, r13, and r23. Values of L are tabulated in Appendix E and have been fitted to an empirical equation * in p and T In (315) K (k'k";n n ,n,,) is the triple integral: * 0c oo r12+ r13 nn *(n1) *(n') K; n f d12 12 13 13 , 0 0 ir12 r131 *(n"l) * x dr23 r23 goL(r12r13r23) x (r12+ r13+ r23) ,R'"(a"l2a3) (318) The function nInii, is given by a spherical harmonic expansion: 1,,(el" 23) = I nm C(Z'rZ";mm'm") Ym(W12) mm m SY'm' (13) YI("m" (23) (319) where C(1V'";mm'm") is a ClebschGordan coefficient [38]. The Ykm are spherical harmonics in the convention of Rose [38]. In (318) and (319) the a. are the interior angles at molecule i in the triangle formed by 1 r12, r13, and r23. Expressions for ^p'i" for multiple interactions are given in Appendix A. F* Note that y3A(AA'A";ss's") given by (314) vanishes if (+'+1") is odd or if the molecules are linear and either (1+'+Z") or F* (2 +i2'+) is odd. Y2B(AA';ss') vanishes unless the anisotropic potential contains =0 spherical harmonics. F* F* F* F* Specific expressions for Y2A' Y2B' Y3A' and Y3B have been evaluated from (312), (313), (314), and (315), respectively, for various aniso tropic potentials. The results are tabulated in Appendix D. F The contributions to y given in (312,13,14,15) are simply related to the corresponding terms in the free energy expansion [33]: n +n ,l A (A;ss') F* T s A Y A(A;ss') = N (320) 2A 4 J NkT n nS, ;nn ,) A (AA';ss') y'F(AA';ss) = 8T 1 snss 2B (321) 2B 8 L (9,1 ;n n ,) NkT S S S F** T* Kn s'+n+n 1 A (AA'A";ss's") YF*(AA'A"s's'")' = T sn' ns3 (322) 3A 4 J n s n sl J ^ ( NkT y* = *T* KY(Zk'k";n n sn s1) A 3(AA'A";ss's") 3B s 8 K (k' ";n n ,n ,,) NkT The L and K integrals in (321) and (323) are defined by: * 00 oo r12+ r13 *(n1) *(n'1) * L(;nn') = dr2 r dr3 r dr r J 12 12r 13 13 *k23( 23 0 0 r 12 r131 x go(r12rl3r23) P(cos cl) (324) * 0o Co r12+ r13 *(n1) *(nt1) K(Z'R";nn'n") E dr12 r12 1) dr3 r13 , 0' 0 r 12 r13 *(n"l) * r23 23 (r12r13r23) A3"(ala2a3) (325) 3.2 Surface Tension Calculations for Model Fluids In this section model fluid calculations using the Pade per turbation theory for molecules obeying potentials of the form (31) are presented. For these calculations, the reference fluid surface F tension y was obtained from the Fowler model expression for a LennardJones fluid: YF 2/ = 3p *2[J5 2J1] (326) 0 5 11 where the J integrals are defined by (316) and are tabulated elsewhere [37]. Figures 4 and 5 show the effect of dipolar forces (A = 1 k 2 = 112) on the surface tension as predicted by the second order theory, the third order theory, and the Pade theory (292). The points on these figures are Monte Carlo calculations of the Fowler model surface tension for the LennardJones plus dipole fluids [39]. Figure 6 shows similar results for quadrupoles (A = 224). The second order and third order terms used in these calculations were determined from the expressions given in Appendix D. The results in Figures 4, 5, and 6 are similar to those found from the corresponding theories for the Helmholtz free energy [29,40]. Including the third order term extends the range of application of the expansion somewhat; however, the third order term overcorrects the second order theory for p* or Q* > 1. The Pade theory, on the other hand, interpolates between the second and third order theories. In the case of the polar fluids, the Pade agrees well with Monte Carlo results up to p* = 1.75. In Figure 7 comparison is made of the effects of various aniso tropies on the surface tension. The dipole and quadrupole curves are the Pade results from Figures 4 and 6, respectively. As in the case of bulk fluid thermodynamic properties, for a given value of the multiple strength (I* or Q ), the quadrupole potential is found to have a larger effect on surface tension than the dipole potential. 2.0 LL 10 0 Figure 4. 1.0 1/2 2.0 )L f(Ecr3)12 Fowler Model Surface Tension for a Fluid of Axially Symmetric Molecules Interacting with LennardJones plus Dipole Model Potential. po3 = 0.85, kT/E = 1.273 q ?I~D~JiC b 'a) r1 u cn '  (n 01% *) * C'4 4l 4 I 0 r cu O 0 0 ,q 0o 0 ( r4 (a C 0 r >il 0 H P4 0 a d oa  C; 0 w LL 0 1.0 Q/(E r5)1/2 Figure 6. Fowler Model Surface Tension for a Fluid of Axially Symmetric Molecules Interacting with LennardJones plus Quadrupole Model Potential. pa3 = 0.85, kT/e = 1.273 2.0 (V b 1.0L 1.0 0 1.0 Strength Constant (p*,Q*,K, or 6) Figure 7. Fowler Model Surface Tension for Fluids of Axially Symmetric Molecules Interacting with LennardJones plus Various Anisotropic Potentials. po3 = 0.85, kT/E = 1.273 The anisotropic overlap and dispersion results are from the second order F F theory, using expressions for Y2A and Y2B given in Appendix D. 3.3 Calculation of the Superficial Excess Internal Energy for Model Fluids To obtain the surface excess internal energy from the Pade expression (296), the temperature derivatives of the anisotropic contributions to the surface tension are required. Equations (312), (313), (314), and (315) give, respectively: F* ,in J y F*(A;ss') F* n +n 1 2A F* s s'1 1 = Y2A(A;ss' ) (327) T A T T F* Y Y2B(AA';ss') a* n L (;nns) 1 = Y2B (AA';ss') * (328) 3T T T YF* (AAA;ss's") [Rn J +n +n 3A(AAA s F/) _____s__s'__s"l 2_ __3A(A_'__';ss_ F* n ,+ns,,1 2 = 3A(AA'A"; ss's") (329) DT T T aY3B(AA'A'; ss 's") Fan KY(Ze''";n n ,n = Y3B (AA'A";ss's") (330) TT T T In Table 2 values for u from (295) and (296) are compared with computer simulation results for LennardJones plus quadrupole model fluids (see Chapters 6 and 7). The reference fluid surface excess internal energy, u was obtained from the Fowler model expression: TABLE 2 Test of the GibbsHelmholtz Equation in the Fowler Model Perturbation Theory for LennardJones Plus Quadrupole Fluids uFo2/ Q/(C5) 1/2 p3 kT/E Pad s MD 0.5 0.85 1.277 1.913 1.923 .016 0.707 0.931 0.765 2.670 2.656 .012 1.0 0.85 1.294 2.505 2.475 .019 u F 2/E = 2np*2 [J Jl (331) so 5 11 F where the J integrals are defined by (316). The values for us were n sa determined in the simulation by evaluating the ensemble average given in Chapter 5. In view of the demonstrated inapplicability of Equation (293) in the Fowler approximation for LennardJones fluids [31], the agreement between the theory and computer simulation shown in Table 2 is surprising. Since the inconsistency in the Fowler expressions for y and u is not limited to spherical potentials [6], the results in Table 2 s must be fortuitous. The agreement may, in part, be attributed to the high density state conditions considered, wherein the Fowler approxima tion is more accurate. Much of the agreement, however, must be due to cancellation of errors between the y and dy/dT terms in (293). 3.4 Surface Tension Calculations for Real Fluids The Pade perturbation theory developed in Chapter 2 has been used to predict pure liquid surface tensions for C02, C2H2, and HBr. In these calculations the reference fluid was taken to be a Mie (n,6) fluid. The anisotropies considered were the multiple interactions up through the quadrupolequadrupole term, as well as anisotropic dispersion and overlap contributions. 3.4.1 Mie (n,6) Reference Contribution to Surface Tension The Mie (n,6) fluid was taken as the reference in the perturbation theory calculations since Twu has determined values for E, a, and n by fitting perturbation theory calculations of liquid densities and pressures to experimental values along the orthobaric line for the fluids considered here [37]. It is felt that the test of the Pade theory for surface ten sion is strengthened by using these independently determined potential parameters. The Mie (n,6) potential u(n,6)(r) is given by [41]: 6 (n6) n 6 (n,6)r\ n 6(n6) ran 6 u (n6 (r) = n) (332) (n6) (n,6) To determine the surface tension y(n,6) for this potential, the Helmholtz 0 free energy of the nonhomogeneous, two phase, (n,6) fluid is expanded to first order in powers of n1 about the LennardJones (12,6) fluid free energy. The surface tension is obtained by applying the thermodynamic definition (21). Then, the Fowler approximation is made in order to obtain a form amenable to calculation. The expansion of the (n,6) free energy about the (12,6) may be done in two ways. In'one method, the values of c and a are taken to be the same for the two fluids. The expression for y(n,6) resulting from this expansion is: F(n,6) 2 F(12,6) 2 o Y0 o= 0 48Tp*2 [n 2{J(12,6) 12,6) + 6 H(12,6) =48w Jen 2{J J + 6 H E [n 11 5 11 x 1 (333) In (333) H12 is the single integral: H(2,6)( dr* *(n'2) (12,6) H (p,T ) = drn r n g (r ) (334) n 0 Values of this integral for n' = 11 are tabulated in Appendix K and * have been fitted to an empirical equation in p and T In the second method of doing the expansion, the values of E and r are taken to be the same for the (n,6) and (12,6) fluids. Here m r is the value of r where u(r) = E. The expression for yF(n,6) m o resulting from this second expansion is: F(n,6) 2 F(12,6) 2 0 0 To *2 (12,6) 1 (12,6) (12,6) = z 48np*2 [(1n2) J J + 6 Hl ] E n 11 2 5 11 x ( ] (335) When the (n,6) fluid is used as the reference, the second and third order terms y2 and Y3 in the perturbation theory involve integrals over (n,6) the (n,6) pair correlation function g (r). These (n,6) correlation o (12,6) functions can be related to LennardJones g (r) functions (for which o there are molecular dynamics results [32]) in the following way [42,43]: (n,6) (12,6),rep u(n,6),rep HS (n,6) g (r) g e y (d ) (336) o o and (12,6) (12,6),rep Bu(126)rep HS (12,6) (3 g (r) g e y (d (337) 0o where the superscripts rep and HS indicate the repulsive and hard sphere potential contributions, respectively. The function y is defined by: y(r) = eu(r) g(r) (338) In (336) and (337) the hard sphere diameter d is taken to be [42]: d = [1 eure(r)]dr (339) 0 Further, assuming yHS(d(n6)) HS (12,6) 340) (336), (337), and (340) give: (n,6) g(12,6) e[u (n,6),rep (12,6),rep] (341) g r) (r) e = (n,6),rep (12,6),rep where [u(n rep u(126)rep] vanishes for r > r . Using (341) with the Mie (n,6) potential in the integrals Jn,' Equation (316), and K(X'Y";nn'n"), Equation (324), Twu has found the resulting values to be negligibly different from those values obtained for the LennardJones (12,6) potential [37], at least for values of n close to 12. Hence, in the calculations reported here, the LennardJones (12,6) pair correlation function has been used in evaluating the terms Y2 and Y3 in the surface tension expansion. In the calculations reported here, Equation (335) was used to determine the reference fluid contribution. Values for the LennardJones F(12,6) term y were obtained from a corresponding states plot of surface tension for simple fluids, Figure 8. For the temperature range 0.6 kT/e < 1.24, the curve in Figure 8 obeys: S V Argon [44] A Krypton [45] O Xenon [46] 1.0  Methane [47] S \,v \ v 00 0 I I I I I 0.60 1.0 kT/E Figure 8. Corresponding States Plot for Surface Tension of Simple Liquids (12,6) 22 =.8950T*2 3.5177T* (12,6)/ = 0.8950T*2 3.5177T + 3.0166 (342) 3.4.2 Anisotropic Contribution to Real Fluid Surface Tension The calculations presented here include anisotropic dispersion, overlap, and multiple contributions up to the quadrupolequadrupole term: u = dis(202) + u dis(022) + udis (224) + u over(202) + u over(022) a dis dis dis over over + uoverdis(202) + u overdis(022) + udisQ(224) overdis overdis disQ + u m(112) + u mul(123) + u (213) + u (224) mul mul mul mul (343) wherein the subscripts dis, over, Q, and mul refer to dispersion, overlap, quadrupole, and multiple, respectively. Appropriate parts of this model have been used by FlytzaniStephanopoulos et al. [33] and Twu [37] in studies of bulk fluid thermodynamic properties. The multiple contribu tions to surface tension for this model potential are obtained by combining (38), (39), (310), and (311) with (312), (313), (314), and (315), respectively: F* F* F* F* 2A = 2A(112) + 2 2A(123) + Y2A(224) (344) F* F* F* 3A = 3 3A(112;112;224) + 6 Y3(112;123;213) F* F* + 6 yA(123;123;224) + Y3A(224;224;224) (345) F* F* F* Y3B = 3B(112;112;112) + 3 Y3B(112;123;123) F* F* + 3 (123;123;224) + y3B(224;224;224) (346) 3B 3B Expressions for the terms on the right side of (344), (345), and (346) F* are given in Appendix D. The y2B term is zero for multipoles since only terms with A # 0 spherical harmonics occur in the multiple potentials. The dispersion and overlap anisotropies in (343) have been included in only the second order term Y2 in calculating the surface tension from the Pade perturbation theory (292). The inclusion of dispersion and overlap contributions to the third order term y3 requires evaluation of difficult multibody terms. Expressions for the Y2A and Y2B terms for anisotropic overlap and dispersion are given in Appendix D. 3.4.3 Results for Real Fluids Figure 9 compares the Pade predictions of surface tension for CO2 with experiment, while Figure 10 shows a similar comparison for C2H2 and HBr. The experimental values for surface tension of CO2 and C2H2 were taken from the compilation by Jasper [48]. The corresponding values of saturated liquid densities for CO2 and C2H2 were taken from Vargaftik [49]. Experimental values of surface tension and saturated liquid density for HBr were obtained from Pearson and Robinson [50]. Values of the potential parameters for C02, C2H2, and HBr were taken from Twu [37] and are listed in Table 3. The deviations in surface tension between theory and experiment are less than 10% for CO2 and less than 12% for C2H2 and HBr for the temperatures shown in the accompanying figures. The consistent deviations between theory and experiment, especially for C2H2 and HBr, suggest that adjustment of the potential parameters would improve the agreement. It is not a very informative test of the theory, however, to adjust potential I I I I m.p. II 1.0 kT/E Surface Tension for CO2 comparing Theory Calculations (points) with Values (line) c.p. I l 1.2 Perturbation Experimental 0.6 (cJ b 0.4 0.21 0 L 0. 3 Figure 9. " HBr C2H2 Figure 10. 0.8 kT/E Surface Tension for C2H2 and HBr Perturbation Theory Calculations with Experimental Values (lines) b b rZ 1.0  0.8 0.6 0.9 comparing (points) TABLE 3 Potential Parameter Values Used in Calculating Surface Tension Fluid c/k (1018) Q(1024) 6 K (K) (A) n (esu cm) (esu cm2) CO2 244.31 3.687 16  4.30 [51] 0.1 0.257 C2H2 253.66 3.901 16  5.01 [52] 0.3 0.270 HBr 248.47 3.790 12 0.788 [51] 4.0 [51]  All values for , u, n, 6 and K are taken from Twu [37]. parameters using surface tension in order to calculate surface tension. The deviations between theory and experiment shown in Figures 9 and 10 do not increase much with temperature, and, in fact, for CO2 the devia tions decrease. This is in contrast to the results found for simple LennardJones fluids wherein surface tension calculated in the Fowler model show rapidly increasing disagreement with experiment as the temperature is increased [2,31]. Further, the Fowler model values of surface tension for LennardJones fluids are generally larger than the experimental values at the same temperature. It may be that, in addi tion to the questionable adequacy of potential parameter values used here, the Pade approximant in some way corrects for errors introduced in using the Fowler model for the interface. There is some evidence for this from the Pade values for the surface excess internal energy presented in Section 3.3. 3.5 Correlation of Surface Tension for Pure Polyatomic Liquids The perturbation theory for surface tension presented in Chapter 2 has been used as a basis for correlating surface tensions of a large number of polyatomic liquids. The perturbation theory gives the surface tension in the Fowler approximation as: F F F F S+ Y + Y + + + Y (347) Y= YO +2A 2B 3A 3B A simple correlation for y may be obtained by making the van der Waals type assumption that the reference fluid radial distribution function is a constant: goL(r 2) = c (348) Using (348) in (287), (288), (289), and (290), together with a similar approximation for the triplet correlation function, go(r12rl3r23), (347) becomes, in reduced form: *2 *2 *3 *4 aPL a2PL a3PL a4PL Y = Y + + + + (349) o *2 *2 T T T T where a' E dr r (350) 1 4 12 12 a 12 2 0 * 2 12 13 a f f * a dr2 r12 dr r dr r 2 2 12 113 23 23 a a Wl2 3 2* ) 0 0 r12 r1 x (r1 + r + r23) (351) 00 a d *3 *3 2 a 2 j dr1 *r1 (352) 3 12 12 a 2 0 r* + r* T2 012 13 S* * a E dr r dr r dr r 4 6 12 12 13 13 23 23 0 0 Ir12 r131 ** * x (r12 + r13 + r23) (353) a a a W l 2 3 12 13 23 Equation (349) can be written in an equivalent form using the critical constants, Tc, Vc, and pc as reducing parameters rather than the potential parameters, E and a: 2 3 2 3 alPR a2PR a3PR a4PR Y = YoR + R + + + (354) R R 2/3 1 where R (Vc/NA)2/3 (kT (355) PR = P/Pc (356) TR = T/Tc (357) In transforming (349) into (354), the proportionality of the potential parameters to the critical constants is obtained by the usual correspond ing states method [41]. Here, however, the potential is for polyatomic fluids and, therefore, contains parameters in addition to the energy and distance parameters and a, e.g., the multiple moments, p, Q, and anisotropic polarizability, K, and overlap parameters 6. For such an intermolecular potential, if the usual derivation of corresponding states theory is followed [41], one finds: kT c * Tc = C1( ,Q,,K,...) (358) 3 * Pc = Pco = c2(H 'Q ,',K,...) (359) PO c * P = = c3(H 'Q ,6,K,..) (360) where cl, c2, and c3 are constants only if the anisotropic potential para * meters p Q 6, K, *** are kept fixed. However, these "constants" can be absorbed in the terms al, a2, etc., as has been done in (354). Thus, 2 al = a'c2/cl, etc. The transformation from the potential parameters E and a as reducing parameters to critical constants as reducing parameters may then be made in the usual way [41]. 61 Equation (354) has been used as a basis for correlating surface tension by treating the parameters a. as semiempirical constants. Various 1 truncated forms of (354), including its Pade form analogous to (292), have been tested against experimental surface tensions for numerous poly atomic liquids. The form giving the best comparison with experiment was found to be that terminated after the Y2B term: 2 OR YR = R + (al + aR) (361) Values for al and a2 for use in (361) have been determined for numerous polyatomic liquids by least squares fitting to available experimental data. The resulting values are listed in Table 4. In these calculations, the reference contribution was obtained from a fit of reduced surface tensions of the inert gases and methane, analogous to the curves in Figure 8. YR = 2.4724T2 7.5918T + 5.0748 (362) oR R R which applies for 0.4 T TR 0.95. The validity of the correlation (361) may be tested by determining 2 whether experimental data give a linear relation between (y YoR)TR/P and pR' as implied in (361). Such a test has been conducted for several polyatomic liquids, and results for carbon dioxide, acetic acid, and methanol are shown in Figures 11, 12, and 13, respectively. These figures are typical results for small to moderately large polyatomics and indicate a satisfactory correlation. The corresponding comparisons TABLE 4 Values for the Parameters a1 Surface Tension Correlation of and a2 in the Equation (361) Equation (361) Substance Range a (102) a2(102) References T P Y r Paraffins Ethane Propane nButane iButane nPentane iPentane nHexane nHeptane nOctane iOctane nNonane nDecane nDodecane nTridecane nTetradecane nPentadecane nHexadecane nHeptadecane nOctadecane nNonadecane nEicosane Cycloparaffins Cyclopentane Methylcyclo pentane Ethylcyclo pentane Cyclohexane Methylcyclo hexane .43.59 .50.77 .50.57 .52.60 .54.67 .59.66 .54.93 .54.95 .48.90 .50.67 .46.63 .44.60 .41.57 .40.55 .41.54 .40.53 .41.55 .41.54 .40.52 .41.52 .40.51 .55.61 .53.59 .50.55 .51.62 .48.65 6.880 14.17 8.548 7.749 6.504 6.250 2.398 1.860 1.335 8.536 7.340 8.699 11.42 12.06 12.21 12.83 11.59 11.10 10.25 11.15 11.15 2.246 5.606 2.469 2.082 1.498 1.451 1.918 1.910 0.7617 1.942 1.230 1.576 2.262 2.360 2.339 2.452 1.893 1.720 1.439 1.596 1.539 5.689 1.339 5.973 1.288 10.65 2.832 6.249 1.351 49 2.885 0.4730 49 Reduced temperature range over which al and a2 were fitted. Sources of experimental data for liquid tensions y. densities p and surface 49 48 49 48 49 48 TABLE 4 (Continued) Substance Range a1(102) a2(102) References T P Y r Olefins Propylene .53.67 7.143 2.173 49 48 1Butene .48.70 4.752 0.9880 49 48 2Butene .51.65 3.606 0.8401 49 48 1Hexene .54.64 8.630 2.142 49 49 1Octene .47.56 8.047 1.876 49 49 Cyclopentene .56.62 5.445 1.195 49 49 Aromatics Benzene .48.93 0.4634 0.5432 53 11 Toluene .46.63 8.130 2.080 49 49 Ethylbenzene .44.60 9.602 2.384 49 49 Isopropylbenzene .46.57 8.175 1.851 49 49 Alcohols Methanol .53.92 11.60 4.971 49 49 nPropanol .55.68 25.56 8.854 49 48 iPropanol .56.60 33.13 10.98 49 49 nButanol .48.54 20.35 6.624 49 49 Organic Halides Methyl Chloride .68.73 1.215 0.1476 53 48 Ethyl Bromide .56.60 15.00 4.565 53 48 Carbon Tetra .49.89 0.9412 0.4869 49 49 chloride Chlorobenzene .43.88 2.138 0.000452 49 49 Oxides Carbon monoxide .61.68 8.502 2.748 49 48 Carbon dioxide .71.95 1.587 1.395 49 49 Water .44.58 34.75 11.90 54 48 Others Acetic Acid .49.86 5.984 2.906 49 49 Acetone .54.69 6.178 1.783 49 49 Ammonia .49.58 6.878 2.269 53 48 Aniline .39.65 13.38 3.976 49 49 Carbon disulfide .51.58 5.968 1.869 Chlorine .47.57 4.103 1.340 49 48 Diethyl ether .59.95 0.4972 1.227 49 49 Ethyl acetate .52.90 0.08686 1.058 49 49 I I I I 18 *0 0 rfor CO2 o 14  10  2.0 2.4 PiPc Figure 11. Test of Surface Tension Correlation (line) for C02 20 OO I cr cL 0 I 2.1 Figure 12. Test of Sur Acetic Acid 2.5 3.0 /face Tension Correlation (line) for *face Tension Correlation (line) for O 2.0 Figure 13. 20. Figure 13. 2.5 P/Pc Test of Surface Methanol Tension Correlation (line) for 3.0 of predicted surface tensions with experimental values for these and several other liquids are shown in Figures 14 and 15. The linear relation suggested by (361) is not obeyed by experi mental data for long chain hydrocarbons, however. Typical plots are shown in Figure 16. In spite of the poor correlation for these substances, the predicted surface tensions using (361) were usually within 3% of the experimental values as shown in Figure 17. Generally, the correlation of (361) reproduced the experimental data for substances tested here within 3% for values of TR < 0.92. In order for the correlation to apply in the critical region, we must have: Rl YoR = 0 (363) C.P. C.P. hence, from (361): a = a2 a (364) Then 2 r R + a (1 p) (365) R in the critical region. This relation was not obeyed by many of the liquids in Table 4. 2.5 1.5 0.5 0.4 Figure 14. 0.6 TT.8 T/Tc Comparison of Surface Tension Calculated from the Correlation (lines) with Experimental Values (points) for Several Polyatomic Liquids 2.0 r 1.0  0 0. 0 CC14 * CO2 A MeOH A BuOH V PrOH V H20 4 0.6 Figure 15. T/Tc 0.8 Comparison of Surface Tension Calculated from the Correlation (lines) with Experimental Values (points) for Several Polyatomic Liquids _I __ 