Surface tension and computer simulation of polyatomic fluids

MISSING IMAGE

Material Information

Title:
Surface tension and computer simulation of polyatomic fluids
Physical Description:
xxxi, 376 leaves : ill. ; 28 cm.
Language:
English
Creator:
Haile, James Mitchell, 1946-
Publication Date:

Subjects

Subjects / Keywords:
Liquids   ( lcsh )
Surface tension   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis--University of Florida.
Bibliography:
Includes bibliographical references (leaves 357-375).
Statement of Responsibility:
by James Mitchell Haile.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
oclc - 04102852
ocm04102852
Classification:
lcc - QD541 .H13
System ID:
AA00011132:00001

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 open-ended 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 Chorng-Horng 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 vapor-liquid interfacial

density profile for Lennard-Jones 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 Shien-Shion 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 VAPOR-LIQUID DENSITY-ORIENTATION 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 Full-Size
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
1l-2m )..










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 Site-Site 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 Density-Orientation
Profile of Polyatomic Fluids....................... 261

8.3 Computer Simulation of Polyatomic Fluids........... 262

APPENDICES:

A EXPRESSIONS FOR THE ANGLE AVERAGES IN EQUATIONS (3-4)
TO (3-7)................................................ 267

B COORDINATE TRANSFORMATION AND INTEGRATION OVER EULER
ANGLES TO OBTAIN EQUATIONS (2-89) AND (2-90)............ 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 (6-77).......................... 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 SITE-SITE 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 (1-1)......... 5

2 Test of the Gibbs-Helmholtz Equation in the Fowler
Model Perturbation Theory for Lennard-Jones 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 (3-61).................... 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 Lennard-Jones + 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 Lennard-Jones plus Quadrupole
Fluid Obtained in this Work and Compared with those
given by Berne and Harp................................... 177

18 Equilibrium Properties for Lennard-Jones plus Quadru-
pole Fluid at pa3 = 0.85, Q/(EC5)1/2 = 1/2............... 181

19 Equilibrium Properties for Lennard-Jones plus Quadru-
pole Fluid at po3 = .931, Q/(Ec5)1/2 = 0.707............. 182

20 Equilibrium Properties for Lennard-Jones plus Quadru-
pole Fluid at po3 = 0.85, Q/(Eo5)1/2 = 1.0............... 183

21 Anisotropic Contributions to Equilibrium Properties for
Lennard-Jones plus Quadrupole Fluid at po3 = 0.85,
Q/i(E5)1/2 = 1/2 ........................................ 188

22 Anisotropic Contributions to Equilibrium Properties for
Lennard-Jones plus Quadrupole Fluid at po3 = .931,
Q/(E~ 5)1/2 = 0.707....................................... 189

23 Anisotropic Contributions to Equilibrium Properties for
Lennard-Jones plus Quadrupole Fluid at po3 = 0.85,
Q/(Ea 5)1/2 = 1.0......................................... 190

24 Equilibrium Properties for Lennard-Jones plus Overlap
Fluid at po3 = 0.85, 6 = 0.10............................ 192

25 Equilibrium Properties for Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 (6-25)
and (6-34) .............................................. 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones 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 Site-Site Correlation Function for Lennard-Jones plus
Quadrupole Fluid with /c = 0.3292, kT/c = 1.277,
po3 = 0.85, Q/(co5)1/2 = 0.5.............................. 350

J2 Site-Site Correlation Function for Lennard-Jones plus
Quadrupole Fluid with /C = 0.2955, kT/E = 0.765,
PO3 = 0.931, Q/(ao5)1/2 = 0.707......................... 351

J3 Site-Site Correlation Function for Lennard-Jones plus
Quadrupole Fluid with i/a = 0.3292, kT/E = 1.294,
po3 = 0.85, Q/(Ec 5)1/2 = 1.0....................... ..... 352

J4 Site-Site Correlation Function for Lennard-Jones plus
Anisotropic Overlap Fluid with /a = 0.3292, kT/E = 1.287,
po3 = 0.85, 6 = 0.10 ..................................... 353

J5 Site-Site Correlation Function for Lennard-Jones 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 Vapor-Liquid 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 Lennard-Jones
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 Lennard-Jones
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 Lennard-Jones
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 Lennard-Jones
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 n-Hexane and
n-Octane................................................ 70

17 Comparison of Surface Tensions Calculated from the
Correlation with Experimental Values for Several
Hydrocarbons... ......................................... 71

18 Interfacial Density Profile for Lennard-Jones Fluid..... 80

19 Interfacial Density-Orientation 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 Density-Orientation 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 Density-Orientation 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 Lennard-Jones Fluid........................... 87

23 Interfacial Density-Orientation 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 Density-Orientation 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 Density-Orientation 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 Density-Orientation 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 Center-Center Pair Correlation Function for a
Lennard-Jones 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 Lennard-Jones plus Dipole Model
Fluid at Primary Pair Orientations....................... 161

35 Pair Potential for Lennard-Jones plus Quadrupole Model
Fluid at Primary Pair Orientations....................... 165

36 Surface of the Lennard-Jones plus Quadrupole Pair
Potential for the Tee Orientation as a Function of
the Quadrupole Strength .................................. 166

37 Pair Potential for Lennard-Jones 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 Lennard-Jones 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 Lennard-Jones plus Anisotropic
Overlap Model Fluid at Primary Pair Orientations ........ 170


xvii










Figure Page

40 Mean-Squared Displacement of Molecular Centers of
Mass for Lennard-Jones plus Quadrupole Fluid............... 173

41 Fluctuation in Temperature for Lennard-Jones plus
Quadrupole Fluid........................................... 175

42 Fluctuation in the Ratio of Translational to Rotational
Kinetic Energy for Lennard-Jones plus Quadrupole Fluid..... 176

43 Effect of Quadrupole Moment on the Center-Center Pair
Correlation Function at pa3 = 0.85.......................... 194

44 Spherical Harmonic Coefficients g2Z for Lennard-Jones
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 Lennard-Jones
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)]r-3 for G2 and Ua'
respectively, for Lennard-Jones plus Quadrupole Fluid...... 214


xviii










Figure Page

54 Angular Pair Correlation Function for the Lennard-Jones
plus Quadrupole Fluid of Figure 44 for the Tee Orienta-
tion (01 = 90, 82 = 0, ( undefined)...................... 217

55 Angular Pair Correlation Function for the Lennard-Jones
plus Quadrupole Fluid of Figure 44 for the Cross and
Parallel Orientations (01 = 02 = 900)..................... 218

56 Angular Pair Correlation Function for the Lennard-Jones
plus Quadrupole Fluid of Figure 44 for a Skewed Orienta-
tion (1 = 82 = = 45 )........... .. ........................ 220

57 Angular Pair Correlation Function for the Lennard-Jones
plus Anisotropic Overlap Fluid of Figure 49 for the Tee
Orientation (81 = 900, 82 = 0, < undefined)............... 222

58 Angular Pair Correlation Function for the Lennard-Jones
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
Lennard-Jones 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
Lennard-Jones 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
Lennard-Jones.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
Lennard-Jones 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
Lennard-Jones 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
Lennard-Jones 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 Lennard-Jones plus Quadrupole Fluid with
Q/(Ca5)1/2 = 1.0 ......................................... 234


xix











Figure Page

66 Surface of the Angular Pair Correlation Function for
Lennard-Jones 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
Lennard-Jones 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
Lennard-Jones 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 Site-Site Pair Correlation Function for Lennard-Jones
plus Quadrupole Fluids................................... 245

73 Possible Square Packing of Lennard-Jones plus Quadru-
pole Molecules for Interpreting g (r) .................... 247

74 Site-Site Pair Correlation Function for Lennard-Jones
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 Lennard-Jones
Molecules in the Volume Element Sampled in the Filmed
Animation.................................................. 256

77 Frame from the Filmed Animation of Lennard-Jones Mole-
cules Corresponding to the Sixth Time-Step in the
Molecular Dynamics Calculation............................ 257

78 Frame from the Filmed Animation of Lennard-Jones Mole-
cules Corresponding to the 101st Time-Step 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 (4-26)

B(zl,Zl2) = Function defined in Equation (4-27)
R
CV = Residual contribution to constant volume heat capacity


2 ;mlm2m) = Clebsch-Gordan 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 (2-32)

F2A,F2B = Functions defined in Equation (2-62) and (2-63),

respectively

F3A,F3B = Functions defined in Equations (2-77) and (2-78),

respectively

GL = Angular correlation parameter

H = Integral defined in Equation (3-34)


xxii










I = Moment of inertia

I(r) = Functions defined in Equations (6-65) and (6-69)

I = Integral defined in Equations (2-40) and (B7)
z

Inn = Integral defined in Equation (Gl)

J = Integral defined in Equation (3-16)
11 2
J = Integral defined in Equation (6-89)
n
K = Coefficient of ellipticity for plane polarized light

K = Integral defined in Equation (3-25)

K = Integral defined in Equation (3-18)

KL,KV = Functions defined in Equations (4-30) and (4-31),

respectively

L = Angular momentum operator

L = Integral defined in Equation (3-24)

L = Integral defined in Equation (3-17)

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 (5-8)

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 (3-50)

to (3-53)

c = Constant defined in Equations (4-21) and (4-22)

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 (3-39)

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)


= Center-center pair correlation function

= Site-site 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 (4-21) and (4-22)

= Exponent of r in repulsive part of Mie potential,

Equation (3-32)

= 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

= x-component of r2
--12
= y-component of r2
--12
= Function defined in Equation (3-38)

= z-coordinate location of molecule 1

= z-component of r12

= Maximum value of zl, defined in Equations (2-42) 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 (4-26) and (4-27)

ai = Interior angle in the triangle formed by rl2, r13, r23,

at molecule i

ai,8i = Constants in the predictor-corrector algorithm; values

given in Equations (6-45) and (6-46), 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 density-orientation 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 (3-19)

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 l-12-23 )

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 Lennard-Jones 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 density-orientation 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 Lennard-Jones 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 Lennard-Jones plus quadrupole and


xxx










Lennard-Jones 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.

Site-site 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,

x-ray, 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 fluid-fluid 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 (1-1)
-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 (1-1) are given in Table 1. Equation (1-1) 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

fluid-fluid 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 Vapor-Liquid Interface










TABLE 1

Examples of Macroscopic and Microscopic
Interfacial Properties Related by Equations

of the Form (1-1)





00
X = k [MB M(z)]dz



0 CO

S(z) L ]dz + [i(z) p ]dz (1-2)
-f 0
000

Y = j [PN pT(z)]dz (1-3)
-0T

7T 2 P (z) P (z)
z x
K = (n + dz (1-4)
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 fluid-fluid

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

Kirkwood-Buff 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 (1-5)
-00










The Kirkwood-Buff 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 Kirkwood-Buff 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) (1-6)



where 6 is the unit step function,



1 if x > 0

6(x) = (1-7)

0 if x < 0


subscript L indicates a bulk liquid property and the negative z

direction is into the liquid. The resulting Fowler-Kirkwood-Buff

expression for surface tension is:

2 m
P u 4 du(r2)
Y =8 dr12 12 dr12 gL(r2) (1-8)
0

As could be expected, the Fowler-Kirkwood-Buff 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) (1-9)
-00 -00 -CO


This relation is more general than the Kirkwood-Buff 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 (1-9) has not been as thoroughly

studied as the older Kirkwood-Buff formula.

The Kirkwood-Buff 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

(1-5) and (1-9). 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 Born-Green-Yvon 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 density-orientation 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 2-4 and is concerned with the theory of vapor-liquid

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

vapor-liquid interfacial density-orientation 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 Lennard-Jones 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 (1-9).

The second is a generalization of the Kirkwood-Buff expression (1-5)

which has been previously derived [10,22]. In this section the deriva-

tion for the general Kirkwood-Buff equation is summarized and its

simplified form obtained when the Fowler approximation is made.


2.1.1 Generalized Kirkwood-Buff 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= (2-1)
Y NVT


Since only the configurational part of the free energy Ac depends on

the interfacial area, (2-1) may be written as:



y = kT nn Z T (2-2)
NVT










where Ac = kT nn Z (2-3)



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 ) (2-4)
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

body-fixed reference frame on the molecule and a space-fixed frame

located external to the system. Substituting (2-4) into (2-2):


kT f rr N N -U(rNWN))
Y = Z dr dw e JNVT (2-5)



The differentiation in (2-5) 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' (2-6)

V ,
z =-z


Performing the differentiation in (2-5) and changing back to the old

variables gives:










N N N NN
Y 1 f f e-8(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

(2-7) 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(N-1)
f( l21 l2) = Z f


N N
-U(r w )
dr 3* dr dw3. dwN e -
-3 -N 3 N


Integrating (2-9) 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, (2-11) becomes:


(2-7)


(2-8)






(2-9)


(2-10)


(2-11)


(2-12)


(2-13)









Y =I dzl


drl2 dwldw2 f(zlrl2l2) rl


au(12)
r12


312 Du(12)
12 8Z
12


(2-14)


Defining the angle average by:


1Wm2


(''')dw dw2


Equation (2-14) can be written:


I dr12

au(12)
Dr12


au(12)
- 3zl2 (z l2 ]>
z12 z 12J 1)2
(2-17)


Equation (2-17) is one form of the general Kirkwood-Buff equation.

Another form may be obtained by transforming (2-17) to spherical


coordinates.


Sau(12)
z 12'Y12


Since [25]


= cos 612


(2-18)


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 (2-19) into (2-17) gives:


12 -
2 j


where


S= f dw


(2-15)


(2-16)


co
f2 f ,
S= 4- jdz1
-0-o


fau(12)
I r12


sin 012
. rl2


(au(12)]
a "12


(12


8u(12)
ar12


(2-19)










Y = YA + YB (2-20)



YA 2 dz drl2 P2(cos 12) r2 (2-21)
-2 ur12 12
-o00
Y = 4 dz I dr12 sin 812 cos 012 (2-22)
B 4 21 1-212 H2 I 2
-00


Equations (2-20,21,22) give the general Kirkwood-Buff 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 Kirkwood-Buff equation (1-5).


2.1.2 Fowler-Kirkwood-Buff Equation for Nonspherical Molecules

The general expression for the surface tension (2-20) is in-

tractable as it stands due to the unknown distribution function

f(zl r2-wl2). 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) (2-23)



where 0 is the unit step function as in (1-7), pL is the bulk liquid

density, and gL is the bulk liquid angular pair correlation function.

Introducing the Fowler model (2-23) into (2-21) allows the integration

over z and w12 to be performed giving:








2 00
F L T 4 u(12)
YA = 8 dr2 r12 (2-24)
A 8 112 1lm2
0

where the superscript F indicates Fowler model. Similarly, (2-22)

reduces to:


Y2 dr r (2-25)
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 Fowler-Kirkwood-Buff

equation for nonspherical molecules is:

2 20
F L 4 au(12)
y f dr12 12 (2-26)
8 1212 L r12 1Wl2
0

For the special case of spherical molecules, (2-26) reduces to

the usual FKB expression (1-8).


2.2 General First Order Perturbation Theory
for Surface Tension

2.2.1 Rigorous First Order Term

The difficulty with using the general Kirkwood-Buff expression

(2-20) lies with the, in general, unknown interfacial distribution

function f(zl-r2wlw2). 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) (2-27)


where X is a perturbation parameter such that when A = 0, (2-27) gives

the reference potential and when A = 1, (2-27) gives the full potential.

Expanding the two-phase system Helmholtz free energy (2-3) in powers of

A and setting X = 1 gives:



A = A + A1 + *" (2-28)


where A1 = drldr2 (2-29)
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 (2-1) to (2-28):



Y = Y + Y1 + *** (2-30)


where Y= I ( dr dr F(r2) (2-31)
as -1 -2 1 1--12) INVT



and F1(zll2) E ( <1f ( 2 u(12;A) > (2-32)
l A=0 wlw2


Applying Green's method as in section (2.1), integrating over xl and

Y,' and transforming rl to r12, (2-31) becomes:





20


1 f f 12 flzs
I = 2 dl r S[f 2NVT



S 3Fl(zr-l2)/@S may be evaluated by considering:


12 ,F (zlrl2) 3 I (zl 1 2 z1
=S + S
as -r12 a1z as


-12 3Fl(z 1-2)
2 r 12
-12


3 z Fl(z l12)
2 12 z12
aFl(Zlr12)
Z F (zIl-
- zz


Thus, (2-33) 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 (2-36) 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 (2-37) 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

(2-48) below, then Y1 vanishes.


(2-33)


(2-34)



(2-35)


(2-36)


(2-37)


SFl(Z-142)
S










2.2.2 First Order Term in the Fowler Model


To make (2-37) amenable to calculation, the Fowler approximation

(2-23) 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

(2-38)


Substituting (2-15) and changing the order of integration:


2
L
2Q4


dr dw dw
-12 1 2


fau(12;A)J
( ax X=O


(2-39)


where Iz I dzl Z1


(2-40)


The integral I may be evaluated by parts to give:

z

I = go,(12) zmax
z oL max


z =
max
-r 12


cos e12


(2-41)


(2-42)


if 6 >2

if 12

if e < 7T
12 2


and 012 is the spherical coordinate polar angle

Figure 3. Equation (2-39) becomes, therefore:


for r12' as shown in


Sdwdw2 u(12;A)
1 2 ( @A ]=


goL(12) zmax (2-43)


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 (2-15) and noting that g(12) in the

bulk liquid is independent of 12, (2-43) can be written as:

2 o
F _L 2 Y = + dr r < g (12)> d9 z (2-44)
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)> (2-45)
1 22 12 12 X =0 oL (242
0

Equation (2-45) 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 (2-30) is in-

creased to higher order, the second order term is found to include a

term containing the reference four-body 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 (2-46)


With this choice of reference, (2-27) becomes:


u(12;X) = u (rl2) + Aua(12)


(2-47)










(2-46) and (2-47) together give the simplification:


a3u(12; A) >
X I--=0 loW2


=
1 2)0


= 0


(2-48)


Then (2-29) and (2-37) have:


(2-49)


Al = Y1 = 0


and the second and third order terms (A2,A3,Y2,y3) simplify. Thus, in

the Pople expansion, (2-28) and (2-30) become, to third order:


AC = AC + A + A (2-50)
o 2 3


Y = YO + Y2 + Y3 (2-51)



2.3.1 Derivation of Second Order Terms, Y2A and Y2B

The second order term in (2-50) for the bulk phase liquid is [27]:


2
A= dr dr
A2 4 -2 a 12


(2-52)


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)


(2-53)










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

(2-53) vanish. Such potential models as anisotropic overlap and dispersion

contain X = 0 spherical harmonics, in which case the second term must be

included.

Equation (2-52) may be written, therefore, as


A2 = A2A + A2B


(2-54)


where


2A dr dr go(r ) 2A 4 -2 (r12) a l> 2


(2-55)


p3
A2B = dr dr2d go(r12r13r23)


r~, r2 and r3 are each integrated

indices and (2-56) may, hence, be


pN3 dr dd
2 drd2dr3 o(r12r13r23)


l2m3

(2-56)

over in (2-56), the indices are

written as:


i23


(2-57)


For a fluid nonuniform in the z-direction, (2-55) and (2-57)

generalize to:



A = dr dr p (zl o(z2) o(12 (2-58)
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

(2-59)

The corresponding terms in the expansion for surface tension (2-51)


may be found by applying (2-1) to (2-58) and (2-59):



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) l2



2B Po(zPP(z2)Po(z3) go(z12-r13) l2w 3


(2-60)



(2-61)


(2-62)



(2-63)


Equations (2-60) and (2-61) have the same form as (2-31); thus they may

be evaluated in a similar manner to give, analogous to (2-37):


Y =
2A 4


Y = f
2B 2


dr12 dz z 2A-2
-L12 1 1z 1


dr12d13 dz
-12 -13 1j


aF2B (Zl112)
z 1 zI


Substituting (2-62) and (2-63) into (2-64) and (2-65), respectively,


and using the relations:


(2-64)


(2-65)









fo(z12) = o(zl)Po(z2) go(zlr2) (2-66)


fo(Zlr12r13) = Po(zl)po(z2)po(z3) go(Zlr12r13) (2-67)


gives,


S d 2 rfo(zlr12)
Y2A dr dzl z (2-68)
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
(2-69)

The derivation of (2-68) and (2-69) 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 (2-53).

2.3.2 Derivation of Third Order Terms, Y3A and Y3B

The third order term in (2-50) for the bulk phase liquid is:


A dr dr (2-70)
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 (2-56), and generalizing to a

fluid which is nonuniform in the z-direction, (2-70) 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.(1z-12-13) x






Applying (2-1) to (2-71) gives:


Y3 = Y3A + Y3B


dr dr


F3A( -12 NVT
NVT


g2
Y3B =6 7 dr rd2d3r F3B(zl-rl2l3)
NVT





F3A(zlrl2 Po(zl)Po(2) go(z 12)

(2-75)




(2-76)


(2-77)


1 o3(12)u (13)u (23)>
F3B(zll2rl3) Po(zl)Po(z2)Po(z3) go(zl23) 12 3

(2-78)

Again, (2-75) and (2-76) have the same form as (2-31) and,

consequently, they may be evaluated in the same manner to give


(2-71)



(2-72)


(2-73)


where


(2-74)


B2
Y = 12
3A 12










analogous to (2-37) (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- (2-79)
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

(2-80)

where fo(z1-12) and fo(zl2E1- 13) are given by (2-66) and (2-67), respec-

tively.

The derivation of (2-79) and (2-80) 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) (2-81)


fo(z!121 = (-z)(-z2) (-z3 goL(r2) (2-82)



Putting (2-81) into (2-68), (2-69), and (2-79) for y2A' Y2B and Y3A'

respectively, and using (2-82) in (2-80) for Y3B' all give terms analo-

gous to (2-39) for y1 in the Fowler model. In each case, the integration

over zI may be done by parts giving, analogous to (2-44):








2 -
F PL r 2 2 f
Y1A dr r2 < (12) g(rL ) d z (2-83)
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 (2-84)

2 2
Y dr 2 (r f d 12 2-85)
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 (2-86)



In the case of the two-body terms, Y2A and Y3A' the values for

z maxare given by (2-42); see Figure 3. Hence, the integration over

S12 can be performed giving:

2 co
2A = drl2 r2 goL(rl2) 2 (2-87)
0

F Lp 3 3
YA = dr r2 g oL(rl2) (2-88)
3A 12 12 12 goL 12 a W1 2(2-88)
0


The integration over 12 and w13 in the three-body 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) (2-89)


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) (2-90)
a1a 12w3 12


F F F F
Computationally convenient forms for 2A' 3A 2B and y3B are

derived from (2-87), (2-88), (2-89), and (2-90), 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 (2-50) in the form of a Pade approximant:


A
Ac = A + (2-91)
o A
3
1
A2


The analogous form for the surface tension expansion (2-51) is:






32





1
Y2



Calculations based on (2-92) 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 Gibbs-Helmholtz equa-

tion:



Us = Y T l(I (2-93)
NV


where u is the surface excess internal energy, U per unit of surface

area S:


U
u = (2-94)
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 (2-93) [30,31]; Freeman and McDonald show

that these two expressions give quite different results for use in the

case of Lennard-Jones liquids [31]. To determine the validity of (2-93)

for the Pade perturbation theroy for surface tension presented above,










we combine (2-92) and (2-93) to obtain:



u = u + u (2-95)
s so sa


where u is the isotropic reference contribution, and
so

DY3 2 3
Y T T 1-2 2
u = (2-96)
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) (3-1)



where u is the Lennard-Jones potential,



Uo(rl2) = 4E[(/rl2)12 (a/r12)6] (3-2)



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) (3-3)


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 (3-1); 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 (2-87,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) (3-4)
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) (3-6)
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)


(3-7)


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 (3-4,5,6,7) have been evaluated [33] and

are tabulated in Appendix A. Substituting those expressions into

(3-4,5,6,7) gives:


F* r F*
Y* = Y2A(A;ss') (3-8)
A ss


F* F*
Y2B Y2B(AA';ss') (3-9)
AA' ss


yA* = YA(AA'A";ss'S") (3-10)
AA'A" ss's"



YB y3B (FAA'A"; ss's") (3-11)
AA'A" ss's"


In each of these equations A = 1 2 ; in (3-10) A' = 1'', A"= 2

while in (3-9) and (3-11) 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 (3-12)

*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


(3-14)


0
0


(3-15)


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 (3-12) and (3-14) J (p*,T ) is the single integral:

00

Jn(* ,T ) E dr 2 r-2 g(r2) (3-16)
0

Values of the J integrals are tabulated elsewhere and have been fitted
*
to an empirical equation in p and T [37].

In (3-13) LY(1 ;n n ,) is the triple integral:
*
0 r12+r13
LYrV/o *-(n-l) 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) (3-17)



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 (3-15) K (k'k";n n ,n,,) is the triple integral:







*
0c oo r12+ r13
nn *-(n-1) *-(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) (3-18)



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) (3-19)


where C(1V'";mm'm") is a Clebsch-Gordan coefficient [38]. The Ykm are

spherical harmonics in the convention of Rose [38]. In (3-18) and (3-19)

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 (3-14) 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 (3-12), (3-13), (3-14), and (3-15), respectively, for various aniso-

tropic potentials. The results are tabulated in Appendix D.
F
The contributions to y given in (3-12,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 (3-20)
2A 4 J NkT
n nS,

;nn ,) A (AA';ss')
y'F(AA';ss) = 8T 1 snss 2B (3-21)
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 (3-22)
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 (3-21) and (3-23) are defined by:
*
00 oo r12+ r13
*-(n-1) *-(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) (3-24)
*
0o Co r12+ r13
*-(n-1) *-(nt-1)
K(Z'R";nn'n") E dr12 r12 -1) dr3 r13 ,
0' 0 r 12- r13

*-(n"-l) *
r23 23 (r12r13r23) A3-"(ala2a3) (3-25)


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 (3-1)

are presented. For these calculations, the reference fluid surface
F
tension y was obtained from the Fowler model expression for a

Lennard-Jones fluid:










YF 2/ = 3p *2[J5 2J1] (3-26)
0 5 11



where the J integrals are defined by (3-16) 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 (2-92). The points on these figures

are Monte Carlo calculations of the Fowler model surface tension for

the Lennard-Jones 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 Lennard-Jones
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 (




r-4




(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
Lennard-Jones 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 Lennard-Jones
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 (2-96), the temperature derivatives of the anisotropic

contributions to the surface tension are required. Equations (3-12),

(3-13), (3-14), and (3-15) give, respectively:



F* ,in J
y F*(A;ss') F* n +n 1
2A F* s s'-1 1
= Y2A(A;ss'- ) (3-27)
T A T T


F* Y
Y2B(AA';ss') a* n L (;nns) 1
= Y2B (AA';ss') *- (3-28)
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") (3-29)
DT -T T


aY3B(AA'A'; ss 's") Fan KY(Ze''";n n ,n
= Y3B (AA'A";ss's") (3-30)
TT T T


In Table 2 values for u from (2-95) and (2-96) are compared

with computer simulation results for Lennard-Jones 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 Gibbs-Helmholtz Equation in the Fowler Model
Perturbation Theory for Lennard-Jones 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 (3-31)
so 5 11


F
where the J integrals are defined by (3-16). 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 (2-93) in

the Fowler approximation for Lennard-Jones 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 (2-93).


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

quadrupole-quadrupole 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 (n-6)- n 6
(n,6)r\ n 6(n-6) ran 6
u (n6 (r) = n-) (3-32)
(n-6) (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 n-1 about the Lennard-Jones (12,6) fluid free

energy. The surface tension is obtained by applying the thermodynamic

definition (2-1). 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 (3-33)

In (3-33) H12 is the single integral:

H(2,6)( dr* *-(n'-2) (12,6)
H (p,T ) = drn r n g (r ) (3-34)
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 [(1-n2) J J + 6 Hl ]
E n 11 2 5 11

x (- ] (3-35)



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 Lennard-Jones 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 ) (3-36)
o o

and


(12,6) (12,6),rep -Bu(126)rep HS (12,6) (3
g (r) g e y (d (3-37)
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)


(3-38)










In (3-36) and (3-37) the hard sphere diameter d is taken to

be [42]:


d = [1- e-ure(r)]dr (3-39)
0

Further, assuming


yHS(d(n6)) HS (12,6) 3-40)


(3-36), (3-37), and (3-40) give:


(n,6) g(12,6) e-[u (n,6),rep (12,6),rep] (3-41)
g r) (r) e =

(n,6),rep (12,6),rep
where [u(n rep u(126)rep] vanishes for r > r .

Using (3-41) with the Mie (n,6) potential in the integrals Jn,'

Equation (3-16), and K(X'Y";nn'n"), Equation (3-24), Twu has found the

resulting values to be negligibly different from those values obtained

for the Lennard-Jones (12,6) potential [37], at least for values of n

close to 12. Hence, in the calculations reported here, the Lennard-Jones

(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 (3-35) was used to

determine the reference fluid contribution. Values for the Lennard-Jones
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



0-0

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 (3-42)



3.4.2 Anisotropic Contribution to Real Fluid Surface Tension

The calculations presented here include anisotropic dispersion,

overlap, and multiple contributions up to the quadrupole-quadrupole

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) + udis-Q(224)
over-dis over-dis dis-Q

+ u m(112) + u mul(123) + u (213) + u (224)
mul mul mul mul
(3-43)

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 Flytzani-Stephanopoulos 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

(3-8), (3-9), (3-10), and (3-11) with (3-12), (3-13), (3-14), and (3-15),

respectively:

F* F* F* F*
2A = 2A(112) + 2 2A(123) + Y2A(224) (3-44)

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) (3-45)

F* F* F*
Y3B = 3B(112;112;112) + 3 Y3B(112;123;123)
F* F*
+ 3 (123;123;224) + y3B(224;224;224) (3-46)
3B 3B










Expressions for the terms on the right side of (3-44), (3-45), and (3-46)

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 (3-43) have been

included in only the second order term Y2 in calculating the surface

tension from the Pade perturbation theory (2-92). 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

Lennard-Jones 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 Lennard-Jones 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 (3-47)
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 (3-48)


Using (3-48) in (2-87), (2-88), (2-89), and (2-90), together with a

similar approximation for the triplet correlation function, go(r12rl3r23),

(3-47) becomes, in reduced form:


*2 *2 *3 *4
aPL a2PL a3PL a4PL
Y = Y + + + + (3-49)
o *2 *2
T T T T

where



a' E dr r (3-50)
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) (3-51)

00
a d *3 *3 -2
a -2 j dr1 *r1 (3-52)
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) (3-53)
a a a W l 2 3 12 13 23


Equation (3-49) 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 + + + (3-54)
R R

2/3 -1
where R (Vc/NA)2/3 (kT (3-55)



PR = P/Pc (3-56)


TR = T/Tc (3-57)



In transforming (3-49) into (3-54), 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,...) (3-58)

3 *
Pc = Pco = c2(H 'Q ,',K,...) (3-59)
PO
c *
P = = c3(H 'Q ,6,K,..) (3-60)


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 (3-54). 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 (3-54) has been used as a basis for correlating surface

tension by treating the parameters a. as semiempirical constants. Various
1
truncated forms of (3-54), including its Pade form analogous to (2-92),

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) (3-61)



Values for al and a2 for use in (3-61) 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 (3-62)
oR R R


which applies for 0.4 T TR 0.95.

The validity of the correlation (3-61) may be tested by determining

2
whether experimental data give a linear relation between (y YoR)TR/P

and pR' as implied in (3-61). 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 (3-61)
Equation (3-61)


Substance Range a (102) a2(102) References

T P Y
r


Paraffins

Ethane
Propane
n-Butane
i-Butane
n-Pentane
i-Pentane
n-Hexane
n-Heptane
n-Octane
i-Octane
n-Nonane
n-Decane
n-Dodecane
n-Tridecane
n-Tetradecane
n-Pentadecane
n-Hexadecane
n-Heptadecane
n-Octadecane
n-Nonadecane
n-Eicosane

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
1-Butene .48-.70 4.752 -0.9880 49 48
2-Butene .51-.65 3.606 -0.8401 49 48
1-Hexene .54-.64 8.630 -2.142 49 49
1-Octene .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
n-Propanol .55-.68 25.56 -8.854 49 48
i-Propanol .56-.60 33.13 -10.98 49 49
n-Butanol .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 (3-61) 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 (3-61) were usually

within 3% of the experimental values as shown in Figure 17.

Generally, the correlation of (3-61) 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 (3-63)
C.P. C.P.


hence, from (3-61):



a = a2 a (3-64)


Then

2

r R + a (1 p) (3-65)
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 __