摘要:

DISCUSSIONS OF THE FARADAY SOCIETYNO. 43 1967The Structure andProperties of LiquidsTHE FARADAY SOCIETYLONDONDistribution arrangements overleaTHE SOCIETY’S PUBLICATIONSTransactions of the Faraday SocietyDiscussions of the Faraday SocietySymposia of the Faraday SocietyPublished monthlyNormally published twice a yearTo be published once a yearMEMBERSof the Faraday Society receive current issues of bothTransactions and Discussions free on publication, and ofSymposia at a reduced price on request. Enquiriesregarding membership of the Society should be addressed to:The Secretary, The Faraday Society,6 Gray’s Inn Square, London WCl (Telephone: 01-242 8101)NON-MEMBERSmay obtain the Society’s publicationseither through their own bookselleror by making application as follows :Annual Subscriptions:Back Issues:A ustralia:Canada:Discussions :Back Issues:FORto current issues of EITHER Transactions, Discussionsand SymposiumOR Transactions onlyComplete Volumes (comprising Transactions and Discussions)from Vol.41 (1945) onwardsThe Aberdeen University Press LtdFarmers Hall, Aberdeen, ScotlandAPPLY TOAND FORAll current issues and back numbersComplete Volumes (comprising Transactions and Discussions)from Vol. 1 (1905) to Vol. 40 (1944)Butterworth $i Co. (Publishers) Ltd, 88 Kingsway, London WC2APPLY TOOVERSEAS ADDRESSESButterworth & Co. (Australia) Ltd.Sydney: 20 Loftus StreetButterworth & Co. (Canada) Ltd.Toronto: 1367 Danforth Avenue 66New Zealond: Butterworth & Co.(New Zealand) Ltd.South Africa: Butterworth & Co. (South Africa) Ltd.Wellington: 4915 1 Ballance StreetDurban: 33/35 Beach GroveU.S.A.: Discussions-urrent and back issues :Complete Volumes 1-40:Enquiries to 88 Kingsway, London, WC2.Johnson Reprint CorporationNew York: 11 1 Fifth Avenue, A GENERAL DISCUSSIONONThe Structure andProperties of Liquidsllth, 12th and 13th April 1967A GENERAL DISCUSSION on The Structure and Properties of Liquidswas held at The University of Exeter, Exeter, on the llth, 12th and 13th April 1967.The President, Prof. C . E. H. Bawn, C.B.E., Ph.D., F.R.S., was in the Chair duringthe opening and closir,g sessions ; Prof. D. H. Everett, Prof. M. L. McGlashan andProf. J. S.Rowlinson acted as Chairmen during the other sessions. Over 273members and visitors were present, including the following 99 from overseas :Dr. Kurt Altenburg GermanyDr. D. Anderson SwedenProf. and Mrs. F. C. Andrews U.S.A.Dr. V. Ardente Ital~7Dr. G. D’Arrigo ItalyProf. T. A. Bak DenmarkProf. A. Bellemans BelgiumDr. H.-J. Bittrich GermanyMr. C. Bruin NetherlandsProf. G. Caglioti ItalyMr. A. Compagner NetherlandsDr. G. Doge GermanyDr. D. Dorner GermanyDr. M. Dupuis FranceDr. B. Eckstein GermanyProf. H. Eisenberg IsraelDr. T. J. V. Findlay AustraliaMr. J. Fischer AustriaDr. D. D. Fitts U.S.A.Prof. M. Fixman U.S.A.Prof. E. Forslind SwedenProf. E. U. Franck GermanyProf. H. S . Frank U.S.A.Mr. S . Fuks BelgiumDr. and Mrs.J. H. Gisolf NetherlandsDr. J. Chr. Gjaldbaek DenmarkDr. D. Henderson AustraliaDr. I. H. S . Henderson CanadaProf. H. G. Hertz GermanyProf. J. B. Hyne CanadaDr. D. Jonah Sierre LeoneMr. 0. Jsns DenmarkProf. M. L. Josien FranceDr. H. Kehiaian PolandProf. and Mrs. J. Koefoed DenmarkDr. J. J. Kozak BelgiumProf. K.-E. Larsson SwedenMr. J. C. Legros BelgiumDr. D. Levesque FranceDr. A. Levi ItalyDr. J. W. Lorimer CanadaDr. W. A. P. Luck GermanyDr. H. Lutje, GermanyProf. M. Magat FranceMiss B. Maijgren SwederrDr. G. Manfre ItalyDr. P. F. Mijnlieff NetherlandsMr. J. H. Misguicli BelgiumDr. L. Mistura ItalyDr. A. H. Narten U.S.A.Mr. L. E. Oden U.S.A.Dr. J. Padova IsraelProf. R. Pecora U.S.A.Mr. E. J. Picard FranceProf.C . J. Pings U.S.A.Dr. W. Pistor GernmnyMr. T. Plesser GermanyProf. J. A. Pople U.S.A.Dr. E. Praestgaard DenmarkMr. C. Reynolds EireDr. F. P. Ricci ItalyDr. H. Richtering GermanyDr. J.-L. Rivail FranceDr. E. B. Robertson CanadaMr. J. Rouchard CongoDr. R. Satterfield U.S.A.Prof. P. Saumagne FranceDr. T. Schneider SwitzerlandProf. R. L. Scott U.S.A.Dr. K. Skold SwedenDr. F. W. Smith CanadaDr. Care1 C. Smitskamp NetherlandsProf. N. S. Snider CanadaDr. St anton BelgiumDr. J. Stecki PolandDr. D. Stehlik GermanyDr. H. Stiller GermanyMr. D. J. Stufkens NetherlandsDr. G. P. Tliomaes BelgiumMr. J. B. van Tricht NetherlandsMr. R. Tufeu FranceProf. E. Tunkelo FinlandDr. G. Urbain FranceDr. and Mrs. F. Vaslow U.S.A.Prof.L. Verlet FranceMr. and Mrs. Vincent FranceMr. W. Waeber SwitzerlandDr. J. Walkley CanadaProf. A. Weiss GeiwzanyProf. B. Widon U.S.A.Mr. D. A. Young U.S.A0 The Faraday Society and Contributors 1967Printed in Great Britain at the University Press, AberdeeCONTENTSpage 716263240506070758997108115128149160Fourteenth Spiers Memorial Lecture-Structural Theories of Fluidsby G. S. RushbrookeEquation of State of a Monatomic Fluid with 6-12 Potentialby S. A. Rice and D. A. YoungExact and Approximate Values of the Distribution Functions of a SimpleFluid by D. Henderson, S. Kim and L. OdenModijied Cell Model for Liquids: Order Defects and IntermolecularPotentials by F. Kohler and J. FischerCalculation of Thermodynamic Properties of Liquid Argon .from Lennard-Jones Parameters by a Monte Carlo Methodby I.R. McDonald and K. SingerGENERAL DIscussIoN-Dr. J. A. Barker, Dr. D. Henderson, Dr. P.Hutchinson, Prof. J. S. Rowlinson, Dr. S. Levine, Mr. A. Moreton,Prof. G. H. A. Cole, Prof. S. A. Rice, Dr. D. A. Young, Prof. J. Walkley,Dr. W. Y. Ng, Prof. U.V. Weber.Random Close-Packed Hard-Sphere Model.I. Eflect of Introducing Holes by J. D. Bernal and S. V. KingII. Geometry of Random Packing of Hard Spheresby J. D. Bernal and J. L. FinneyHidden Variables in the Critical Region by M. FixmanGENERAL D~scussro~-Mr. G. Mason, Mr. R. H. Beresford, Prof. J.Walkley, Dr. I, M. Hillier, Mr. R. Collins, Prof. D. H. Everett, Dr. J.Finney, Prof.C. Domb, Dr. B. L. SmithStructure of Liquids. Part 7.--Determination of Intermolecular PotentialFunctions and Correlation Functions in Fluid Argon by X-Ray DiflractionTechniques by C. J. PingsX-Ray Diffraction Study of Liquid Water in the Temperature Range4-200°C by A. H. Narten, M. D. Danford and H. A. LevyInfra-Red Absorption of HDO in Water at High Pressures and Temperaturesby E. U. Franck and K. RothSpectroscopic Studies Concerning the Structure and the ThermodynamicBehaviour of H20, CH,OH and C2H50H by W. A. P. LuckGENERAL DIscussIoN.-Dr. M. Davies, Dr. J. N. Sherwood, Dr. G.Caglioti, Dr. M. Corchia, Dr. G. Rizzi, Dr. D. I. Page, Mr. J. W. Perram,Dr. S. Levine, Dr. W. A. P. Luck, Prof. H. S. Frank, Dr. J. Padova,Prof. M. L. Josien, Prof.M. Magat, Dr. B. Eckstein, Prof. J. B. HyneRadiation Scattering Studies of the Structure and Transport Properties ofLiquids by P. A. EgelstaffBrillouin Scattering of Neutrons from Liquidsby B. Dorner, Th. Plesser and H. StillerCONTENTS 6169184192196205212216223231235243248Neutron Scattering Spectroscopy of Liquidsby B. K. Aldred, R. C. Eden and J. W. WhiteGENERAL DiscussroN.-Dr. M. Davies, Prof. K.-E. Larsson, Prof. J.Stecki, Dr. A. LeviNuclear Magnetic Resonance i?i the Study of Liquids by J. A. PopleStudy of Molecular Motion in Liquids by Measurement of Nuclear Relaxa-tion by R. A. Dwek and R. E. RichardsAngular Correlation in Liquids by A. D. BuckinghamFrequency- Dependent Direct Correlation Function by J. SteckiInfluence of Molecular Rotation on Some Physical Properties oj’ Liquidsby D. B. Davies and A. J. MathesonViscoelastic Relaxation in Supercooled Liquidsby A. J. Barlow and J. LambEuclidean Geometry and the Flow of Generalized Liquids by F. W. SmithGENERAL DrscussIoN.-Dr. M. Davies, Prof. A. Bellemans, Prof. H. G.Hertz, Prof. A. D. Buckingham, Dr. G. H. Findenegg, Dr. A. J. Matheson,Dr. R. A. Dwek, Dr. I. Henderson, Dr. F. W. Smith, Dr. B. CleaverSummarizing Remarks by J. S. RowlinsonAUTHOR INDE

ISSN:0366-9033    

DOI:10.1039/DF9674300001

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

FOURTEENTH SPIERS MEMORIAL LECTUREStructural Theories of FluidsBY G. S. RUSHBROOKEDept . of Theoretical Physics, The University, Newcastle-upon-TyneReceioed 13th April, 1967To have the honour of giving the Spiers Memorial Lecture at this General Dis-cussion of the Faraday Society on the Structure and Properties of Liquids is a privilegewhich I prize : but to respond responsibly I must keep to the one aspect of the subjectwith which I have some familiarity, namely, the equilibrium structural theories asthey relate to monatomic fluids. In speaking on these I shall hope to give a convenientbackground to some of the papers with which we shall later be concerned in theearlier half of this meeting.The last Faraday Society Discussion concerned with the structure and propertiesof pure liquids was that held at Edinburgh in 1936.The date coincides, ratherprecisely, with the start of the transition from theories of interpretation to theoriesof prediction as regards the structure of simple fluids. Of course, there is still, andalways will be, room for theories of interpretation : it would be foolish, for example,to attempt to predict the structure of water from first principles. I am glad thatlater we have papers on so important a fluid. But for the really simple fluids,in particular (perhaps exclusively) for the inert gases, the last thirty years have seenincreasingly successful attempts to predict scattering intensities, X-ray or neutron,directly from assumptions about the interatomic forces. And for liquid metalsalso, there is hope that the attempt will be profitable.If we ignore multiple scattering, the scattering function i(s) is essentially theFourier transform of the pair correlation function h(r), or g(r) - I , where g(r) is theradial distribution function.For the customary assumption that the interatomicforces may be represented by additive pair potentials, all thermodynamic propertiescan be derived from this correlation function h(r), if it is known over an appropriaterange of density and temperature. The theory of the structure of simple liquidsis thus intimately enmeshed in the theory of the equation of state and, more par-ticularly, in the theory of phase-changes. This is very satisfactory for those theo-reticians, like myself, to whom the most remarkable thing about the liquid state is thatit exists.The three theories, customarily called Born-Green, hyperchain, and Percus-Yevick,of which any discussion of the theory of fluids must take account, had very differentorigins : though the last two are sufficiently alike in structure to be classed together.It is the Born-Green theory which dates back to the mid-thirties, its two equations,and9(1,2,3) = 9(1,21g(2,31g(1,3), (2)having been given, by Yvon and Kirkwood respectively, in 1935.The first,in which g(1,2) and g(1,2,3) denote the pair and triplet distribution functions and8 STRUCTURAL THEORIES OF FLUIDS4( 1,2) the interatomic potential, is an exact equation : one of a heirarchy of equationslinking successive distribution functions, it is a consequence of Boltzmann’s distribu-tion law.The second is an approximation, the superposition approximation,expressing the additivity of potentials of mean force, which can be argued on physicalgrounds. They were brought together by Yvon in 1937. Yvon also gave (1935)the pressure, or virial, equation,p = pk-T- +p2fr#’(r)g(r)d3r, (3)enabling us to pass from the distribution function g(r) to the equation of state.But the time was not yet ripe for the exploitation of these equations. Some tenyears later (1946) the theory was redeveloped, independently of Yvon’s work, bothby Bogolyubov and by Born and Green.5 Within four years (1950), Kirkwood,Maun and Alder had solved the equations, and those of the very similar theorydeveloped by Kirkwood himself, for hard spheres, on an I.B.M.computer : and withthis we enter the modern era. For without machine calculations, and machineexperiments, no progress in the structural theory of fluids is possible.Although quantitatively, the Born-Green theory is very much inferior to itslater rivals, perhaps no theory has yet had greater consequences. For it was thediscovery by Kirkwood, Maun and Alder that the equations did not have integrablesolutions for densities greater than pa3 = 0.95, where p is the number density andc the sphere diameter,* which led Alder and Wainwright to embark on machinecalculations, by the method now known as molecular dynamics, to produce an experi-mental equation of state for a hard sphere fluid.The results, first published in1957, and supported by independent Monte Carlo calculations by Wood andJacobson,* gave rather clear evidence of a transition from a fluid to a solid phase ata density approximately that predicted by the Born-Green theory. The Born-Greentheory gives appreciably toG low a pressure, perhaps 30 % too low at the limit pa3 =0.95 ; but the later theories, while they do justice to the experimental fluid isotherm,would not have predicted the fluid-solid transition.One may legitimately question the reliability of these machine calculations,inevitably based on the dynamics of a comparatively small number of particles. Butour confidence is restored by the close agreement between the fluid branch of thisexperimental equation of state and that calculated by the quite independent methodof evaluating successive virial coefficients.It is a sobering thought that whereasB~ltzmann,~ in 1899, using some calculations of von Laar, could writeP+ ap2 = k-Tp[l+ bp + 0-625b2p2 + 0.2869b3p3 + . . .(where b = 2na3/3), and regard it as a more precise form of van der Waals’s equation,even now, with the indispensable advantages of the Ursell-Mayer theory and high-speed computers, we can add only three more terms lo+0.1103b4p4+ 0.0386b5p5 +0*0138b6(with some uncertainty in the last digits). Nevertheless, extrapolation of the extendedseries by the method of Pad6 approximants, gives nice agreement with the fluidbranch of the experimental isotherm : which thus constitutes a firm, if artificial,base line for testing theories at high temperatures.And probably few of us todaywould disbelieve that freezing, or melting, is essentially an excluded volume problem :as was indeed foreshadowed by the considerations of Bernal at the Edinburgh Con-ference in 1936, in one of the earliest papers of the predictive type.* In these units, close packing corresponds to p03 = 1.414G . S . RUSHBROOKE 9But we must turn to the other theories. Perhaps their most attractive featureis that both of them incorporate the Omstein-Zernike l 1 equationh(1,2) = ~(1,2)+p[c(l,3)h(2,3)d3, (4)relating the total correlation function h(r), or g(r)- 1, to the direct correlation functionc(r). Essentially this equation provides the definition of one function in terms ofanother, and nothing more.It has the consequence, of course, that the compressi-bility equation of Zernike and Prins (1927) l2.Taplap = i+pJh(1,2)d2 = qo), (5)where I(0) is the structure factor, 1 + i(s), for limitingly small scattering angles, andwhich is a second route to thermodynamic behaviour, can be written in the alternativeform,-!- 2 = 1-p c(1,2)d2KT ap s This shows that for a highly compressible fluid, near its critical point, c(r) is short-ranged compared with h(r). But whether we have gained anything, whether c(r) isa useful physical concept, can be judged only by whether c(r) is more simply, ortransparently, related to the interatomic potential than is h(r). A definition is nota theory.The simplest possible theory, of hyperchain or Percus-Yevick type, is to assumethat c(r) is simply the Mayerf-function, i.e., to writeThis gives the direct correlation function the range of the interatomic potential, andincorporates a temperature dependence of the correct Boltzmann form.The assump-tion is intuitively appealing : and both Ornstein and Debye, in their work on criticalopalescence, argue very much along these lines. In the theory of liquids, however,it is not good enough. It leads,13 essentially, to the crudest form of the linearizedBorn-Green equations. Nevertheless, it corresponds precisely to what Montroll andMayer l4 did in the theory of ionic solutions, when they derived the Debye-Huckeltheory by summing over a certain class of interaction-graphs (in the Ursell-Mayerformalism) : namely the chain, or ring, diagrams. And it was the achievement ofMontroll and Mayer in showing that physically interesting theories could result fromsumming over subsets of interaction-graphs in a low density expansion which ultimatelyled to the hyperchain equation : though we must not lay the defects of this theoryat their door.The set of graphs over which we sum to obtain the hyperchain equation is dictatedfaute de mieux, rather than on physical grounds.It is just the largest subset, completewithin its topological class, which we can handle conveniently. But the result,obtained independently by many people some eight years ago,15 is a theory whichretains the Ornstein-Zernike equation together with a second equation for c(r),namely,c(r) = exp [ - #(r)/lcT] - 1.(7)c(r) = h(r) -In [ 1 + h(r)] - 4 ( r ) / ~ T . (8)The Percus-Yevick theory,16 although of the same structure, had an entirelydifferent origin. It could have arisen, as Stell later showed, out of the graphicalanalysis that led to the hyperchain theory : retaining a different, but equally plausible,set of interaction graphs. In fact, it arose out of an attempt to define collectiv10 STRUCTURAL THEORIES OF FLUIDScoordinates for a fluid, in a paper (1959) which I am probably not alone in findingobscure. Its best derivation is undoubtedly that which Percus later gave (1962)using the mathematical technique of functional differentiation recently introduced,in the theory of liquids, by Yvon.” This, which examines the change in densityat one point produced by a change in force field at another, brings us closest tothinking physically, and furthest both from the dilute gas and from arbitrary assump-tions of mathematical convenience. Fortunately, as Percus showed, all our theoriescan in fact be derived this way: though the Born-Green theory then appears assomething of a tour deforce.But however it is derived, the Percus-Yevick theoryjoins to the Ornstein-Zernike equation the second prescription for c(r),c(r) = [ 1 - exp ( 4 ( r ) / ~ T ) ] [ 1 + h(r)]. (9)We must now turn to the testing of these theories. They can be used in threerather different ways. Assuming a knowledge of 4(r), we can calculate the thermo-dynamic properties : or we can calculate the correlation fiinctions and structurefactors.Alternatively, we can use the theories to determine 4(r) from scatteringexperiments.Not unnaturally, the earliest calculations, with any theory, are the simplest :in this case calculating the predicted values of successive virial coefficients for hardspheres. These are not exact theories, nor are they entirely internally consistent.Different routes to any thermodynamic property will lead to different answers.The two routes most commonly used are based, respectively, on the virial equationand the compressibility equation. Ultimately we may have to be content with themost satisfactory prescription (preferably understanding why it succeeds) ; buta serious lack of internal consistency is certainly a defect in any theory.At first, for hard spheres, the test was mainly one of internal consistency, sincethe true virial coefficients were not known beyond the fourth.But now that we knowthe rigid sphere isotherm with some confidence, it is more satisfactory to comparethe theories directly with it rather than with low density expansions.18 Judged inthis way, the Percus-Yevick theory is incomparably the best: in respect of bothinternal consistency and agreement with the true isotherm. There is the addedadvantage, as Thiele l9 and Wertheim 2o have shown, that, for hard spheres, thePercus-Yevick equations can be solved exactly, in closed form. The compressibilityprescription accords with the known isotherm to within perhaps 2 % over the wholefluid range.But the thermodynamic properties show no evidence of a phase transi-tion : indeed the Thiele-Wertheim solutions show no singularity until we reach thequite unphysical density pa3 = 1.9, far greater than that of close packing (pa3 = 1.4).The Thiele-Wertheim solution, however, enables us also to find the correlation func-tions ; and we now know 21 that for pc3 >1-18, g(r) takes on negative values. Sincethis is physically impossible, it gives us an excuse for disregarding the more implausiblepredictions of the Percus-Yevick theory.A hard sphere gas, however, though it provides a useful high-temperature baseline, is far removed from a real liquid. To separate the fluid phase, at low enoughtemperatures, into gas and liquid regions, we must add an attractive potential:and we should soften the repulsive core.Quite appropriately, most further workhas been based on the Lennard-Jones potential4(r) = 4&[(a/r)12 - (a/r)7.For this potential, all three theories have been compared with a high-temperatureMonte Carlo isotherm, for T*(= KT/E) = 2-74, obtained ten years ago by Wood andParker 2 2 : there is quantitative agreement, or disagreement, of the same order as iG . S . RUSHBROOKE 11the hard sphere case. But this is still at twice the critical temperature ; and we mustnext turn to the critical constants themselves.These critical constants have now been found for all three theories, notably byL e v e s q ~ e , ~ ~ whose recent paper is the source of the numbers I shall quote.For thehyperchain theory, hitherto the most studied, Levesque’s estimates agree, withintheir limits of error, with those of Klein and Green,24 and de Boer et ~ 1 . ~ ~ Table 1shows T:, pccr3, and the ratio (p/prcT),, for the Lennard-Jones potential,TABLE 1B.G.h.c.P.Y.“ argon ”T:1 -451 -251.251-26PcQ30.400.260.290.32( P I P K T ) ~0-440.350.300.30together with values customarily quoted for argon, assuming the Lennard-Jonesparameters Elk = 120°, cr = 3.4A. Bearing in mind that Tz, for the theoreticalpredictions, is uncertain to perhaps 2 x, and that the uncertainty in p c is greater,we might well conclude that the modern theories were by no means inadequate.But this is not the whole story : first, because I have deliberately chosen the valuesobtained from the virial equation, and not from the compressibility equation ;secondly, because we do not really know what the true experimental values are fora Lennard-Jones gas.To take the second point first, Levesque has also determined the critical constantspredicted by the hyperchain and Percus-Yevick theories for the hard core with squarewell potential, +(r) = - E , o<r<1*5a, and the corresponding results are shown intable 2.This table shows also the predictions of Verlet’s P.Y.11 equation 26 (aTABLE 2h.c.P.Y.P.Y. I1M.D.Tr1.161.161 -281.28Pc030.250-2 10.320.33(PIPKTIc0.320.300-320.3 1natural extension of the Percus-Yevick theory, when approached through functionaldifferentiation), together with recent “ experimental ” results, for this potential,found by Alder from molecular dynamic studies (and quoted by Levesque).It nowlooks as if both hyperchain and Percus-Yevick theories do only scant justice to thecritical region, yielding T‘ to within perhaps 10 %. The P.Y.11 equation, however,seems to be appreciably more accurate.I shall not discuss the corresponding figures obtained by use of the compressibilityequation, save to say that the internal inconsistency in the theories is of the sameorder as their inaccuracy. It is already clear both that first-order theories give onlya moderately satisfactory account of thermodynamic properties in the critical region,and that there is need for caution before drawing conclusions for real fluids.At greater densities and, more particularly, at lower temperatures, in the trueliquid region not far from the triple point, the thermodynamic situation is muchworse: certainly if the virial equation is used to find the pressure.To take anextreme example, Gaskell 27 has shown that were we to use the hyperchain equation12 STRUCTURAL THEORIES OF FLUIDSto infer the pressure from accurate scattering experiments on liquid argon, essentiallyby solving (8) for 4(r) and substituting in (3), we should overestimate the pressure atthe triple point by a factor of about lo3. But it does not follow from this that thetheories are useless. We are far from being able to calculate a vapour pressure curvefrom first principles : but we may yet have a good theory of liquid structure.Thevirial equation is exceptionally sensitive to the precise form of g(r) in the liquid region.If we look at the internal energy the results are much more encouraging ; indeed, asLevesque shows, hyperchain and Percus-Yevick predictions are in very close agree-ment with Monte-Carlo results for the Lennard-Jones potential. And there is theobservation by Ashcroft and Lekner,28 that if the neutron diffraction data on alkalimetals are fitted, as regards the first peak, by the first peak of the Percus-Yevickhard sphere theory, which for Na at 100°C means choosing the density value po3 =0.85, then the compressibility from zero-angle scattering on Percus-Yevick theory isextraordinarily close to the experimental, thermodynamic, value (within 5 %).Whilethis may tell us more about metals than about our theories, it at least suggests thatthe Percus-Yevick theory does good justice to the excluded volume problem in theliquid state, even at effectively high densities.Perhaps our greatest present need is to be able to answer for the structural pro-perties the questions we have already largely answered for the thermodynamic ones,namely, over what regions of the fluid phase, gas, liquid, or supercritical, equationssuch as (8) and (9) adequately interrelate the structural properties, h(r) and c(r),or their Fourier transforms, and the interatomic potential, #(r). We may ask, ofcourse, adequately for what ? Here, surely, the most physically important questionis whether current theories enable us to infer properties of the potential 4(r) fromthe best experimental measurements of the structure factors of real fluids.I am notthinking of the oscillatory potentials found by Johnson, Hutchinson and March 29in this way from the structure factors of liquid metals, oscillations which have beenonly partially confirmed by the later work of Ascarelli 30 : though it is to the creditof the theories if they can point to significant differences between metals and inertgases. The issue, rather, is whether we can attach quantitative significance to suchpredictions of $ ( r ) for the inert gases themselves.This question can be answered only by computer studies of the kind on whichVerlet 31 has recently embarked. For a chosen 4(r), h(r) and c(r) must be determinedby machine experiments over a wide range of p and T : then equations such as (8)and (9) can be tested by seeing to what extent, and over what range of p and T,they successfully reproduce 4(r).A preliminary report from Verlet is encouraging,for although at large densities or low temperatures the shape of $(r) is not wellreproduced, nevertheless for temperatures of the order of T,, and for a range ofdensity extending well beyond pc, the depth of the potential minimum in +(r) is re-produced by (9) to within 2 %. The analysis, by Mikolaj and Pings,32 of theirexperimental measurements on argon, suggests that the hyperchain equation, (8),may have an approximately equal validity in this region.More importantly, if theseresults can be trusted, the finding by Mikolaj and Pings that the depth of the potentialminimum so deduced for argon varies linearly with p must afford evidence of theinadequacy of the assumption of additive pair potentials, and demonstrate that forreal liquids we are concerned also with multibody forces. If, and when, our theoriesof dense fluids enable us to say this with complete confidence, and with some quantita-tive precision, they will have been justified.Attempts to improve on the present theories usually bring us back to the problemof the triplet distribution function, g(1,2,3). This is inevitably true of the Born-Green theory, where it must be the superposition approximation that is inadequatG. S. RUSHBROOKE 13since the other, Yvon, equation is exact : and it is the very essence of Born-Greentheory to use this exact equation.Prof. Rice will later be telling of his work withLekner to replace eqn. (2) by a more adequate approximation. But it is true also ofthe hyperchain and Perms-Yevick theories, which successfully by-pass the tripletdistribution problem in their original forms. If we proceed systematically along thepath indicated by the functional differentiation approach of Yvon, Percus and Verlet,the Ornstein-Zernike equation is retained, but the correction to the second expressionfor c(r) involves g(1,2,3) : and the theory is not self-contained without an additionalassumption. Though a crude approximation may suffice for calculating a smallcorrection, and good though Verlet’s P.Y.11 theory appears to be, we do not knowhow sensitive the predictions are to the element of arbitrariness here, and may doubtif these second-order theories are yet in their final form.We are also faced with the triplet distribution function if we attempt to bringthree-body forces, or triplet potentials, into the theory of fluids.For although thecompressibility equation ( 5 ) or (6), is unmodified, the pressure equation (3), mustnow include the virial of the three-body forces: which implies a term involvingg( 1,2,3). Nevertheless, if we are content with predicting the scattering functions?triplet potentials can be introduced at the level of the hyperchain and Percus-Yevicktheories : though it is less certain that the theories are quite strong enough to standthe strain.On the hyperchain theory, the inclusion of triplet potentials is particularly simple.33For eqn. (8) is just a truncated form of the exact equationh(1,2) - c( 1,2) = In [ 1 + h(1,2)] + 4( 1,2)/rcT- x( 1,2) (10)where, in the graphical analysis? x(1,2) corresponds to the so-called “ elementary ”graphs : which are neglected in the hyperchain theory. With only two-body forces,the first elementary graph has two field-points, and is proportional to a term in p2.With triplet potentials? however, there is a new class of elementary graph havingonly one field-point? and summing to giveM , 2 ) = P [ f ( 1,2,3)e(ly3)e(2,3)d3where e(i, j) is the Boltzmann factor for the pair potential, exp [---4(i, j)/rcT] andf(i,j,k)is the Mayer function for the triplet potential, exp [ - 4(i,j,k)/rcT] - 1.Thus, at thelevel of hyperchain theory, the use of (8) to infer +(r) from h(r) and c(r) will lead notto the true pair potential 4(1,2) but to an effective pair potential4*(1,2) = 4(1?2) - KT.l(1?2),I would expect this result to have the same range of validity as the hyperchain theory :but shall not pursue the matter further since Prof. Pings will be describing an alterna-tive, and possibly sounder, approach to what I think are essentially the same con-clusions.But even if we refine our theories, by going to a higher order of approximationor including triplet potentials, there remain two outstanding problems : the problemof asymptotic behaviour and the problem of freezing.It is ironic that whereasOrnstein and Zernike introduced the direct correlation function c(r) to argue intui-tively that this is short-ranged compared with h(r), and consequently that for large rh(r) - exp ( - ccr)/v,where a+O at the critical point (or wherever the compressibility becomes infinite1 4 STRUCTURAL THEORIES OF FLUIDSit is just this classic conclusion that is today most open to attack.34 All our approxi-mate theories suggest that c(r) decays like 4(r), (more accurately, - $(r)/rcT), whichfor real fluids, in the absence of triplet potentials, means like l / r 6 : i.e., more slowlythan does h(r) on the Ornstein-Zernike theory. Even if we try to salvage, and refine,the Ornstein-Zernike result by supposing that perhapsh(r)-A&)+B exp (-ar) cos @+y)/r,(where the coefficients are functions of p and T, and a and p vanish at the criticalpoint), we know that the two-dimensional analogue of this is disproved by the workof Onsager and Kaufman on the Ising One technical comment may bepermitted.SinceI($) = 1 +i(s) = 1 + P i ( $ ) = 1/[1 -pE(s)] (1 1)[where &) is the Fourier transform of h(r), and in scattering experimentss = (4n/R) sin (0/2)] it is clear that 1 -pe'(s) cannot become negative for real s :we cannot have a negative intensity of scattered radiation. On the Ornstein-Zerniketheory, for p = pc and T = Tc the equation, 1 -pi.(s) = 0, has a double root at theorigin. It would seem that any heuristic theory must consider the roots of 1 -pE(s) = 0in the complex plane, without assuming that ?(s) is an analytic function of s2 nearthe origin : since this last assumption leads to the Ornstein-Zernike result.Fisher'shypothesis, that (in three dimensions)h(r)-exp (-ar)/rl+q, where O<q<lavoids this assumption : and its two-dimensional analogue can be made to cover theknown behaviour of the Ising model. Nevertheless, for van der Waals forces, itdoes not avoid the conceptual difficulty of a correlation function which decays morerapidly than do the interatomic forces whxh produce it.36 While this may not beinherently impossible, certainly the theory is not yet so firmly rooted as to force usto so surprising a conclusion.Writing eqn. (5) and (6) aswe see that near freezing, indeed anywhere in the true liquid region, when the fluid iscomparatively inconipressible, 1 - pZ(0) is large (though we are not now concernedwith infinities) and it is 1 +pz(O) which is small.The question whether our theoriesof liquids give us any real indication of the onset of freezing, or are capable in principleof so doing, is, I think, the most important with which we are at present faced. Is itaccidental, or significant, that the Born-Green theory suggested a phase-change fora hard-sphere gas? It was, of course, l(0) which ceased to exist : though withoutthe divergence which would indicate infinite compressibility. If freezing is essentiallyan excluded volume problem, can we approach a theory of freezing against a con-tinuum background, or must we discuss the excluded volume problem against thebackground of a lattice-gas ? Are there other, physically acceptable, solutions ofthe Percus-Yevick equations for hard besides that found by Wertheimand Thiele? I think we do not know the answers to these questions : and althoughthe non-equilibrium, or dynamic properties of fluids offer a richer experimental field,there is certainly still work to be done in the development of a satisfactory equilibriumtheoryG. S .RUSHBROOKE 151 J. Yvon, Actualites Sci. Ind., (Hermann et Cie, Paris, 1935), p. 203.2 J. G. Kirkwood, J. Chem. Physics, 1935, 3, 300.3 J. Yvon, Actualites Sci. Ind. (Hermann et Cie, Paris, 1937), 542, 543.4N. N. Bogolyubov, J. Phys. U.R.S.S., 1946,10, 257, 265.5 M. Born and M.S. Green, A General Kinetic Theory of Liquids (Cambridge University Press,6 J. G. Kirkwood, E. K. Maun and B. J. Alder, J. Chem. Physics, 1950,18,1040.7 T. Wainwright and B. J. Alder, Suppl. Nuovo Cimento, 1958,9, 116.B. J. Alder and T. Wainwright, J . Chem. Physics, 1957, 27, 1209 ; 1959, 31, 459.8 W. W. Wood and J. D. Jacobson, J. Chem. Physics, 1957,27, 1207.W. W. Wood, F. R. Parker and J. D. Jacobson, Suppl. Nuovo Cimento, 1958,9, 133.9 L. Boltzmann, Proc. Roy. Acad. Sci. Amst., 1899, 1, 398.10 F. N. Ree and W. G. Hoover, J. Chem. Physics, 1964, 40, 939 and private communication.11 L. S. Ornstein and F. Zernike, Proc. Acad. Sci. Amst., 1914, 17, 793.F. Zernike, Proc. Acad. Sci. Amst., 1916, 19, 1520.12 F. Zernike and J. A. Prins, 2. Physik, 1927, 41, 184.13 G. S. Rushbrooke and H. I. Scoins, Proc. Roy. SQC. A , 1953,216, 203.14 E. W. Montroll and J. E. Mayer, J. Chern. Physics, 1941, 9, 626.15 M. S. Green, J. Chem. Physics, 1960, 33, 1403.1949).J. M. J. van Leeuwen, J. Groeneveld and J. de Boer, Physica, 1959, 25, 792.E. Meeron, J. Math. Physics, 1960, 1, 192.T. Morita and K. Hiroike, Progr, Theor. Physics, 1960, 23, 1003.G. S. Rushbrooke, Physica, 1960,26,259.L. Verlet, Nuovo Cimento, 1960, 18, 77.G. Stell, Physica, 1963, 29, 517.J. K. Percus, Physic. Rev. Letters, 1962, 8, 462.16 J. K. Percus and G. J. Yevick, Physic. Rev., 1958, 110, 1.17 J. Yvon, Suppl. Nuovo Cimento, 1958, 9, 144.18 J. S. Rowlinson, Report Prog. Physics, 1965, 28, 169.19 E. Thiele, J. Chem. Physics, 1963, 39, 474.20 M. S. Wertheim, J. Math. Physics, 1964, 5, 643.21 G. Throop and R. J. Bearman, J. Chem. Physics, 1965,42, 2408.22 W. W. Wood and F. R. Parker, J. Chem. Physics, 1957, 27, 720.23 D. Levesque, Physica, 1966,32, 1985.24 M. Klein and M. S. Green, J. Chem. Physics, 1963, 39, 1367.25 J. de Boer, J. M. J. van Leeuwen and J. Groeneveld, Physica, 1964, 30, 2265.26 L. Verlet, Physica, 1965, 31, 959.27 T. Gaskell, Proc. Physic. SOC., 1966, 89, 231.28 N. W. Ashcroft and J. Lekner, Physic. Rev., 1966, 175, 83.29 M. D. Johnson, P. Hutchinson and N. H. March, Proc. Roy. Soc. A , 1964, 282, 283.30 P. Ascarelli, Physic. Rev., 1966, 143, 36.31 L. Verlet, Preprint, 1966.32 P. G. Mikolaj and C. J. Pings, Physic. Rev. Letters, 1965, 15, 849.33 G. S. Rushbrooke and M. Silbert, Mul. Physics, 1967, in press.34 M. E. Fisher, J. JAath. Physics, 1964, 5, 944.35 B. Kaufman and L. Onsager, Physic. Rev., 1949, 76, 1244.36 B. Widom, J. Chem. Physics, 1964, 41,74.37 H. N. V. Temperley, Proc. Physic. Soc., 1964, 84, 339.R. J. Baxter, Physic. Rev., 1967, 154, 170

ISSN:0366-9033    

DOI:10.1039/DF9674300007

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Equation of State of a Monatomic Fluid with 6-12 PotentialBY STUART A. RICE AND DAVID A. YOUNG*Institute for the Study of Metals and Dept. of Chemistry,The University of Chicago, Chicago, Illinois 60637Received 24th October, 1967The modified Yvon-Born-Green equation of Rice and Lekner is solved for a monatomic fluidwith 6-12 potential at kT/c = 2.74. The equation of state and excess internal energy are calculatedand compared with Monte Carlo calculations ; the agreement is good. Calculations with hard-spheremodifications and with the Kihara potential are also discussed. The Kirkwood limit of stability forthe fluid phase is calculated.A study,’ by Rice and Lekner, of the equation of state of a rigid-sphere fluid, hasshown that the Yvon-Born-Green equation may be so improved as to yield anequation of state in close agreement with the molecular dynamics data of Alder andWainwright.2 In this paper we present an analogous study of a fluid with moleculeswhich interact through a 6-12 Lennard-Jones potential.Monte Carlo calculations 3of the equation of state corresponding to this potential have been reported for asupercritical isotherm at kT/E = 2-74, where E is the depth of the potential well.Our calculations, based on the improved YBG equation, will be at this temperatureso as to permit comparison of theory and “experiment”.Recent work 4 has provided the solutions of the Yvon-Born-Green (YBG),Percus-Yevick (PY), and the Convolution Hypernetted Chain (CHNC) equationsfor the 6-12 potential at kT/E = 2.74. As for the hard sphere potential, the PYtheory yields the best agreement with experiment over a wide density range, but givesno explicit evidence for a liquid-solid phase transition.Our motive for choosing toextend the YBG theory is that it predicts a transition to occur at high density, whilethe PY and CHNC theories do n0t.t Hence we regard the YBG equation to containthe qualitative features necessary to describe the limits of stability of the fluid phase.TRIPLET DISTRIBUTION FUNCTIONThe exact, formal, expression for the triplet distribution function of a fluid isg(3)(123) = g(2)(12)g(2)(13)g(2)(23) exp ~(123, p),~(123, P ) = C pngfl+3(123).(1)(2)wherecofl=lHere p is the fluid number density. The 6, are expansion terms developed and dis-cussed by Meeron 5 and Salpeter.6 The superposition approximation may be obtainedby placing z = 0.This approximation yields the correct third virial coefficient fora fluid, but incorrect higher virial coefficients. Inclusion of terms up to 6, yieldscorrect virial coefficients up to Bm.*Fannie and John Hertz Fellow.?Note added in proof: The PY equation admits solid-like solutions above a critical density butdoes not prohibit simultaneous fluid-like solutions.1S . A . RICE AND D . A . YOUNG 17For the hard sphere potential, Rice and Lekner (referred to as RL henceforth)found that 6 4 and 8 5 could be calculated with a reasonably small amount of computertime, but that 8 6 could not. For the 6-12 potential we find that 84 can be calculatedaccurately, that 6s can be estimated fairly well, and that 156 cannot be calculated at all,again owing to the computer time involved.Starting with 64 and 65, we shall usethe technique of the Pad6 approximant to achieve an estimate of the entire 6, series.The derivation of the diagrams corresponding to the 6 is fully explained in RL andhere we only point out those aspects of the calculation which are special to the 6-12potential.The calculation of 64(123) was carried out numerically by means of Gaussianquadrature. Now0'= p f 1 4 f 2 4 f 3 4 ,I6,( 123) = a4/ \(3)2 0 ' 0 3where hj = exp( - u(ij)/kT) - 1 and d4 is the element of volume swept out by particlenumber 4. There are no boundary problems such as arise in the hard sphere casebecause the 6-12 potential is continuous.Since for moderately large distancesr,j&(r) is small, the integration in (3) may be carried out over a small region. Takingadvantage of the symmetry of the plane of the triplet (123) with respect to the per-pendicular z axis, we find that a hemisphere of radius 3.0 molecular diameters servesthe purpose adequately. All distances will be given in units of the molecular diameter,the diameter being defined as the distance CT for which the potential u(a) = 0.The integration (3) was divided into three parts. In the interior region(O-O<r<l-O), where the integrand is large, the integral was evaluated by using Carte-sian coordinates and 16 point Gaussian quadrature. Then the contributions fromthe two spherical shells (1.0 GrG1-5) and (1.5 <r G3.0) were evaluated in sphericalpolar coordinates using 12 point Gaussian quadrature.The results comparefavourably with those obtained by using larger numbers of Gaussian points, viz., 20,16, and 16. The average error is about 1 %. The integrations required in (3) werecarried out for all triplet configurations (R,S,T) such that 0.0 <R,S,T<2*6 andAR = A S = AT = 0.20. R = 2.6 was found to be satisfactory cutoff since64(2.6,2.6,2.6) = 0.0015, which is close enough to zero to be negligible.The 65( 123) diagrams were grouped and evaluated using the same techniques as inRL. There are 13 diagrams to evaluate. By using the functionspij = d5fi5 fj5; tijk = S,(ijk); eij = fij+ 1 = exp [I -u(ij)/kT],the calculation may be reduced to three 3-fold integrations and one 6-fold integration ;-5.J d4f 1 4 f 2 4 f 3 4 I d5f 1 5 f 2 5 f 3 5 f 4 5 * (4)The first three integrals may be evaluated by using values of 64 tabulated with aninterval of 0.1. Since the storage of 64 in this state in the computer required approxi-mately 19,700 locations, it was not feasible to increase the state of subdivision of 84.Hence the tijk required in the integrand had to be truncated to the nearest tabulatedvalue, at the expense of accuracy.This difficulty was avoided in the work with hard spheres because 84 is known interms of elementary functions and can be rapidly computed for any desired tripletconfiguration. Exactly how much accuracy is lost by the truncation is not known18 EQUATiON OF STATE OF A MONATOMIC FLUIDHowever, very accurate values of pij were calculated and these were often much largerthan the tijk which appeared with them in the integrands. Hence the inaccuraciesin tijk are masked to some extent by this discrepancy in size.Again the individual integrals were divided into three parts, this time using 12, 12,and 8-point Gaussian quadrature in a sphere of radius 3.0.Comparison with moreaccurate calculations shows an average error of 5 % due to imprecise integration.The total error is not known because of the truncation of the ?i& terms. Once again65(R,S,T) was obtained for all configurations up to R = S = T = 2-6 andAR = A S = AT = 0.2.64 and 8 5 are thus obtained as three-dimensional arrays with a spacing of 0-2 diam.units.These values were interpolated by using a linear Taylor expansion. In theexamples below, R, S, and Tare the tabulated triplet lengths (multiples of 0.2) and x, y ,.3.2. I0 -. I-.2-.3-.4-. 5-. 6-.7FIG. 1 . 4 4 and 85 for the 6-12 potential at ~ T / E = 2.74 in the equilateral configuration (R,R,R).and z are the interpolated lengths (odd multiples of 0.1).in this way. There are three cases to be considered :Both d4 and 8s are tabulatedI.11.111.G(R,S,z) = .5(6[R, S, ( Z + el)] + S[R, S, ( Z - .1)]>,G(R,y,z) = .5(6[R, ( y + *1), ( z - .1)] + 6[R, (y - .l), (2 + el)]},G(x,y,z) = .5(6[(x + -1), (y - *l), ( z - d)] + 6[(x - -1), (y + *1), ( z - .1)] +6[(x - -1), (y - .1), (2 + -1)] - 8[(x - .1), ( y - .1), ( z - ml)]}.( 5 4(5b)(54In each case the 6 in the parentheses on the right side are tabulated values, or they arenonexistent owing to the violation of the restriction that R, S, and Tmust be the sidesof a triangle.In the event that one of the 6 necessary for the interpolation does notexist, a simpler approximation is used. The result is some 2700 pieces of data foreach of 6 4 and 8s.Several special values of 64 were calculated and checked against interpolatedresults in order to test the interpolation scheme. The error ranged from 3 to 20 %(for case 111), with an overall average error of 5 %. Hence the overall picture of thetabulated 84 and 6s values is one of fair but not high accuracy.Both 84 and 8 5 take on positive and negative values, and both are large andnegative near (O,O,O).In this sense they differ strongly from the hard sphere terms inwhich 64<0 and 85>0 over their respective ranges. 8 4 and 8 5 for the equilateraltriplet configuration are shown in fig. 1S. A . RICE AND D . A . YOUNG 19YVON-BORN-GREEN-EQUATlONThe YBG equation has been reduced to a form which can be used for numericalcomputation. Details of the reduction are described by Hill.7 We write, followingthe same arguments,kT In g(R) = - zl(R) + npIW .'(y)g(y)dy"y y 2 -(R-X)2[g(x) - Ilxdx +0 R-yu(R) = 4&( - A).Here R, x, y and z are reduced distances with units of molecular diameters. Eqn. (6)is exact, but since only the first two terms of z(x,y,z) are known, we must resort toapproximation if a solution is to be found.A good approximation to the entireseries is the Pad6 approximantwhich we henceforth adopt.NUMERICAL INTEGRATION OF THE MODIFIED YBG EQUATIONWe have obtained numerical solutions over a range of densities for(i) z = 0 (superposition) ; (ii) z = p64 ;(iii) z = pS4+p265 ; (iv) z = p64/(1 -p&/64) = z(Pad6).(v) z = p64(HS) ; (vi) z = p64(HS) +p265(HS) ;(vii) z = Pad6 (HS) ;Also we have tried mixing terms for hard spheres and the 6-12 potential:(viii) z = p64(6-12)+p265(HS).Of the terms (v)-(viii) only (v) turns out to be of real interest.In each of the above calculations the parameter p was selected in the range0.1 < p G1.2, and kT/& was fixed at 2.74. Solutions of the YBG equation wereobtained by the iteration procedure commonly used in solving nonlinear integralequations.This is discussed in RL. The input for the (k + 1)th iteration is taken to beFor low densities, we found that the optimum CI can be as high as 0.75, and that forhigh densities CL must be reduced to 0.2, where convergence is correspondingly slower.If a is too large for a given value of p the iterative solution undergoes a divergentoscillation. For a given z and p the first input used was the convergent solution of thepreviously chosen density in the series of calculations.Convergence was observed by calculating the reduced pressure for each iteration :g r ( r ) = d n ( r ) + aiX"t(r) - gt(r)l- (8)u'(r)g(r)r4nr2dr. (9)The approach to a stable solution could be either monotonic or oscillatory.Thelarger the value of a the greater the tendency for oscillatory behaviour. Convergenc20 EQUATION OF STATE OF A MONATOMIC FLUIDwas usually taken to mean that the input and output g(r) differed by not more than0.5 % at the highest density. The accuracy of the solutions depends on : (i) theinterval Ar used in the numerical integration ; (ii) the maximum value of R to whichthe integration is carried ; (iii) the number of iterations performed ; (iv) the accuracyFor the superposition approximation, some solutions were obtained for Ar = 0.05.The pressures calculated from these results differed from results interpolated fromAr = 0.1 by about 3 % for p = 0.8. The superposition pressures at high densitycoincide with those of Broyles 4 to within 2-3 %.Hence it may be assumed that theresults obtained for Ar = 0.1 are reasonably accurate.Interpolations of g ( r ) were quadratic (i.e. second order in the Taylor expansion)and the pressures obtained from interpolation (Ar = 0.05) differ from those of theoriginal results (Ar = 0-1) by about + 5 %. Most of the contribution to the pressureoccurs near r = 1.0 where the function g(r) is steep. A special interpolation techniquewas used to achieve an increment of 0.05 in the integration of the triple integral ineqn. (6). For z = p&+pZ&, for densities above p = 0.6, divergent results wereobtained with the increment 0.1. These divergences disappeared when the smallerincrement was used, showing that they were numerical in origin.Comparison of thepressures obtained for the two increments showed that atp = 0-6 the accuracy obtainedis approximately 1 %.Condition (ii) was significant only for the high density solutions in the super-position approximation. All solutions were obtained with R(max) = 4.0, and inaddition R(max) for p = 1.1 and 1.2 was extended to 5.0. The discrepancy inpressures between the two limits in the latter cases are 5 % for p = 1-2 and 2 % forp = 1.1. Condition (iii) was met by iterating until the desired convergence wasobtained, viz., a self-consistency of better than 0-5 % in the reduced pressure.of ‘C(x,r,d.TABLE 1.-YBG REDUCED PRESSURES PplpkT FOR A SERIES OF ‘C FUNCTIONS.PARENTHESES ARE NOT EXACT OWING TO SLOW CONVERGENCE.VALUES INTHE SUPERPOSITION PRESSUREAT p = 1-2 IS PROBABLY DIVERGENT.*1-4-6-7-8-91.01.11.2superposition 6 4-975 -9791.154 1-2151.537 2.013 - 2.9462.259 (5.06)27243.2293-6834,074obs.P64+ P265-979 -979 -979 - -9761.212 1.206 1 -226 1.30 1.211.973 2.079 1.912 2.03 1-892.803 - 2.55(4.44) 3.3604.3725.626Pad6 p&(HS) Kihara (Levelt)The accuracy of the z functions may be checked by evaluating virial coefficients.According to the theory, inclusion of terms in 1: up to S, should result in correctvirial coefficients up to Bm.For the superposition approximation B3 was calculatedto be 0.365 B2. The exact value is 0.365 2.32. B4 was found to be 0.104 B3. Thecorrect value is found, by interpolation of Barker’s results,s to be 0.116 B3.Similarlyfor BS we calculate 0.056 B4 and Barker’s value is 0.049 234. Here B = 2m3/3. Thediscrepancies cited reflect the limited accuracy of 84 and 8s.Table 1 contains values ofp/pkTfor 7 = 0 (superposition), ‘C 5 pS4, z = p&+p26~and z = PadC. AIso z = p&(HS) is included. These results are shown in fig. 2 anS . A . RICE AND D . A . YOUNG 213. The excess internal energy of the fluid may be calculated from the pair distributionfunction, i.e.,-!- NkT = LJrn 2kT u(r)g(r)4nr2dr. (10)The calculated excess internal energy U/NkT is given in table 2 for each of the casesIi' 8.06.0 '''1 DATAMONTE CARL0h45.0\4.03.02.01 .o -0.0 --uPO3.I .2 .3 .4 .5 .6 .7 .8 .9 1.0 1.1 1.2 1.3FIG.2.-Equation of state for the 6-12 potential fluid at kT/c = 274. The Monte Carlo data aretaken from ref. (3).4.0 -3.5 -3.0 -bl 2.5- %.52 . 0 -1.5 -1.0 -I I I I-I.4 .5 .6 .;I .8PO3FIG. 3.Detail of the fluid branch of the equation of state at kT/E = 2-7422 EQUATION OF STATE OF A MONATOMIC FLUIDTABLE 2.-YBG EXCESS INTERNAL ENERGIES U/NkT FOR A SERIES OF z FUNCTIONS. VALUESIN PARENTHESES ARE NOT EXACT OWING TO SLOW CONVERGENCE. THE S~PERPOS~~ION ENERGYATP = 1.2 IS PROBABLY DIVERGENT.obs.P superposition p84 p6,+,28, Yade p&dHS) Kihara (Levelt)--218 *1 -*223 -*222 --222 --222 --222 -*4 --877 -a859 --865 -*858 -a889 --744 --.829.6 -1.310 -1.267 -1.286 -1.269 -1.357 -1.118 -1.205-7'8 - 1.758 (- 1.625) -(1.630) - 1.863*9 -1.989 -2.1691.0 -2.234 -2.6161.1 -2.5171.2 (-2.84)- - 1 -455 - 1-477 -TABLE TH THE PAIR DISTRIBUTION FUNCTION Q(R) FOR Z =pdq+p2&.R0.10.20.30.40.50.60.70.80.91.01.11-21.31.41-51.61-71.81.92.02 12.22.32.42-52.62.72.82-93.03.13-23.33.43.53.63.73.83-94-0p = 0.40.00.00.00.00.00.00.00.00.1231.2361.6011 -4221.2191 -0780.99 10-9440.9260.9330.9600.9951 -0201 -0271.0231.0151 -0061 -0000.9970.9950.9960.9981 -0001 -0021 -0021 -0021 -0021.0011 -0001 *ooo1 -0001.000p = 0.60.00.00.00.00.00.00.00.00-1731 -5691.8541.5111 -2051,0080.8940.8440.8400.8800.9551 -0341 -0761 so731 -0481 ,0180.9940.98 10,9750.9790.9850.9991 a0071.0101 -0091 *OOG1 -0020.9980.9970.9970.9981 a000p = 0.70.00.00.00.00.00.00.00.00.2161.8352.0321 -5601.1830.9530.8270.7790-7880.8530.9641 -0721.1221.1071 -0621.0150.9800-9650.9590.9670 9841 -0021-01 51.0191.0161 -0091 -0020.9960-9940.9940.9961.000p = 0.80.00.00.00.00.00.00.00.00-2902.3372.3801.6401.1 120.8150-6680.6280.6660.79 11.0001 -2001 ~2721 -2091.0930-9890-9200.8940.8940-9260.9781 -0281 -0581 -0601 -0421.0150.99 10.9770.9730.98 10.9951.01S.A. RICE AND D. A. YOUNG 23cited, and is plotted in figs.4 and 5. For 0 <p G0.8, the superposition excess internalenergy is nearly linear in density, with an equation U,/NkT = -2.1975~03. Pairdistribution functions are listed in table 3.t20 MONTE CARLO DATA-. 2-.4-.69 8-1.0$ -1.4-1.6-1.8-2.0-2.2-2.4-2.6-2.8i-( -1.2.I .2 -3 .4 .5 .6 3 .8 .9 1.0 1.1 1.2 1.3p03FTG. excess internal energy of the 6-12 potential fluid at kT/E = 2-74. The Monte Carlo dataare taken from ref. (3).I 1 1 1 1 1 1 1oMONTE CARLO-1.1-1.2 ---1.6---1.7 SUPERPOSITION\-.8-\-.9-1.0-1.3- 2 -1.4-----1.5-1.7\-Is8 1.9 ir -*.+ , , , , , , , , , _I.4 .5 .6 .? .6p03FIG. 5 . D e t a i l of the excess internal energy of the fluid at /cT/E = 2-74.COMPARISON WITH EXPERIMENTThe 6-12 potential, although qualitatively sound, is inadequate to describequantitatively potential functions in real gases. One of the best pair potential func-tions is that due to Kihara.In particular, we have chosen the form used b24 EQUATION OF STATE OF A MONATOMIC FLUIDBarker :9U(R) = 4 e [ ( _ 9 > L 2 - ( z r ] R-‘1 R-‘1 .This potential function was inserted into the YBG equation with z = pS4(6-12)+p2a5(6-12) and the pressures and energies were calculated from the solutions of thisequation. We believe this to be a reasonable approximation since the main effects ofthe change in potential should be in the term zi(R), and not in the several 6,.Comparison with the argon data of Levelt 10 and the results for the pure 6-12potential show that the Kihara pressures are too large, and the energies too small (inmagnitude).In fact, the 6-12 results are better than the Kihara results. This is anillustration of the importance of the repulsive part of the potential in determining theequation of state of a dense fluid. The Kihara potential has a steeper repulsive partthan the 6-12 and this is reflected in the higher pressures. Evidently in a dense fluidthe Kihara potential is a less accurate description of the real effective potential thanthe 6-12.DISCUSSIONFrom the entries in the tables it is seen that the approximations z = 0, z = pS4,and z = pS4+p2& yield increasingly accurate pressures. For the internal energy,it appears that z = p64 is closer to the Monte Carlo curve than z = pS4+p265.Inboth the calculated pressure and energy, the Pad6 results lie nearer to the z = pS4curve than to the z = p64+p265 curve. Theoretically, the Pad6 approximant shouldlead to better pressures and energies than either of the other z functions. We believethat it is the inaccuracies in the 64 and S5 functions that are responsible for the failureof the Pad6 approximant to provide good agreement with experiment.The results with 64(HS) represent a “mixture” of the 6-12 and hard-spherepotentials. These calculations were motivated by the conjecture that the very steeprepulsive part of the 6-12 potential causes it to behave much like a hard core. Hencethe first-order correction using the hard sphere 6 4 should yield an improved equationof state.The results are indeed better than the superposition approximation, andsolutions may be obtained for the full range of densities. However, the excessinternal energies for 64(HS) are worse than for the superposition approximation. Thisemphasizes the qualitative difference between the 6- 12 and hard-sphere potentials.Also, the higher approximations z = p64(HS)+p265(HS) and z = Pad6 (HS) yieldpressures which are poorer than for p64(HS).From studies of the hard sphere and disc fluids,l, 11 the Kirkwood criterion for thestability of the fluid phase is not exact,l2 but coincides with the limit of stability forthe fluid in the superposition approximation (and for short-range modifications ofthis approximation). The limit of stability is given by the first real root of1 +?I: [kR cos kR - sin kR]u’(R)g(R,Q,p)dR = 0.By evaluating the left-hand side with a series of g(R), we obtain a function of kwith damped oscillations.This function resembles the Bessel functions obtainedfrom this equation for hard spheres and disks. At some density p,, the first minimumof the function will become zero, defining the upper limit of stability of the fluid.Since we do not have g(R,p) for a continuous range of densities, we must use inter-polation to locate the density cut-off. For the 6-12 potential, the superpositionapproximation leads to k , = 6.0 and pc = 1.12. For z = p64+p265, k, = 6.0 anS . A . RICE A N D D. A . YOUNG 25pc = 0433. Although these results are approximate, they are supported by the factthat we have found the solutions of the superposition approximation for p> 1.2 to bedefinitely divergent.At the highest fluid densities the YBG equation of state diverges markedly from theMonte Carlo equation of state.Although some of this discrepancy is probably dueto inaccuracies in the calculations, it seems likely that 8 4 and 85 are not adequate toaccount for all of the important long-range interactions in the 6-12 fluid. Clearly, adifferent approach is needed.We conclude that an adequate theory of the fluid must : (a) treat the short-rangerepulsions accurately, i.e., more accurately than by using the hard-sphere potential ;(b) take into account relatively long-range correlations, i.e., of the order of 3-5 atomicdiameters. A mean field treatment, superimposed on the 64, 8 5 fluctuation termsdescribed herein, is indicated.We thank the Directorate of Chemical Sciences, AFOSR for financial support.We have also benefited from the use of facilities provided by ARPA for materialsresearch at the University of Chicago. We are indebted to John Lekner for helpfuldiscussions. One of us (D. A. Y.) thanks the Hertz Foundation for generous financialassistance .1 S. A. Rice and J. Lekner, J . Chem. Physics, 1965, 42, 3559.2 B. J. Alder and T. E. Wainwright, J. Chem. Physics, 1960, 33, 1439.3 W. W. Wood and F. R. Parker, J. Chem. Physics, 1957, 27, 720.4A. A. Broyles, S. U. Chung, and H. L. Sahlin, J. Chem. Physics, 1962, 37, 2462.5 E. Meeron, J. Chem. Physics, 1957, 27, 1238.6 E. E. Salpeter, Ann. Physics, 1958, 5, 183.7 T. L. Hill, Statistical Mechanics (McGraw-Hill Book Co., Inc., New York, 1956). * J. A. Barker, P. J. Leonard, and A. Pomp, J. Chem. Physics, 1966,44,4206.9 J . A. Barker, W. Fock, and F. Smith, Physics Fluids, 1964,7, 897.10 J. M. H. Levelt, Physica, 1960, 26, 361.1 1 D. A. Young and S. A. Rice, to be published.12 J. G. Kirkwood, in Phase Transformations in SoZids, ed. R. Smoluchowski, J. E. Mayer, andW. A. Weyl (John Wiley and Sons, Inc., New York, 1951), chap. 3

ISSN:0366-9033    

DOI:10.1039/DF9674300016

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Exact and Approximate Values of the DistributionFunctions of a Simple Fluid*BY DOUGLAS HENDERSON tDivision of Applied Chemistry, C.S.I.R.O. Chemical Research Laboratories,Melbourne, Victoria, AustraliaSUNGWOON KIM AND LYNN ODENfDepartment of Physics, University of Waterloo, Waterloo, Ontario, CanadaReceiued 29th December, 1966The terms which are proportional to the cube of the density in the expansions of the pair distribu-tion function g(r), and the direct correlation fuiiction c(r), have been evaluated for the 6 : 12 potential.In addition, these terms have been evaluated using the Percus-Yevick (PY), hyper-netted chain (HNC),PY2, and HNC2 theories. Comparison of these approximate and exact values for g(r) and c(r) andof the resulting approximate and exact values of the fifth virial coefficient shows that, to this order inthe density, the PY2 and HNC2 theories are very reliable.Extensions of the PY and HNC theorieswhich are simpler thean the PY2 and HNC2 theories are also considered.If only pair interactions are considered,l the equation of state of a fluid of A7molecules in a volume V can be calculated from the pair distribution function g(r)by means of the pressure equation, or from the direct correlation function c(r) bymeans of the compressibility equation. These two functions are related by theequation :where c12 = c(r12), etc., p = N/V, and h(r) = g ( r ) - 1.The problem is to calculate g(r) and c(r). An approximate theory may be obtainedby postulating an equation, in addition to (l), relating g(r) and c(r).If such a theoryis tested by comparing the resulting approximate equation of state with experiment it isdifficult to make conclusions because of inexact knowledge of the intermolecularpotential u(r) and because of the perturbations resulting from the presence of three-body forces in real systems.An alternate procedure is to expand g ( r ) and c(r) in a power series inp and calculatethe coefficients both exactly and on the basis of the approximate theory. A conipari-son of these coefficients and the resulting virial coefficient Bh involves no assumptionsother than those inherent in the theory itself. The conclusions reached from such acomparison are only valid at low densities but they do not appear to be sensitive to thenumber of terms included in the expansion and are of interest in a discussion of theliquid state.hl2 = c12 +P.fh13C23dr3, (1)* This work has been supported by grants from the National Research Council of Canada and theOffice of Saline Water, United States Department of the Interior.t Alfred P. Sloan Foundation Fellow and Ian Potter Foundation Fellow ; permanent address :Department of Physics, University of Waterloo, Waterloo, Ontario, Canada.$ present address : Fresno State College, Fresno, California, U S A .2I>. HENDERSON, S . KIM A N D L . O D E N 27EXACT RESULTSWe denote the coefficients of pla in the expansions of y ( r ) = exp (pu(r))g(r) andc(r) by gn(r) and cn(r), respectively. The lowest-order coefficients are go = 1 andcg(r) = f(r) = exp ( -pu(r)) - 1.In previous publications 2, 3 we have presentedresults for gl, 9 2 , c1, and c2 using the 6 : 12 potential. In this paper we present our3 -- 1r*and 4 terms in the expansion for g(r).FIG. 1 .-g(r) for p* = 1 and T* = 2.5. The curves labelled 1, 2, 3, and 4 are obtained using 1 , 2, 3results for 93 and c3. Values for g(r), obtained using a four-term expansion, areplotted in fig. 1 for T* = kT/& = 2.5 and p* = pa3 = 1. This high density waschosen so as to magnify the contributions of the higher-order terms. In fig. 2 asimilar plot of c(r) is given. Unfortunately, our calculations for 93 and c3 are notcomplete so that it is not yet possible to show results for lower temperatures.PY A N D HNC THEORIESand 5 that c(r) can be written in the form,Eqn.(2) and (3) are little more than definitions of d(r) and E(r), and are useful only ifd(r) or E(r) can be approximated.The Percus-Yevick (PY) approximation is d(r) = 0 and the hyper-netted chain(HNC) approximation is E(r) = 0. Both approximations give correct results for goand g1 and thus give exact B2 and B3. However, they give approximate results forc(r) = f(r)y(r) +y(r) - 1 -In y(r) + E(r). (328 FUNCTIONS OF A SIMPLE FLUIDthe higher gn. In fig. 3 and 4 the errors in the PY and HNC values for g; = g2/a3 andg: = 93/06 at T* = 2.5 are shown. At this temperature the HNC theory gives theleast error for r* = r/a small but the PY theory gives the least error in the importantregion r * - 1.r*and 4 terms in the expansion for c(r).FIG. 2.-c(r) for p* = 1 and T* = 2.5.The curves labelled I, 2, 3, and 4 are obtained using 1,2, 3,I I 9 I O11 5 2.0 2.5 1.0r*FIG. 3.-A& = gz (approximate)-gg for T* = 2-5.In tables 1 and 2, B4 and B4, obtained from these theories, are compared with theexact values. The parameter b = 2na3/3. The necessary cluster integrals have beentabulated previously.39 6, 7 The compressibility equation yields the best B4 and B5 aD. HENDERSON, S . KIM AND L . ODEN 290.402*m 20- 0.2!!!!\.\ T * : 2.5\i;HNCZI /I1 I I1.0 15 2.0 2.5r*FIG. 4.-Ag4 = gf (approximate)-& for T* = 2-5.TABLE 1 .-FOURTH VIRIAL COEFFICIENTS FOR 6 : 12 POTENTIAL(in units of b3)0-6 -0.81.01.21.41.62.02-53.05.010.0182.17 -- 9.327- 0.2650,33930.27050.1 8980.12300.11310.1 1980.13420.1 156130-89- 7.179-0,31840.21700.20280-15910.12050.11600.12100.12990.1 120- 194.87- 11.699- 0.721 30.27890.30870.25290-18720.16580.16210.15340.1 177- 193.40- 12.597- 1.2473- 0.05020.07770.07650.06660.07750.09 120.1 1090.0937- 196.33- 10.768- 0.47760.33470.30430.23870- 18070-17390.18150- 19490.1681TABLE 2.-FIFTH COMPRESSIBILITY VIRIAL COEFFICIENTS FOR 6 12 POTENTIAL(in units of b4)T* exact PY PY2 HNC HNC2 v 0,0.60.81.01.21.41.62.02.53.05.010.0-- 76.6- 2.770-023- 0.0021- 0.041 5- 0.00990.03650.05790.06290.03901845.4 -- 44.8- 1.930.0020.03410.01 190-021 80.04410.05470.05500.03653818.1-81.6- 3.090.0090.0062- 0.0348- 0.00670,03870.05980.06420.0393-3689.8 -- 105.0- 6.61- 0.502- 0.01 300.02220.03700.05490-06090.04720.02 143836.8 -- 80.8- 3.020.00 1- 0*0040- 0.041 3- 0.00860.03870.06070.06650.04223614.0 -- 85.9-4.14-0.191- 0.0227- 0.02900.00330-04360-061 80.06450.041 8'3580.6- 99.2- 6.08-0.478- 0.0370-040030.02540.051 10-06140.055 10.03330 FUNCTIONS OF A SIMPLE FLUIDhigh temperatures while the pressure equation yields the best B4 and B5 at lowtemperatures. This is true of all theories considered in this paper.The PY theorygives ihe best overall results while the HNC theory gives the best results at lowtemperatures.PY2 AND HNC2 THEORIESVerlet S has proposed extensions of the PY and HNC theories.His expressionsfor d(r) and E(r) are too complicated to be given here. Oden el aZ.9 have given thefirst two terms in the density expansions of the PY2 ar,d HNC2 approximations tod(r) and E(r). Both theories give go, 91, and 92 exactly. In fig. 4 the errors in the PY2and HNC2 values for g; at T* = 2-5 are shown. The HNC2 theory gives the bestoverall values but the PY2 values are best in the most important region r* - 1. Valuesfor B5 which result from these theories are listed in table 2 ; they are very good.However, in the form given byVerlet, they involve an asymmetric approximation for the triplet distribution function,g(3)(rl,~,r3).This asymmetry does not affect 9 2 or 9 3 but will affect the higher-ordergn. A symmetric approximation which could be used is the superposition approxima-tion.10 However, if this approximation is used, the results for 93 are not as good asthose obtained from Verlet's asymmetric approximation. An alternate symmetricexpression isThe PY2 and HNC2 theories are promising.Eqn. (4) yields the same 9 3 as does Verlet's approximation to g(3) and should be moresatisfactory for the higher-order gn. We have not yet tested this approximation for the6 : 12 potential. However, our calculations 11 for Gaussian molecules and hardspheres indicate that eqn. (4) is a considerable improvement over the superpositionapproximation.SIMPLER EXTENSIONSThe PY2 and HNC2 theories are complicated. It would be useful to have somesimpler extensions of the PY and HNC theories.Each of these extensions yields anexact 9 2 and B4. Verlet 12 and Rowlinson 13 have proposed the approximation :We refer to this as the V approximation. In table 2 values for B5 which result fromthis approximation are listed. In general, the results are worse than PY values.However, at low temperatures the results are good and eqn. (5) may be of some valueat low temperatures. Green 14 has proposed an approximation which is similar to (5)but which yields poor values for Bg.An approximation would be to truncate the series expansion for E after one term.ThusEl2 = +P2J f13f14f23f24f34dr&4*This approximation, which we call the El approximation, corrects the HNC expres-sion for g2 but does not greatly affect the higher-order gn.This can be seen from table2. The El values for B5 are slightly improved over the HNC results. The advantageof this approximation is that it can be continued. When our calculations of 9 3 arecomplete it will be possible to add a second term to (6) and obtain the E2 approxirna-tion. The E2 approximation will yield g3 exactly and should give reasonably goodvalues for the higher-order gnD . HENDERSON, S. KIM AND L. ODEN 31A similar series of approximations could be based on the PY theory and could bereferred to as d l , d2, etc. approximations. However, in view of the superiority of theHNC results at low temperatures it would probably be better to use the HNC theoryas a basis.Van Leeuwen et aZ.5 have proposed an approximation where E(r) is givenby (6) with eachf(r) replaced by h(r). However, the 93 and Bs which result from thisapproximation are disappointing.The PY2 and HNC2 theories are promising theories of the liquid state. However,they are complicated and it is harder to obtain numerical values for g(r) from thesetheories than from the simpler PY and HNC theories. On the other hand, the E2and d2 theories, although probably less satisfactory than the PY2 and HNC2 theories,are no more difficult to use than are the PY and HNC theories.One of us (D. H.) is grateful to the staff of the C.S.I.R.O. Chemical ResearchLaboratories for their hospitality during his visit.1 J. S. Rowlinson, Reports Prog. Physics, 1965, 28, 169.2 D. Henderson, Mol. Physics, 1966, 10, 73.3 D. Henderson and L. Oden, Mol. Physics, 1966, 10,405.4 G. Stell, Physica, 1963, 29, 517.5 J. M. J. van Leeuwen, J. Groeneveld and J. de Boer, Physica, 1959, 25, 792.6 S. Kim, D. Henderson and L. Oden, J. Chem. Physics, 1966,45.7 J. A. Barker, P. J. Leonard and A. Pompe, J. Chem. Physics, 1966, 44, 4206.8 L. Verlet, Physica, 1964, 30, 95.9 L. Oden, D. Henderson and R. Chen, Physics Letters, 1966, 21, 420.10 J. G. Kirkwood, J. Chem. Physics, 1935, 3, 300.11 D. Henderson, J. Chem. Physics, 1967, (in press).12 L. Verlet, Physica, 1965, 31, 959.13 J. S. Rowlinson, Mol. Physics, 1966, 10, 533.14 H. S. Green, Physics Fluids, 1965, 8, 1

ISSN:0366-9033    

DOI:10.1039/DF9674300026

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Modified Cell Model for Liquids: Order Defectsand Intermolecular PotentialsBY FRIEDRICH KOHLBR AND J. FISCHERinstitute of Physical Chemistry, University of Vienna, AustriaReceived 16th January 1967On the basis of a cell model, an explicit introduction of coordination defects is possible. In binarymixtures, consideration of coordination defects gives only a small correction to the thermodynamicexcess properties of mixing. The excess free energy is lowered, and its concentration dependenceis influenced in about the same way as by using the quasichemical approximation. Furthermore,the modified cell model has been used to calculate the thermodynamic properties of solid and liquidargon, employing a Kihara pair potential. Comparison is made with previous calculations using aLennard-Jones pair potential.The parameters of the Kihara potential giving the best fit at highdensities differ markedly from those giving the best fit for transport properties and second virialcoefficients. The influence of the exact form of the pair potential on the thermodynamic properties isnot very important.The modified cell model 1 9 2 combines some development of Kirkwood’s varia-tional method 3 with an explicit consideration of order defects. The latter areresponsible for density fluctuations, by analogy with Kirkwood’s 3 undeterminedparameter for multiple occupancy of the cells. Therefore, we have carried out thedevelopment of Kirkwood’s variational method for single occupancy of the cells,and then consider order defects.In developing Kirkwood’s variational method, we have had two aims: (i) toaccount for the distribution of the locations of the neighbouring molecules withintheir cells, (ii) to allow a certain correlation of the motion of the central moleculewith the motion of the neighbouring molecules, at least near the edges of the cells.The first point has been considered by Mayer and Careri.4 As we want to employthe Lennard-Jones 6 : 12 pair potential or the Kihara potential, we have used asquare-well approximation for the probability density in the neighbouring cells ratherthan the Gaussian approximation used by Mayer and Careri.Like them, we haveadjusted the size of the square well consistently to the free volume available in thecentral cell. The second point has been raised by Hirschfelder, Dahler, and Thacher,swho have pointed out that without such a correlation the cell model will always givetoo positive an energy and too small a free volume, especially at low densities.Where-as Kirkwood 3 has based his variational method on the assumption of independentmotions of molecules in their cells, our method shows a certain degree of dependenceupon these motions. Therefore, it is necessary to introduce the correlation functionbefore minimizing the free energy. Our starting equation isN k=N- 1Here P(r1,rz. . . rN)drldr2. . . drN denotes the probability of finding the first particlea the dislocation rl of the centre of the first cell, the second particle at thedislocation r2 of the centre of its cell, etc.+(rl)drl is the probability that a particle is3F . KOHLER AND J . FISCHER 33at a distance r from the centre of its cell, provided that no correlation with the motionof neighbouring molecules is taking place. The effect of the correlation functionx(rik), which is supposed to be a function only of the mutual distance of pairs ofparticles, is to overweigh pairwise configurations which are strongly attractive and toexclude pairwise configurations which are strongly repulsive. The correlation islimited to small mutual distances. For higher distances, the x are set arbitrarilyequal to zero. The limitation is made in such a way that the product n(l + x i k ) canbe developed into various sums, each containing factors with different subscriptsonly, analogous to the treatment of real gases at low densities :n(l + X i k ) = +&ik + C & i k h n + * * (2)i > kInserting eqn.(1) and (2) into the expression for the free energy, and minimizing itwith respect to each 4 and x, these functions can be calculated. We have approxi-mated the probability density 4 by a square well, and the correlation function x by (x = -1) for small pairwise distances, and by (x = 0) for larger distances. Byusing the minimizing procedure a consistent size for the square well is establishedand a value of the pairwise distance where the step should occur is obtained.An order defect is thought of as the removal of a nearest neighbour into the shellof next-nearest neighbours. In this process, the distance between the centres ofneighbouring cells has been stretched from the nearest neighbour distance a to thedistance A = 23a.Let Fa be the free energy for a nearest neighbour distance aaccording to the model above, without order defects. Let FA be the same for anearest neighbour distance A . Then the free energy for an assembly with nearestneighbour distance a, but with a fraction x of stretched contacts, is assumed to begiven byHere [&)IN denotes the number of possible arrangements due to order defects.The fraction x is again determined by a minimization process. For densities lowerthan the critical density, x becomes equal to about 0-5 ; R In g(x) is then the communalentropy per mole of gas. Including a correction for the void of a face-centeredcubic packing, In g(x) has been written1 :In eqn.(4), z denotes the coordination number.F = (1 -x)F,+xF,-RT In g ( x ) (3)In&) = -(2/2)(xln x+(l -x) In (1 -x)+1,88x(1 -x)} (4)BINARY MIXTURESThe model for the order defects may be combined with simpler (and poorer)approximations of the cell model than the one described in the introduction. We nowinvestigate the effect of our coordination defects on the thermodynamic excess pro-perties of mixing of simple liquids. For this, we use the cell model for single occupancyin the approximation given by Prigogine and coworkers6,7. Again, a square well isassumed for the probability density, but the square well is arbitrarily extended tothe distance (a- 0 ) / 2 from the centre of the cell (a being the diameter of the cell, andD being the mutual distance between two molecules where the pair potential passesthrough zero).Furthermore, the energy is calculated as if all molecules would belocated at the centres of their cells.In mixtures, a given molecule A is surrounded by both A and B molecules.Therefore, the energy parameter A,, = ZE& (E& . . . minimum of pair potential betweenA molecules) applicable in the pwe A liquid has to be changed to a parameterAa = EAaa + (1 - X)Aab in the mixture. In the same way, we have for cells occupied by34 MODIFIED CELL MODEL FOR LIQUIDSa B molecule an energy parameter A b = EAab + (1 - x)l\bb. Here E denotes the overallmole fraction of A, and random mixing is implied. For non-random mixing, themole fraction applicable for the construction of A, would differ somewhat from themole fraction to be used in the expression for &.Introducing the concept of stretchedcontacts, we define the following quantities : x, and x b denote, for A-cells and B-cellsrespectively, the fraction of neighbouring molecules which have been removed intothe shell of next-nearest neighbours (i.e., the fractions of stretched contacts).Starting from an A-cell, a fraction y of the normal contracts, and a fraction 7 of thestretched contacts will lead to other A-cells. Starting from a B-cell, the correspondingfractions are called 5 and F. The various quantities are interrelated by the followingrelationships : (a) mass balance equations of the form,y( 1 - xu) + yx, = x, ( 5 )z ( l - x a ) ( l - y ) = ( l - z ) ( l - x b ) t , (6)(b) identity between AB and BA contacts like(c) an assumption about the molecular distribution.We have tried as alternativesthe statistical distribution,and the quasichemical approximation 8m$/(maambb) = 4r7*Here N,b,nm, and N b b denote the numbers of AB, AA, and BB normal contacts,respectively, and q stands for exp(&b - Am- I\bb)/(ZkT).The starting point is the equation for the free energy of the mixture :F = - kT In [Ul,".C' -xa)QNaxaYNb(l -Xb)aNbxba b b 1 + ( N a / 2 ) [ ( 1 -xa)Ea+xaE,l fNa and Nb are the numbers of A and B molecules. Y, and Ea denotes the translationalpartition function within the free volume and the energy of an A molecule in a cell ofdiameter a,, a, and E, stands for the same quantities in a cell of diameter aa2*The double bar means that the number of contacts corresponds to the statisticaldistribution.The free energy is essentially a function of the size of A- and B-cells, and of thefraction of stretched contacts F(a,,ab,xa,xb). In order to determine the parameters,F has to be minimized subject to the condition that the volume, which is a functionof the same parameters, remains constant : V(aa,ab,xa,xb) = const.This leads toequations of the typeaF av -+I-- = 0,aab dab (9)where the Lagrange multiplier A can be evaluated from one of these expressions, e.g.F . KOHLER AND J . FISCHER 35Then each of the partial differentials can be expressed by one of them, say, by a,F/da,,and the equation of state can be written asAs long as P can be assumed to be practically zero, each partial derivative of the freeenergy with respect to the parameters vanishes, as the derivatives of Y with respectto the parameters are all different from zero.Thus we get-by analogy with thetreatment of Prigogine and Bellemans6 -the conditionsd F / W = -Y = 4 = (i3Fjaaa)(t3aa/aV). (11)i?F/i?a, = aF/aa, = dF/dxa = dF/dx, = 0. (12)The computation proceeds as follows. First, estimated values of xu and xb areused to calculate the values of y,y,c and f . Then, with a first guess of ab, values forQ, and xu are determined in small iteration cycles. On the basis of these values, thesame iteration cycles are used for the determination of a b and xb. Now a new set ofthe y and 5 is calculated, and the whole procedure is repeated. The most sensitiveparameters are xu and x b .In our example, about five cycles have been necessaryto arrive at values constant to 0-001 %.n 3 2500E5I1d0lg 2450.4 W'5. wc52 4 0 00.0 0.1 0.3 0 . 5 0.7 0.9 -XFIG. l.-GE/[X(l-x)] as calculated for 293-15°K with the parameters given in the text. Upperpair of curves : no coordination defects considered, statistical and strictly regular distribution, respec-tively. Lower pair of curves : the same with coordination defects.The example presented here has been calculated with constants appropriate tothe mixture 1,2,4-trichlorobenzene + n-hexane 99 10 (A, = 6623 1.5 ; Aab = 53 124.7 ;-Abb = 43425.7 J/mole; RZ, (the distance at the minimum of the pair potential of Amolecules) = 6.6466 x 10-8 ; R& = (R&+ R&)/2 ; Rib = 6-6782 x 10-8 cm).Thetemperature variation of Am, caused by the dipole-dipole contribution, has beenneglected. Fig. 1 shows the excess free energy of mixing GE, divided by the productof the mole fractions. The upper pair of curves is calculated without coordinationdefects (Le., xa = x b = O), for statistical and strictly-regular distribution, respectively.The lower pair of curves show the results of the consideration of coordination defects.Their effect on the absolute value of GE is not very marked, which is satisfying fromthe point of view of current theories of mixtures. The main effect of the coordinationdefect is the enhanced cilrvature of the GE/[-F(I -X)] function.For this detail, theconsideration of coordination defects is at least as important as the effect of thequasichemical approximation. As the curvature of the GE/[X( 1 - X2] function isrelzted to the shape of consolute curves,ll one might expect that the consideration o36 MODIFIED CELL MODEL FOR LIQUIDScoordination defects will facilitate the explanation of the consolute curves foundexperimentally. On the other hand, our calculations-which are based on theassumption of central force fields of the molecules-have given no indication of aninversion of the curvature of the GE/[X( 1 - E)] function at high Z values, as has beenfound e~perimentally.~ Fig. 2 gives the concentration dependence of the fractions ofI I0.1 0.3 0.5 0.7 0.9 - XFIG.2 . T h e fraction of stretched contacts for A- and B-cells, respectively, corresponding to thelower pair of curves in fig. 1.stretched contacts, starting from A- and B-cells. Finally, the order of magnitudeof the TSE/[X(l -X)] function for statistical distribution is - 150 J/mole deg without,and + 70 J/mole deg. with, consideration' of coordination defects. Introduction ofthe quasichemical approximation makes this function more negative, the value forX = 0.5 being -20 J/mole deg.COORDINATION I N DENSE SUPERCRITICAL VAPOURSAnother possible application of the model of coordination defects is based on thefact that volume expansion may be caused primarily either by diminishing thecoordination or by extending the average distance to the nearest neighbours.Fig. 3shows the result of calculations for argon at O"C, which demonstrates that up to avolume, about twice that of the liquid at the boiling point, expansion is causedprimarily by diminution of the coordination. The nearest neighbour distanceincreases only about 10 %. One might speculate and correlate this with the findingsof Hensel and Franck l2 on the strong increase in conductance with pressure foundin the dense, supercritical vapour of mercury. The assumption seems to be reasonablethat electron transfer occurs only between atoms which are below a certain minimumdistance. Fig. 3 shows that there might be a large range of volumes, where normalcontacts are less than this minimum distance and stretched contacts are greater.Then the conductance would be determined by the number of interruptions whichoccur in a row of normal contacts.By interruption we mean that the row of normalcontacts has either no continuation or the continuation would result in an opposingpotential. Therefore, the last atom must be the starting point for nine stretchedcontacts at least. The specific resistance would be essentially proportional to thP. KOHLER AND J . FISCHERninth power of the fraction of stretched contacts. In detail,13x9( 1 - x)( 1 + 2x2)'37Here IC is the specific conductivity, which is referred to the nearest neighbour distance13 0 4 0 5 0 6 0 7 0 8 0V(cm3)FIG. 3.-The nearest neighbour distance a and the mean coordination number 2, as calculated forargon at 273.15"K.T4 12 2 0 28PIP,FIG.4 . T h e r.h.s. of eqn. (13) for argon at 273.15"K.a, and the subscript 1.d. stands for "low density", where transferable electrons areabsent. The quantity on the r.h.s. of eqn. (13) depends on the pressure accordingto fig. 4, where the calculated results of the 0°C-isotherm of argon are used. It isquestionable to transfer the results from argon to mercury, although it demonstrateshow important the distinction can be between volume expansion by coordinationdefects and volume expansion by increasing the nearest neighbour distance38 MODIFIED CELL MODEL FOR LIQUIDSPAIR POTENTIAL AND THERMODYNAMIC PROPERTIESIn the previous paper,2 the modified cell model has been applied to argon at higherdensities, using a Lennard-Jones 6 : 12 pair potcntial with the minimum coordinatesEO = 128 k, R" = 3-83 A.In the present paper, a Kihara potential&(r) = &0[(-)"-2(->3 R"-s R"-sr - s r - shas been employed with the parameters of Barker, Fock, and Smith 13 (i.e., E' =The free energy of the solid is calculated without order defects, but by consider-ing coupling of the vibrations (e.g., by employing a Debye approximation).Previously,2 the calculated melting point was too low (55°K in comparison with theexperimental value of 84°K). We thought that possibly the Lennard-Jones potentialwas responsible for the free energy of the liquid being too negative, as the Lennard-Jones potential overemphasizes the attraction at distances larger than R".142.9 k, R" = 3.7338 A, s = 0-90069 R").1.0 1.2 1.4 1-6PFIG.5.-The difference of the Gibbs free energy between solid and liquid, as calculated with aLennard-Jones potential (L.J.) and with a Kihara potential (K.).FIG. 6.-The gas-liquid coexistence curve of argon at high densities.The present calculation shows that the effect of the form of the pair potential onthe thermodynamic properties is small in comparison to the effect of the choice of theminimum coordinates, and in comparison to statistical approximations. Thedifference between the free energy curves of liquid and solid is not greatly changed byusing the Kihara potential, and the small change is in the wrong direction. Fig. 5gives the difference of the Gibbs' free energy between solid and liquid according to thepotential used.Entropy and volume difference between solid and liquid phase areabout correct at the temperature of the experimental melting point. Furthermore,fig. 6 shows part of the coexistence curve between liquid and gas at high densities.The parameters of the Lennard-Jones potential used correspond to a ratio E"/R",which is rather high compared to the values most frequently given. The parametersof the Kihara potential are supposed to fit the measurements of gas viscosities andsecond virial coefficients. Nevertheless the use of the Kihara potential gives a smallervolume and a larger energy of solid and liquid than the use of the Lennard-JonespotentialF . KOHLER AND J . FISCHER 39Fig. 7 shows the entropy of the solid as function of the temperature.This functionseems to be most seriously influenced by the square-well approximation for the proba-bility density, which leads to characteristic temperatures rising with temperature.If a constant characteristic temperature 0" is used (taken from the curvature of25h8 -23\ 3 eFIG. 7.-The entropy of solid argon; (L.J.) and (K.) calculated with a square-well approximationfor the probability density in the cells, using a Lennard-Jones and Kihara potential, respectively;0" calculated with a characteristic temperature evaluated directly from the hnnard-Jones potential.FIG. 8.--The volume of solid argon ; descriptions of curves as in fig. 7.the Lennard-Jones potential around the potential minimum), then too high an entropyis obtained. In the square-well approximation, the Kihara curve should differ fromthe Lennard-Jones potential around the potential minimum), then too great an entropyand this is the case. Finally fig. 8 shows the volume of the solid as function of thetemperature. Qualitatively, the situation is the same as for the entropy, but thevolume according to the Kihara potential is much too small.1 F. Kohler, Ber. Bunsenges., 1966, 70, 1068.2 F. Kohler and F. Weissenback, submitted to J. Chem. Physics.3 J. G. Kirkwood, J. Chem. Physics, 1950,18, 380.4 J. E. Mayer and G. Careri, J. Chem. Physics, 1952, 20, 1001.5 J. S. Dahler, J. 0. Hirschfelder and H. C. Thacher, J. Chem. Physics, 1956, 25, 249.6 e.g., I. Prigogine and A. Bellemans, Disc. Fcruday Soc., 1953, 15, 80 ; I. Prigogine, MolecularTheory of Solutions (North Holland 1957).7 see also, F. Kohler, Monatsh, 1957 88, 857.8 E. A. Guggenheim, Mixtures (Oxford, 1952).9 R. J. Munn and F. Kohler, Monarsh, 1960,91, 381 ; A. Neckel and F. Kohler, Monatsh, 1956,10 F. Kohler, Chem. Technik, 1966,18,272.11 F. Kohler, J. Chem. Physics, 1955, 23, 1398.12 F. Hensel and E. U. Franck, Ber. Bunsenges, 1966, 70, 11 54 ; we are indebted to the authors13 cited from J. S. Rowlinson, Disc. Furaduy Soc., 1965,40, 19.87, 176.for communicating their results prior to publication

ISSN:0366-9033    

DOI:10.1039/DF9674300032

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Calculation of Thermodynamic Properties of Liquid Argon fromLennard-Jones Parameters by a Monte Carlo MethodBY I. R. MCDONALD AND K. SINGERDept. of Chemistry, Royal Holloway College (University of London),Englefield Green, Surrey.Received 16th January 1967Extensions of the Monte Carlo method of Metropolis et al. are described which permit bothisochoric and isothermal extrapolations of Monte Carlo data. Thermodynamic properties arecalculated for liquid argon from Lennard-Jones parameters for a wide V, T-range. The agreementbetween calculated and experimental values is, on the whole, satisfactory.The Monte Carlo (MC) method devised by Metropolis et aZ.1 has been applied tothe calculation of equilibrium properties of a variety of systems.1-12 Of particularinterest in the present context is the work of Wood et al.who calculated the radialdistribution function, internal energy, pressure and specific heat of argon over a largerange of densities 4 at 55”C, and for one density 5 of the liquid at - 147°C fromLennard-Jones (LJ)-parameters.This paper explores the possibility of calculating the principal thermodynamicfunctions of a simple classical liquid from a LJ-potential over an appreciable range ofvolume and temperature. The project has become practicable (a) because the greaterpower of computers now available permits the execution of computations of the typereported in ref. (4) without excessive cost in computer-time ; and (b) because the MCmethod has been modified so that properties for a certain range of V, T-values can beobtained from MC data for one V, T-point.Argon is chosen for this work becausereliable LJ-parameters 13 and excellent P- V-T data 14-16 covering the entire liquidrange are available.MONTE CARL0 METHODThe essential features of the method of ref. (l), which has been discussed,lp 49 17are these. A fluid “sample” is represented by N coordinate triplets confined withina cube of volume v. The potential energy @, assumed to be a sum of pair energies,is calculated for the initial configuration. One of the N “particles” receives a randomdisplacement and the resultant change A@ of the potential energy is determined. IfA@ is negative, the new configuration is accepted ; if it is positive, the move is acceptedonly with the probability p = exp(-PA@), l/kT.If the move is accepted, thenew configuration and the new potential energy become the “current” properties ;if it is rejected the system returns to the previous state. The sequence of configurationsgenerated by repetition of this procedure has been shown to tend towards a Boltzmann-weighted ensemble. The mean value of any function of the 3N coordinates over asufficiently long MC sequence is therefore an approximation to the ensemble averageof this quantity. Periodic boundary conditions are imposed to minimize the error4I . R . MCDONALD AND K . SINGER 41resulting from the small size of the sample, which, for reasons of economy, usuallycannot exceed N - 102.The mean potential energy (@> so obtained determines the molar configurationalinternal energy Ut ; the mean virial,( Y ) = <xii.via>idetermines the pressure according toIsothermal changes of the Helmholtz free energy can be calculated from (">-valuesand isochoric changes of the configurational entropy from (@)-values byPV = R T -(No/3N)(Y>AA = -(NOIN) Pdv (1) and ASt = CLdT/T, CJ = (No/N)(a(@}/aT),. (2)In these and subsequent formulae, N is the number and v the volume of the MCsample ; molar quantities are denoted by NO and V ; the dagger 7 means "configura-tional" and the angular brackets (> indicate ensemble averages.The present extension of this method is designed to evaluate the classical con-figurational partition function in the formQt(N,v,T) = (VN/N. l)rmaxg(0) exp ( - pO)d@,f:: 1::Q, minwhere g(CD)dCD is the fraction of all accessible configurations for which the potentialenergy has a value between @ and @+do.Boltzmann-weighted sequences ofconfigurations and averages of @ and Y are calculated as in ref. (4) ; but in addition,a histogram is accumulated of the numbers of configurations with potential energieswithin small specified ranges of width 6. The entry for the n-th interval is the numberof configurations with potential energies between CD, - 6/2 and Qn + 612, denoted byFor sufficiently small 6, Fn is assumed to converge statistically toFn F(Qn - 6/2,mn + 612).where K depends on v but not on T. Sincethe histogram determines the partition function apart from the normalizing constantK .Qt at another temperature, T' can be obtained by reweighting the histogram:where p' = l/kT'.In practice, the range of temperatures which can be reached by thisprocedure is limited by the width of the parent histogram.Since K cannot be determined by the present method, isothermal changes of Qfare calculated indirectly, by means of (1). To compute changes along a neighbouringisotherm T', a reweighting procedure is applied. Let {a), be the set of G, con-figurations for which the potential energy lies between @,-a12 and @,,+6/2, Yl.n,orthe virial of one of these configurations, and Tn the average for the set (a)n, i.e.42 THERMODYNAMIC PROPERTIES OF LIQUID ARGONThe ensemble average of the virial at Tis(yy>T = CF,(T)qn/(xFn(T))n nand the value (Y)Tt is obtained by substituting Fn(T) = Fn(T)exp[-(D’-Q)@,J forFn(T) in (6).Tables of qyr,-values are calculated during the MC runs.Another extension of the range of application of the MC method arises from thepossibility of calculating from MC data obtained for the LJ-parameters ( E , a) propertiescorresponding to different parameters (d, a’). This can be done if the differencesbetween the “new” and the “old” parameters are small.The LJ-potential is defined by$(r) = 4~[(a/r)’~ - ( ~ / r ) ~ ] (7)modified by a hard core, i.e., $(r) = 00 for r < rc, introduced to increase the speed ofcomputation. Since the potential function and the corresponding virial function,$(r) = rd4/dr, are linear combinations of r-12 and r-6, the potential function 4’(r)and the virial function $’(r) for the parameters (E’, a’) can be written as linear com-binations of $ and $.The same relationships hold for the total potential energy@n,a and the total virial Yn,a of any configuration a belonging to the set ( c l I ncalculated with different LJ-functions, i.e.,= a@,,,+ b v n , a , (8)(9) y ; , a = c@n,a + d y n , a ,wherewith z = E’/E and A = d/a.The ensemble average for the “new” potential energy is<af> = CC x (a@n+byn,J ~ X P [-D(aQn+byn,a)I)/(C C ~ X P [-D(aan+byn,a)J>,n taln n {ahwhere @n,= is replaced by an since @n,ag@n. The further substitutionsN -Yn,a f Yn+&, a@,+ b!Pn f @n,lead to<Q‘) = {C c ( G n + bAn,a) ~ X P [ - B ( G n + b ~ n , a ) l ~ / ( ~ c ~ X P C - B ( G n + b ~ n , a > l ) * (12)n Ialn n {ahEqn.(lOa) shows that b -g 1 if 10-0‘1 $a ; in this case exp ( -#lbAn,a) may be expandedand terms of higher than the second order in pbAn,a neglected. If the substitutionsare made in the resulting expressions, one obtains(1 3) (a’) = { E F ‘ n [ 6 n ( l + p 2 b 2 ~ / 2 ) - B b 2 ~ l ) l ( C ~ n ( 1 +B 2 b 2 2 An/2))n nThe analogous formula for the virial iswhere yn=c@n+dYn. Eqn. (13) and (14) are exact if the distributions of A%,a arI . R. MCDONALD AND K. SINGER 43Gaussian. If j P b 2 z < 1 for all n, one may use the simpler expressions,(Y‘) E (CFn[‘u,-pbdii3)/CFn. n (16)n -The values of YYn may be calculated during MC runswithout loss of computing speed.These formulae can be used (a) to investigate the effect of changes of the LJ-parameters on the equilibrium properties, and (b) to extrapolate data calculated foru,T to u‘,T’. Eqn.(1 3)-( 16) are applicable to isothermal changes because substitutionof the LJ-parameter 0 by 0’ = a/A is equivalent to the scaling of all linear dimensionsof a sample configuration by A, i.e., r i p d r i j for all i, j and u-u’ = A321. Changes ofE to E’ are equivalent to changes of Tand can be dealt with by the reweighting procedureaccording to (4) and (6). The isothermal extrapolation is limited to small changes ofvolume because, apart from the approximation made in the derivation of eqn. (13)and (14), the reweighting of the histogram Fn to becomes inaccurate if the twohistograms do not overlap sufficiently.MC data for the calculation of thermodynamic properties may also be obtainedby a method not based on ref.(1). In this, a sequence is generated by a random walkin 3N-dimensional configuration space without Boltzmann-weighting. A move isrejected only if it raises the potential energy above a preset upper bound. Thehistogram G(Qn - 6/2, Qn + 6/2) = Gn so obtained is assumed to tend to g(Q) in (3), andthe partition function is calculated fromi.e., the variances ofi.e., by Boltzmann-weighting of the histogram according to the temperature required.The limitation of this method of “unweighted configurations” arises from the veryrapid increase of g(@) with a. A considerable number of MC runs with differentupper energy bounds is therefore required to cover a @-range sufficient for the calcula-tion of properties over a wide T-range.Isothermal changes and extrapolations canbe dealt with as described above. The disadvantages of the “unweighted” methodincrease with density and with N ; but they are not too serious in the gaseous range,even at high pressures. Results obtained by this method for gaseous argon over awide P-V-T range will be reported elsewhere.18 The two methods have not beencompared systematically .RESULTS AND DISCUSSIONCalculated values of selected thermodynamic properties of argon at five tempera-tures between the triple point (8343°K) and the critical temperature (150-7°K) and arange of densities are shown in table 1. The experimental results for the same V, T-points are given where these are available ; properties of the orthobaric liquidobtained by graphical interpolation are shown in table 2.Experimental P-V-T dataat 86.3 and 97.0 OK are derived from the empirical equation of state given in ref.(15) ; those for 108.2, 136.0 and 14802°K by interpolation of the data of ref. (16).All calculations are based on samples with N = 108 and on the LJ-parametersElk = 119*76”K, 0 = 3-405&13 with rc = 2-7A (cf. eqn. (7)). Not less than 168,000configurations were used to obtain the data for a given u,T-point. The mutualconsistency, the broad agreement with experiment and the magnitude of the statisticalfluctuations of the calculated mean values indicate that MC sequences of this lengt44 THERMODYNAMIC PROPERTIES OF LIQUID ARGONusually yield reliable statistical data.The standard deviations (calculated fromsub-averages over 4000 configurations) are m0-3 x 10-15 ergs for (<DIN) and =2 to4 x lO-.15 ergs for ('FIN). The relative error in the virial is therefore largest whenT"K86.397.0a97*Ob108.2c136*Od148-2av <;>cm3 ergx 101424-78 11-9025.60 11.512644 11.1526.90 10.9327.50 10.5627.50 10.3028.03 10.1228.48 9.9528.95 9.7629.68 9.5030.92 9-1428.48 9.7529.68 9.3627-50 10.1428-03 9.9728.48 9.8128-95 9.6229.68 9.3828.48 9.5929.68 9-2433.51 7.9736.58 7.3340.00 6.7233.51 7.8735-00 7.5836.58 7.2538.19 6.9540.00 6.64TABLE 1 .-THERMODYNAMIC PROPERTIES OF ARGON-1 avcm2 dyne-1 dyne cni-2 deg-1F(B)T (.E)V=(E)TY - c;x 10-7-ui 12 R $P('> dyne cm-2 x 10-8 xefop14 MC a p t .MC expt. MC expt. ~ 1 0 1 0MC expt. MC expt.1.64 1.57 9.98 0.91 1.60 1 *744.12 -0.43 9.66 0.96 1.66 1.82668 -2.36 9.36 1.27 1.71 2.447.32 -2.79 9.17 1 *236.36 -2.03 8.85 1 *280.66 2-13 2.10 864* 1*35* 1.87 1.73 2.30 1.942.51 0.76 0.92 8.49 1.03 2.17 1.91 1.73 2.073-34 0.17 0.09 8.35 1.09 2.34 2-06 1-77 2.154.18 -0.42 -0.67 8.19 1.02 2.50 2.19 1.72 2.225.55 - 1.34 7-97 0-85 2-76 1 *407.58 -2.60 0.68 3.10 1.05-0.14 2.93 2.16 7-28* 1*00*2.91 0.81 0.35 6.99 0.9 1-2.21 4.55 4.24 7.57* 0*88* 1.88 1.94 1.52 1.210.01 2.87 3.01 7.44 0.93 2.18 2.14 1.47 1.390.85 2.23 2.16 7.32 0-83 2.31 2.29 1.42 1-501.60 1.68 1-39 7.18 0.91 2.45 2.43 1.50 1.593.43 0.40 0.35 7.01 0.71 2-71 2.64 1-18 1.69-2.87 5.61 6-41* 1*13*0.78 2.54 229 6.18 0.781.42 2.52 2.37 4.24 4.15 0.56 0.58 4.73 4.53 0.83 1.214.20 0.79 0.95 3.90 3.82 0.47 0.57 6.70 8.30 0.60 0.975.88 -0.12 0.28 3.58 0.46 0.58 17-31 23.50 0.56 0.61-0.35 3.89 3.85 3.75 0-61 0.61 2.89 0.881.68 2.49 2.79 3.69 3.60 4.27 4.722.80 1.93 2.01 3.54 3.46 0.50 0.57 6-50 6.58 0.64 0.893.38 1.45 1.52 3.39 3.32 9.15 9.664.54 0.80 1.13 3.25 3.18 0.44 0.57 15-42 16.73 0.54 0.65MC = present work ; expt = experimental results ; * see text.deg-1 x 103MC expt.2.783.024.17430 3.362.75 3.954-14 4.434.30 4863.863.322.86 2.353-20 2.973.28 3.443-68 3.863.20 4.463.91 5.264.03 8.089.76 14.292.554-16 5.868.25 10.87direct MC calculation ; reweighted from 86.3"K data ; reweighted from 97.O"K data ; reweighted from 148~2°K data.( Y / N ) is small, i.e., when PVMRT and the relative error in P is largest when P%O( ( Y / N ) M 3kT).Estimated probable statistical errors are : 6U+* 5 cal/mole and6P M 2 x 107 dyne cm-2. The isochoric reweighting gives reliable results for variationsof up to 15 %in T.PRESSUREThe variation of PV/RT is shown in the figure ; with the exception of one point(28.48 cm3, 97-0"K) the agreement with experiment may be regarded as satisfactorythroughout the range of stability of the liquid. At 86-3"K the calculated pressuresform distinct "solid" and "liquid" branches. Although the pressure calculated forV = 27.50 cm3 from an initial face-centred cubic arrangement remained much belowthat calculated from an initial disordered liquid-type configuration, the pressure ofthe "solid" is anomalous, being higher than that calculated for Y = 26.90 cm3.This arises almost certainly from the presence of liquid-type configurations or regionsin the MC sample, indicating a fluctuation between two pha~es.4~10 At 148.2"Kand Y = 40.00 cm3 the calculated pressure falls well below the experimental valueI .R . MCDONALD AND K . SINGER 45This is interpreted as the onset of a van der Waals’ loop, a view supported by moreextensive computations 18 in this region with N = 32. The values of P at 9700°Kobtained by isochoric reweighting of the 86.3”K-data are in particularly good agree-ment with experiment.The reason for the discrepancy between these results and1.51.00 5FY w s4 0-- 0 5- 1.0-----35I 1 1 1 1 1 1 I I I24 26 28 30 32 34 36 38 4 0molar volume cm3FIG. 1 .-D- V isotherms for argon. 0, data from direct MC calculation ; 0 , data from isochoricreweighting. Full curves : expt. results ; dashed curve : sketched through MC points. Curve 1 :86~3°K (“solid” branch) ; curve 2 : 86.3”K (“liquid” branch, expt. equation of state extrapolated tonegative pressures) ; curve 3 : 97.O”K ; curve 4 : 136-0°K ; curve 5 : 148.2”K.the value obtained directly from MC runs for 28.48 cm3 is not known. Inspectionof the statistical data reveals signs of long-term fluctuations which suggests thatsubstantially longer sequences than 170,000 configurations are occasionally required.INTERNAL ENERGY AND LATENT HEATSThe agreement between observed 14 and calculated values of Ut at 148.2 andComparison with experimental data for the lower temperaturesThe latent heat AH,,, of vaporization at 86~3°K can be136~0°K is good.is made in table 2.TABLE 2.-pROPERTIES OF ORTHOBARIC LIQUID ARGONT P V - Uf/RT C p OK dyne cm-2x 10-6 cm3calc.expt. calc. expt. calc. expt.86.3 0.92 28.62 28.52 8.34 8-19 1 a07 0.8597.0 2-59 29.91 29-95 6.97 6.86 0.65 0.78108.2 5.97 31.13 31-71 5-92 5.69 0.35(?) 0.65calc = graphical interpolation of MC data of table 1 ; expt = linear interpolation of experimentalresults.21calculated from the experimental vapour pressure 20 (0.91 x 106 dyne cm-2) and theMC data.This gives a value of 1543 cal/mole compared with the measured value 20of 1564 at 857°K. The discontinuity of the pressure at 86.3”K and I/ = 27.50 cm46 THERMODYNAMIC PROPERTIES OF LIQUID ARGONindicates a phase transition in the MC model ; in a real system the phase equilibriumis characterized by the equality of the Gibbs free energy and of P at different densities.Properties relating to the liquid-solid transition calculated from MC data and theTABLE 3.-LATENT HEAT OF FUSION AND VOLUMES OF MELTING SOLIDAND FREEZING LIQUID FOR ARGON AT T = 86~3°Kmelting pressure 19 1.03 x 108 dyne cm-2Monte Carlo expt.Vsol (cm3) 25.0 24,619Kiq (c&) 27.9 27.919AHfu, (cal/mole) 246 284 (T = 84"K)zoexperimental melting pressure 19 are shown in table 3.excellent agreement with experiment but VsOl is markedly too high.AHf,, is-by about 14 %-lower than the observed value.20AHf,, are insensitive to the choice of pressure.The calculated VIis is inAs a consequenceBoth AHvap andHEAT CAPACITYThe computed and observed 14 values of Ct are in good agreement at 136.0 and148.2"K.The agreement, however, becomes significantly poorer as the volumeincreases, i.e., as the critical density is approached. The MC method is expectedto underestimate near the critical point because the periodic boundary conditionsuppresses all large density fluctuations. At 86~3°K the MC values for Ci/R showa clear maximum in the neighbourhood of the pressure discontinuity (V = 27.50 ~1113).Experimental and calculated values for the orthobaric liquid are compared in table 2.Ref.(16) tabulates Cy as a function of P at 90, 100 and 110°K but the data show asmuch scatter as those given for 97.0 and 108.2"K in table 1 and interpolation isdifficult. The values of Ct/R measured at 100°K between P = 4 x 106 and 300 x 106dyne cm-2 lie between 1-09 and 1.78, suggesting that the MC results are too low.OTHER PROPERTIESThe compressibilities (column 10 of table 1) are calculated by least-squares fittingof a parabola to the MC P-V data for each temperature. The thermal pressurecoefficients (column 1 1) are determined by isochoric reweighting (eqn. (6)), AT = 1°K).The expansivity (column 12) is the negative product of the two quantities. The goodagreement at 97-0°K is probably fortuitous.Elsewhere the agreement is mostly fair.ISOTHERMAL EXTRAPOLATIONValues of 2 have not yet been reliably determined over a range of densities. Theresults obtained for V = 28.03 and 28- 95 cm3 from the MC data for V = 28.48 cm3and T = 86.3"K by means of eqn. (15) and (16) are shown in table 4. The differencesbetween the extrapolated and the independently calculated values of (@/iV) and(Y/N} are smaller than the statistical error of the MC calculation. For moderatelysinall extrapolations of (<DIN) (AV/Y-2 %) values of ," are not required. This isalso true for small extrapolations of (!PIN) (AV/V-0-5 %) when the second term ineqn. (16) can be neglected. Extrapolation from 28.48 to 28-38 cm3 yielded a com-pressibility of 2.24 x 10-10 cm2 dyne-1 compared with 2.34 x 10-10 cm2 dyne-1obtained from the MC isotherm. These results are regarded as encouragingI .R . MCDONALD A N D K . SlNGER 47TABLE d.-ISOTHERMAL EXPRAPOLATION OF MONTE CARL0 DATA FORV = 28-48 cm3, T = 86~3°Kdirect calculation extrapolationV <;> <;> @ (;)6cm3) (erg x 1015) (ergx 1015)28.03 - 10.12 2.5 1 - 10.08 2-6428.95 - 9.76 4-18 - 9.81 4.34CONCLUSIONSCalculations based on an MC method and an LJ-potential have been shown toyield quantitatively significant values of thermodynamic propel ties of a simple liquid(argon) over a considerable V,T range, including phase transitions. Satisfactoryresults have been obtained for the internal energy. For the pressure, the agreementwith experimental values is less good but still, on the whole, satisfactory.For thederivatives of UP and P the agreement with experiment is better than qualitative.In a number of cases the discrepancies between calculated and experimentalresults are larger than the statistical sampling errors inherent in the MC method.It is hoped that future work will show to what extent the discrepancies are due to theapproximations of the MC method, e.g., small sample size and the periodic boundarycondition, or to inadequacy of the potential function.By calculating not merely the mean values but also the frequency distributions forthe potential energy and the virial, it has been possible greatly to extend the T-rangecovered by a given number of MC calculations. Further economies in computationbecome possible by means of an isothermal extrapolation procedure for which onlylimited results have so far been obtained.This method will also permit the systematicinvestigation of the effects of changes of the LJ-parameters on the thermodynamicproperties.The computing time required to obtain results of the type reported in this paper islarge but not prohibitive. MC calculations for more complex systems, e.g., mixturesof simple fluids or systems governed by non-central forces, need no longer be regardedas impractical.APPENDIXAPPROXIMATIONSEFFECT OF DISTANT PARTICLES. TO reduce computing labour the sum Of all pairenergies is replaced by one in which the potential due to a uniform density is substitutedfor the effect of the distant particles,4 i.e.,@ = 4(rij)+(A7/2u)[m +(r)4nr2dr = @*(r,).In this work r , = 6.5 A.MC runs yield ensemble averages of the “effective”potential energy ( W ( r , = 6.5 A)). These values are further corrected by computingthe values of @*(r, = 6-5 A) and 4i)*(rz = 8-4 A) for a small number of configurationsat each density and adding the mean difference of (@*(rm = 8.4 A) - @*(rm = 6.5 A)]to (@*(rm = 6.5A)) in all cases. The same procedure is applied to the virials. Forthe potential energy the “effective” long-range potential amounts to 15-20 % of thetotal; for the virial, the long-range contribution is of opposite sign and of the sameorder of magnitude as the exactly calculated short range contribution. The “effective”rij <rm r48 THERMODYNAMIC PROPERTIES OF LIQUID ARGONlong range contributions and their corrections vary smoothly with density. It isestimated that the error resulting from this procedure does not exceed 0.5 % for 4and 5 % for Y.THE SAMPLE SIZE.In all calculations N = 108. The mean values of all macro-scopic equilibrium properties are assumed to be determined by the mean values ofthe MC sample. This involves errors because the periodic boundary conditionexcludes configurations which in reality occur, and because density fluctuations aresuppressed. Density fluctuations in small volume elements can be calculated fromthe isothermal compressibility. For N = 108 such fluctuations would not contributesignificantly to the macroscopic entropy if the microscopic subsystems are treated asindependent.This assumption is, I: ever, unjustified because adjacent subsystemsshare a large number of pair interaciiwns which must be expected to lead to a correla-tion of their fluctuations. No analysis of these errors has been carried out so far.Properties calculated for some u,T-points with samples of size N = 32 and N = 108are in good agreement.18THE COMPUTER PROGRAMMEA large part of the computation is done by means of integer arithmetic (i.e.,24-bit words). This reduces the execution times of transfers, additions and sub-tractions and makes an extensive use of tables possible. The sample cube has thelength L = 1000; the 3 x 108 coordinates are integers between 0 and 999. At theprevailing densities a hard core radius re = 2.7 A corresponds to > 100 units. Noserious distortions are likely to result from the discreteness of the variables. Changesof volume are taken into account by scaling of the LJ parameters. The periodicboundary condition requires that in the calculation of a pair distance rfj2, the com-ponents x ~ j , yij, zfj must be “modified”, i.e., if xi-.q>L/2, it is replaced byxg +L- xj, etc.To determine the interaction potential between particles i and j ,the computer calculates xi - xj, yi -J)], zi - zj from the listed particle coordinates ;it then obtains the squares of the “modified” differences from a table, adds these andlooks up the values of the potential and of the virial functions for the rig2 from tables.Lists are kept of all coordinates, of all pair potentials in the “present” configuration(Wfj) and of the contribution of each particle i to the “present” value of iD (Ai =+cW& ‘&4g = a).When a move is attempted by a displacement of particle k,a list is formed of the provisionally altered pair potentials (Vkj, all ji:k), and theprovisional new contribution of k to @ ( A = Eva,). The provisional new @ isaccepted or rejected according to the criterion of ref. (1) (the exponential functionsrequired for exp (-PA@) are tabulated). If the move is accepted, the lists of co-ordinates, pair interactions, particle contributions and the value of @ are updated(e.g., V k p Wkj, A+&) and an entry corresponding to the new value of @ is made inthe histogram. If the move is rejected an entry corresponding to the “present” valueof @ is made (again) in the histogram.The treatment of the virials by means oftables is analogous.The programme generates about 3000 configurations per minute on the Mark 1Atlas computer of the University of London. Calculations for a v,T-point werebased on not less than six ten-minute runs (168,000 configurations) of which the firstone was discarded to allow for “equilibration”. The formation of tables for thepotential and the virial functions at the beginning of each run does not add signifi-cantly to the computing time; calculations with any other type of pair potentialj4i ij 4 I . R . MCDONALD AND K . SINGER 49instead of the LJ-potential would not present any problem (though eqn.(13)-(16) forthe isothermal extrapolation would not be applicable).We thank the Institute of Computer Science (University of London) for theallocation of sufficient computer time, Mr. K. R. Knight for his contribution to thecomputer programme, and the Science Research Council for financial support. Weare grateful to Prof. J. S. Rowlinson for bringing to our attention some importantexperimental data.1 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J, Chem.2 M. N. Rosenbluth and A. W. Rosenbluth, J. Chem. Physics, 1954, 22, 881.3 W. W. Wood and J. D. Jacobson, J. Chem. Physics, 1957,27,1207.4 W. W. Wood and F. R. Parker, J. Chem. Physics, 1957, 27, 720.5 W. W. Wood, F. R. Parker and J. D. Jacobson, Nuuuo Cimento, supp. 1, vol. 9, series X, 1958,6 Z . L. Salsburg, J. D. Jacobson, W. Fickett and W. W. Wood, J. Chem. Physics, 1959, 30, 65.7 E. B. Smith and K. R. Lea, Trans. Faraudy Soc., 1963,59, 153.8 A. Rotenberg, J. Chem. Physics, 1965, 43, 1198.9 A. Rotenberg, J. Chem. Physics, 1965,43,4377.1oB. J. Alder, Physic. Rev. Letters, 1966, 16, 88.11 A. A. Barker, Austral. J. Physics, 1965, 18, 119.12 S. G. Brush, H. L. Sahlin and E. Teller, J. Chem. Physics, 1966, 45, 2102.13 A. Michels, Hub. Wijker and Hk. Wijker, Physica, 1949, 15, 627.14 A. Michels, J. M. Levelt and G. J. Wolkers, Physica, 1958, 24, 769.15 A. van Itterbeek and 0. Verbeke, Physicu, 1960,26,931.16 A. van Itterbeek, 0. Verbeke and K. Staes, Physica, 1963, 29, 742.17 W. W. Wood and J. D. Jacobson, Proc. Western Joint Computer Conference, (San Francisco,18 I. R. McDonald and K. Singer, to be published.19 F. Din, Thermodynamic Functions of Gases, vol. 2, (Butterworths, London, 1956), p. 146-201.20 P. Flubacher, A. J. Leadbetter and J. A. Morrison, Proc. Physic. SOC., 1961, 78, 1449.21 J. S. Rowlinson, Liquids and Liquid Mixtures, (Butterworths, London, 1959), p. 64.Physics, 1953, 21, 1087.133.1959), 261

ISSN:0366-9033    

DOI:10.1039/DF9674300040

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

GENERAL DISCUSSIONDr. J. A. Barker and Dr. D. Henderson * (C.S.I.R.O., Melbourne) said: At present,reliable values for the radial distribution function and the equation of state are avail-able only for the hard-sphere fluid. Thus, a promising method for dealing withsystems of molecules with attractive forces is to consider the attractive potential asa perturbation on the hard-sphere potential. We have recently used perturbationtheory to calculate the equation of state of a fluid of molecules interacting accordingto the square-well potential. This is a particularly good potential because the effectof the attractive forces is not complicated by the “ softness ” of the repulsive partof the potential. In addition, Monte Carlo and molecular dynamics calculationshave provided quasi-experimental data with which our calculations can be comparedwithout the uncertainty due to presence of three-body forces and the lack of know-ledge of the intermolecular potential which is inevitable in applications to realfluids.We divide the intermolecular potential into the “ unperturbed ” hard-spherepotential uO(R) and the “ perturbation ” u,(R).In our calculations, u,(R) = --8for 0 < R < 1-50 and is zero otherwise. If N , is the number of intermolecular distancesin the interval a<R< I.%, then the Helmholtz free energy can be written :where Fo is the free energy of the unperturbed system and the angular brackets mean“ average over the configurations of the unperturbed system ”.The average ( N , ) is given by an integral over the radial distribution function ofthc unperturbed system, go(R).Thus, the first-order term is identical with that givenby Z ~ a n z i g . ~ The second-order term is equivalent to that given by Zwanzig but issimpler and more suggestive. For example, it is intuitively plausible that the devia-tions of N1 from its mean value should be least important at high densities so thatthe perturbation expansion (1) converges best at high densities. Also N1 can beregarded as representing the number of molecules in a spherical shell surrounding acentral molecule. If this shell were a large macroscopic volume, then<m - <w2 = <N,>kT(aP/aP), (2)where (dpldp) is the macroscopic compressibility. An even better approximationis obtained if we use the local compressibility instead of the macroscopic compressi-bility, i.e., if we replace (dp/dp)ogo(R) by (d/dp)[pg(R)],.Our results are based onthis local compressibility approximation. However, similar results are obtained ifthe macroscopic compressibility approximation is used. We are presently makingMonte Carlo calculations of ( N f ) and ( N , ) which will eventually supersede (2).At present we have only a few scattered results. However, these results do indicatethat (2) is a reasonable approximation.* permanent address : University of Waterloo, Canada.’ B. J. Alder, unpublished results.A. Rotenberg, J. Chern. Physics, 1965, 43, 1148.R. W. Zwanzig, J. Chern. Physics, 1954, 22, 1420.5GENERAL DISCUSSlONUsing (2), the free energy is51where q = npo3 /6.In obtaining (3) we have used the Percus-Yevick (PY) compressi-bility isotherm to obtain (aplap),. In evaluating F we used the PY expressionfor go(R). To order T-2, (3) gives the exact second virial coefficient.8-40 0.2 0.4 0 0Vol vFIG. 1 .--Equation of state for the square-well potential. The points are the Monte Carlo values ofRotenberg and the curves are isotherms calculated from eqn. (3). The points given by 0 0 , $ and 0were calculated using 256 molecules at E/kT = 0,0-3, 1 and 2, respectively, whife the points given by0 were calculated using 864 molecules at E/kT = 1.In fig. 1 we have compared the equation of state which results from numericaldifferentiation of (3) with the Monte Carlo calculations of Rotenberg.' The agree-ment is good-even at the lowest temperature, kT/& = 0.5, which is far below thecritical temperature.The agreement with Alder's molecular dynamics calculationsis even better.We have generalized this procedure to include potentials with a soft core suchas the 6 : 12 potential. For an arbitrary potential u(R) we define a modified potentialM. S. Wertheim, Physic. Rev. Letters, 1963,10, 321 ; J. Math. Physics, 1964,5, 643.J. Chem. Physics, 1963, 39, 474.A. Rotenberg, J . Chem. Physics, 1965, 43, 1148.B. J. Alder, unpublished results.E. Thiele52 GENERAL DISCUSSIONby the equations, ( " i d ) , d + - - a , R - dt'(a,y,d,Q; R) = u d+-UR - d Q - i (= O , a<d+-<d+-,U &= yu(R), R > Q. (4)When cc = y = 0, v becomes the hard-sphere potential with diameter d, while foru = y = 1 the original potential u(R) is recovered ; u is an inverse steepness parameter1 I I I I I0.0 0.2 0.4 0.6 0,B 1.0No '1 VFIG.2.-Equation of state for the 6 : 12 potential, kT/c = 2-74. The solid line gives the resultsobtained from eqn. (5) and the broken line gives the results. of McQuarrie and Katz. The pointsgiven by 0 and 0 were calculated by Wood and Parker for 32 and 108 molxules, respectively.for the repulsive region and y a depth parameter for the attractive region. By amethod based on the ideas of Rowlinson and Zwanzig the first derivatives of theconfiguration integral with respect to u and y can be calculated in the limit u, y+O.Thus, we find for the Helmholtz free energyF, -?go( y, ,I)[ JI exp (- u( R)/kT)dR -(a - d ) +z7Jom go ( Nd3 7'2 ' r m R 2 d R kT + higher order terms, ( 5 )1 F -- -NkT - NkTJ.S. Rowlinson, Mol. Physics, 1964, 7, 349; 1964, 8, 107.R. W. Zwanzig, J. Chem. Physics, 1954, 22, 1420GENERAL DISCUSSION 53where Fo and go are the free energy and radial distribution function respectively forthe hard-sphere system. For a = y = 1 we obtain an expression for the free energycorresponding to the potential u(R). We choose CT so that u(a) = 0 and determine dso that the second term in (5) is zero. These choices should minimize the contributionof the higher-order terms and are a generalization of Rowlinson's pr0cedure.l Inaddition, these choices ensure that (5) reduces to (3) if the square-well potential is used.Thus, d is an effective core diameter which may depend on temperature but not ondensity. Pressures calculated by differentiation of (5), using the Percus-Yevickradial distribution function, the 6 : 12 potential, and neglecting the higher-orderterms, are compared with the Monte Carlo results of Wood and Parker for~ T / E = 2.74 in fig.2. The contribution of the term in y 2 calculated by our " localcompressibility " approximation is very small at this temperature. Also shown areresults given by the equation of McQuarrie and K a t ~ . ~ The results given by eqn. (5)appear to be satisfactory.Dr. P. Hutchinson (A.E.R.E., Harwell) said: I report the results of some calculationsin which the Percus-Yevick (P.Y.) and hypernetted chain (H.N.C.) approximationsare tested against the molecular dynamics calculations of Rahman at densitiescorresponding to the liquid phase of argon.From thepair correlation function h(r), computed by Rahman, the direct correlation functionC(r) was calculated. It is then possible from the P.Y. and H.N.C. equations tocalculate the potentials 4p.y. and #H.N.C. which would give such forms of C(r) andh(r) in these approximations. These can then be compared with the potential 4(r>which was the starting point of Rahman's calculation. This is a 6-exp Buckinghampotential truncated at 7.65 A ; the temperature was T = 85.5"K, and the densitywas 0.02103 atm/A3 (equivalent to an argon density of 1.407 g/cm3). The temperatureis thus significantly below the critical temperature for argon.The numerical procedure is to calculate C(r) via a double Fourier transform.Fourier transformation of the Ornstein Zernike equation giveswhereThe procedure is that suggested by Johnson, Hutchinson and March.5a k ) = L(k)/(l +pJi(k)), (1)(14 h{k) = J exp (-ik .r)h(r)d3r.Hence to calculate C(r) we first evaluate h(k) by Fourier transformation. Fromthis we may obtain C(k) from (I), and c ( k ) may be transformed to obtain C(r).Great care is necessary to obtain accurate results from a numerical Fourier transformsmethod. The overall accuracy of the method may be assessed by evaluatingh,(r) = - [exp (ik . r)h"(k)d3k.Formally, h,(r) = lz(r), so that the difference between the two functions is an indicationof the error in the method.Comparison of these two functions showed that thedifferences were about 1 % except at the origin where they may rise to -5 %. How-ever, errors near the origin are unimportant as this region contributes little to the(W3J. S. Rowlinson, Mol. Physics, 1964, 7, 349 ; 1964, 8, 107.W. W. Wood and F. R. Parker, J. Chern. Physics, 1957, 27, 720.D. A. McQuarrie and J. L. Katz, J. Chem. Physics, 1966, 44, 2393.A. Rahman, Physic. Rev. Letters, 1964, 12, 575.M. D. Johnson, P. Hutchinson and N. H. March, Proc. Roy. Soc. A , 1964,282,28354 GENERAL DISCUSSIONtransformation (la). An additional check is available by direct solution of theOrnstein Zernike equation for C(r). Given h(r) the 0-Z equation is simply aninhomogeneous linear integral equation which may be solved by direct inversion.While this last method is subject to greater numerical error it produces the samegeneral result for C(r) as that described below.The results of the Fourier transform method are shown in fig.1 where C(r) isplotted with h(r) and - 4(r)/kT. The values for C(r) are reliable for r <4 A, in thatthe values obtained there are insensitive to the truncation point of h(r) and the choice1.7 2.0 3.0 4.0 5 . 0 6 . 0 7.0 &O 8.5r, (A)- , &I; - - -Y C(4; 0, cp(r)/kT.Fro. i.-Comparison of the direct and total correlation functions and the potential.of mesh for the integral over k. Beyond 6A the results are strongly sensitive to theabove variations and decrease in magnitude when C(k) is smoothed.This indicatesthat the last part of the curve is due to cancellation error (noise) and that C(r) isvanishingly small in this region. The most notable feature of C(r) is that it has thesame general shape as the potential but dies away rather more rapidly. This iscontrary to the expectation of both the P.Y. and H.N.C. approximations thatC(r)- -4(r)/kT for r large. In fig. 2 are shown the two approximate potentials4p.y.(r) and $H.N.C.(r) in comparison with &), given in dimensionless units. The fitobtained by 4H.N.C.(r) is much the better of the two. The well depth of q5p.y.(r) is muchtoo large and displaced too far to the right. The comparison between well depths is+p,ye(4.6)/kT = -4.01, #H.N.C.(4*0)/kT = - 1.28, +(3*8)/kT = - 1-40GENERAL DISCUSSION 55It has been shown by Gaskell that when the compressibility is small the H.N.C.approximation gives values of the pressure which are several orders of magnitude toolarge the virial equation is used. This situation would hold here as calculation giveskT(ap/ap), N 0.05.2.0I .51.0k5BI 0 ' 5nW0.0- 0 5- 1 .01 :0i ' iI Ir, (4FIG. 2.Xomparison of exact and approximate potentials.-? -rg(r)/kT; - - -$ -WXNc(r)/kT; 0, -cpF.Y,(r)/kTThe virial equation is2n O0 P = pkT - ,p2 [ r34'(r)g(r)dr.J oThis minimum value of $H,N.C,(~) is shifted just far enough to bring the large positivevalues of -@(Y) into coincidence with the peak in g(r), and this produces the grossover-estimate of the pressure.The virial equation is thus a very sensitive test ofany theory in the low-pressure, high-density state of a liquid. Finally we may con-clude that despite the poor results obtained from the virial equation the H.N.C.approximation is superior to the P.Y. equation for a liquid well below its criticaltemperature. Grateful thanks are due to Dr. A. Rahman for making his resultsavailable for this calculation.Prof. J. S . Rowlinson (Imperial College, London) said: In his lecture Rushbrookementioned that he and Silbert have been able to show that the inclusion of a tripletT. Gaskell, Proc. Physic. SOC., 1966, 89, 231.G. S. Rushbrooke and M. Silbert, Mol. Physics, 1967, 12, 50556 GENERAL DISCUSSIONpotential in the hyper-netted chain theory is " particularly simple " since it leads onlyto the introduction of a new elementary graph, x1(1,2) in his notation, which is oflower order in the density than those arising from the two-body potential.There isan analogous result for the Percus-Yevick the0ry.l The new graph is here alwaysof the form [e(l,2)xl(l,2)], which leads at once to the physically reasonable resultthat the triplet forces are of importance only outside the repulsive regions of the two-body potential, i.e., outside the regions in which e(1,2) = 0. Moreover, the simplegraphical expansion that is so characteristic of the Percus-Yevick theory is stillpreserved when the triplet potential is included. It can be shown that g(1,2) is stilla sum of all graphs in which the interiorf(i,j) bonds do not cross.The only changeneeded is the inclusion of those graphs in which one or more of the peripheral bondsare notf(i,j), but [e(i,j)xl(i,j)].Dr. S. Levine (Manchester University) (partly communicated) : The Percus-Yevicktheory has been described in terms of " turning on " an external field which enablesone to vary arbitrarily the local number density of molecules in the fluid in the grandcanonical ensemble. This density is regarded as an independent variable and thena suitable dependent variable must be so chosen as to give a rapidly convergentfunctional Taylor expansion. Are there not other and indeed simpler and moremeaningful explanations giving the physical basis of this theory ?Prof. J. S. Rowlinson (Imperial College, London) said: I shall attempt to answerthis question in my summarizing remarks.Mr.A. Moreton and Prof. G. H. A. Cole (University of Hull) (communicated):We too would stress the importance of the Born-Green-Yvon equation as an exactequation which shows a fluid-solid transition. The need to improve upon theKirkwood closure procedure is then of central importance in the theory of liquids.We have recently used the procedure of functional Taylor expansion to express g ( 3 )in terms of gt2) and the pair potential function u(r). The general method of functionalexpansion has been set down by Verlet amongst others.For our purposes we expand the functional p(l)(3/u12) [exp (Pu,(3)) - 11 as a seriesin terms of p(')(4/u1)[exp (Pu1(4)) - 11. Here p ( l ) is the single particle distribution,uI2 is an external potential due to the presence of two extra particles at positions rland r2.The expansion is made about the condition u1 = 0; u12 = u1 +u2 reducesto u12 = u2 in this limit. Taken to the linear term, our expansion gives the Kirkwoodsuperposition statement, which is (1) of the Rice-Young paper but with z = 0.The inclusion of the quadratic term in the expansion gives9'3'(1,2,3) = g'2'(1,2)g'2'(l,3)g'2'(2,3) 1 + p dr4y(1,4)f(1,4) x 1 1Here we useY(r> = g'2'(r) exp p 4 9 ,f ( r ) = exp (- pu(r)) - 1.In the law-density limit the various functions in (1) can be expressed as a simpledensity expansion using an iterative method and (1) is then exact up to the linear termJ. S. Rowlinson, MoZ.Physics, 1967, 12, 513.L. Verlet, Physics, 1966, 32, 304GENERAL DISCUSSION 57in the density. When used in conjunction with the Born-Green-Yvon equation,(1) then provides as an exact expression for the fourth virial coefficient. The tripletterm in the functional Taylor expansion has also been obtained. Although it is toolengthy an expression to set down here we have found that, in the low density limit,it leads to an exact expression for g ( 3 ) up to the quadratic density term, and so givesan exact expression for the fifth virial coefficient.For higher densities than those where g ( 2 ) can be expressed as a simple densityexpansion, the above expansion procedure is still applicable and the expression (1)for g ( 3 ) can be used as a closure expression for the Born-Green-Yvon equation.The solution of this more general situation, for given u(r), could apply to higher fluiddensities. We are at present calculating g(2) by using (1) and the Born-Green-Yvon equation, for a function u of the Lennard-Jones type but with a hard core atvanishing interparticle separation.Our method could be extended to includecorrelations higher than the triplet. In this way it can be considered with the methodof Fisher.'Prof. Stuart A. Rice and Dr. David A. Young (University of Chicago) said: Themediocre quality of the thermodynamic functions which are obtained from thesolutions of the Yvon-Born-Green equation seem to exclude this equation as a usefultool in the study of fluids. However, by modifying the YBG equation we mayobtain useful and interesting results.The YBG equation arises from the introductionof the Kirkwood superposition approximation into the second member of an infinitesequence of integral equations. This latter equation, which involves the three-particle distribution g( 129, is exact, and we may construct increasingly good approxi-mations to this exact equation by multiplying the superposition product by a seriesof triplet correlation terms. By introducing the first two correlation terms, we getgood results for rigid spheres and discs, and for 6-12 and square-well potentials.We may even expect good results from the simple YBG equation for long-rangepotentials of the Coulomb type.2However,Kirkwood has shown that periodic, solid-like solutions may exist at densities abovethis singular point.Also, by introducing an approximation to the long-rangebehaviour of g(123) into the modified equations we may be able to correlate theYBG singularity with the fluid-solid phase transitions predicted by molecular dynamics.Despite the limitations imposed by the speed of present-day computers, we hope thatthese studies will add a little to our understandirlg of the statistical theory of fluidsand phase transitions.Prof. J. Walkley and Dr. Wing Y. Ng (Simon Fraser University) said: In twopapers use has been made of the Lennard-Jones 12-6 potential. In each case anagreement between theory and experiment is found, despite the inadequacy of thispotential function to express the two-body second virial coefficient behaviour.Inan N-body situation the total potential experienced by any one particle cannot beregarded as the simple pair-additive sum of the interaction with all other particles.For the 3-body case the perturbation of the pair-additive s u n coming from the mutualdisposition of the three particles has been calculated for the dispersion r6 attractionterm.4 The use of the 12 : 6 pair potential in a strictly pair-additive manner appearsto generate an (apparently) acceptable N-body potential. In calculations pertainingI . Z . Fisher, Statistical Theory ofLiquids (Univ. of Chicago Press, 1964), p. 152.C. W. Hirt, Physics Fluids, 1967, 10, 565.J. G. Kirkwood and E. Monroe, J. Chern. Physics, 1941,9, 514.R. J. Bell and A.E. Kingston, Proc. Physic. Soc., 1966, 88,901.The exact meaning of the singularity in the YBG equation is unclear58 GENERAL DISCUSSIONto the solid phase of the inert gases, it was similarly observed that the simple bi-reciprocal 12 : 6 potential allowed theoretical volume-temperature and heat capacity-temperature dataare in good agreement.l For arange of two parameterm : 6 potentials(viz., 4(rl1) = &f(rij/u), with usual rotation) the characteristic parameters E and Q canbe obtained by fitting the theory (a harmonic Einstein theory) to certain propertiesextrapolated to OK. The total lattice energy U, and the molar volume V , wereoriginally used. The " best '' set of parameters were then obtained by a plot of thetheoretical reduced zero-point energy (A:) and the reduced " experimental " zero-point energy (9R0, /8N&, with 0, the limiting Debye frequency) simultaneously asfunctions of the de Boer parameter (A = h/(m&)%o).The intersection of the twocurves give the " best " value of E and Q without specifying the " best " potentialfunction. However, this set of values allowed both second virial and viscosity datafor the inert gases to be reduced to a corresponding states pattern of behaviour.2TABLE 1.-nz : n KIHARA CORE POTENTIAL PARAMETERS FOR ARGONV 20 .(A1 B (cal mole-1) m : n E/kCK)9 : 6 123.2 3.3932 0.157 0 18710 : 6 123.0 3.3948 0.1056 0 18711 : 6 122.0 3.3960 0.0567 0 18712 : 6 120.5 3.3987 0.0098 0 187Experimental data : Uo = - 1846 cal mole-' ; Ro = 3.7549 (A) ;KT = 0-375 cm dyne-2 x lolo ; Lo = = (9/80,) = 187 cal mole-'.TABLE 2.-m : n KIHARA CORE POTENTIAL PARAMETERS FOR XENONB (cal mole-1) V A0 m : n ~lk("K) d A )9 : 6 240 3.9483 0.157 0 13310 : 6 238 3.9510 0.106 0 13311 : e 236 3.9521 0.0567 0 13312 : 6 234 3.9599 0-0098 0 133Experimental data : Uo = - 3856 cal mole-' ; Ro = 4.3356 (A) ;KT = 0-28 cm dyne-2 x 10'' ; = (9180,) = 123 cal mole-'.We have extended these calculations to include a non-additivity term in the staticlattice potential sum.In keeping with theory we consider only the largest termcoming from the 3-body interaction and written as Av/a$, where A is the appropriatelattice constant, a i j the nearest neighbour distance and v an unknown non-additivecoefficient. Calculations were also made more flexible by considering a series ofrn : 6 Kihara core pair potentials, viz.,where r = rlj/a tc.The calculation involves four unknown coefficients, the reductionparameters e and a, the coefficient v and the reduced hard-core cut off 6. Fourzero-point properties are required for their evaluation, and besides the three usedabove, U,, Yo and Ro, an accurate value of the zero-point isothermal compressibilityKT could be used for both argon and xenon. The parameters obtained by fittingthe theoretical model (again a harmonic Einstein model) to this data are given intables 1 and 2.I. H. Hillier and J. Walkley, J. Chem. Physics, 1965, 43, 3713.J. Walkley, J . Chem. Physics, 1966,44, 2417GENERAL DISCUSSION 59From these tables the following observations may be made: (i) for all m : 6potentials examined the non-additivity coefficient is effectively zero ; (ii) for bothargon and zenon, for the same value of rn the resulting value of 6 is identical, inti-mating an essentially corresponding states behaviour pattern, and (iii) for both argonand zenon each of the rn : 6 Kihara potentials gives, through the different 6 values,curves which are virtually identical and which are exceedingly close to that of a 12 : 6Lennard-Jones potential.The calculations are, then, consistent with the accepteduse of the Lennard-Jones 12-6 potential as an acceptable '' many body " potential.Prof. U. v. Weber (Rostock University) said: We studied the deviations fromideality of vapours by an isochoric method to get the second virial coefficients andby means of them the intermolecular potentials.Density was varied to take careof higher coefficients. To make the adjustment of parameters and powers in runambigously, the range of temperature must be as large as possible. The bestfitting powers of attraction obtained by the method of Buckingham are r9 forbenzene and n-octane and r-l0 for hexane and n-octane with reasonable diameters.I restrict myself to quasispherical molecules.Only r6 is theoretically justified and is indeed a fair potential for rare gases, butnot satisfying for large molecules. In our Institute, W. Lichtenstein investigatedvarious potential models. Kihara allows for extension to hard cores, but his assump-tion of effect only between the two nearest points is not justified. Laeuger consideredr6 potentials between the corners of two cubes. Hoover and Rocco regardedspherical shell potentials with centres of attraction and repulsion on the shell, butthey integrated only numerically. Lichtenstein considers two spherical models :attraction in -6 power between every volume element of the first particle to everyvolume element of the second and between every surface element of the first particleto the second. Repulsion was assumed in either model between the two centres inpowers -24, to obtain an analytically integrable expression. The shell model fittedbetter. The argon, tetrachlormethane and benzene data were fitted as well asdesirable by this potential with reasonable diameters. To avoid singularity in the caseof contact, recession of the points of attraction up to a distinct amount was supposed.This seems to me a reasonable step towards a dynamic model. The potential well ismuch deeper and more narrow than for a Lennard-Jones potential, a fact whichmight be taken into account in the more complicated theory of liquids

ISSN:0366-9033    

DOI:10.1039/DF9674300050

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Random Close-Packed Hard-Sphere ModelI. Effect of Introducing HolesBY J. D. BERNAL AND S . V. KINGBirkbeck College, University of London, Malet Street,London, W.C. 1.Received 19th December, 1966The structures have been examined of lower density random models, obtained by removing atomsfrom the close-packed, random-model built here. Radial distribution functions and coordinationdistributions for models containing different proportions of holes have been calculated.Radial distribution functions calculated from X-ray or neutron diffraction patternsof liquids show that the effect of temperature is less to increase the interatomicdistances between atoms than to reduce the number of neighbours of each atom.That is, the position of the first peak in the radial distribution function changes littlebut the total area of the peak which is proportional to the number of first neighboursdecreases.In the random close-packed sphere model of a liquid the increase in volumecorresponding to increase in temperature 1 could be effected in two ways.First, theproportions of the basic polyhedra could be changed, i.e., there would be an increasein the proportions of the larger polyhedra so that there would be more archimedeanantiprisms and tetragonal dodecahedra and fewer of the smaller polyhedra, theoctahedra and tetrahedra. Secondly, atomic-sized holes could be introduced intothe body of the liquid.By using the first method alone it is only possible to achieve a maximum increasein volume of about 17 %. It is of interest therefore to investigate the effects of intro-ducing atomic sized holes at random into the hard sphere model (or rather to removeatoms at random from the model).Atoms were removed randomly from the model.Initially, the removal of atomswas restricted so that no two holes could be neighbours, a " neighbour " being definedto be within J2 diameters of another atom. There will be a limit to the number ofholes which can be introduced in this way. A preliminary calculation was carried outto ascertain the maximum possible number of holes which could be removed fromthe model so that no two holes were neighbours.This was done as follows: the central atom in the mass was removed, then,looking at the atoms in order of their distance from the centre of the mass, each atomwas tested if it was a neighbour of a previously chosen hole. The atom remainedin the assembly if it was not a neighbour of any previously chosen hole ; then it wasremoved from the assembly.This process continued until all atoms in the assemblyhad been visited. The maximum number of holes which could be introduced in thisway was 18-6 % (a neighbour being defined as within 1.414), roughly correspondingto an increase in volume of 18.6 %. If a neighbour was defined as being within 1-1atomic diameters (position of potential = 0 for Lennard-Jones potential) the number6J . D . BERNAL AND S . V . KING 61of holes which could be introduced was 26 %. For a neighbour defined as within1.05 diameters the number of holes which could be introduced was 28-3 %.Two sets of calculations were carried out on introducing holes at random into thehard sphere assembly : (A) atoms were removed only if they were not neighbours ofdistance from centre in atomic diametersFIG.1.-Radial distribution functions. Random model with no holes and 35 % holes.coordination number (for neighbours within 1-1 diam.)FIG. 2.-Coordination distribution functions. Random model with no holes and 35 % holes.holes already introduced (i.e.? within 1.414); (B) atoms were removed at random,holes could be neighbours. In the first (A) set of calculations after about 16.6 %of holes had been introduced it was impossible to introduce any further holes unlessholes were allowed to be neighbours.In both sets of calculations, holes were introduced in increasing proportions of5 %.First 5 % holes were introduced then 10, 15, 20, 25, 30, 35,40 and finally 45, 62 RANDOM CLOSE-PACKED HARD-SPHERE MODELthis being well above the critical volume. In the first (A) set of calculations after thesaturation point of 16.6 % had been reached all further holes were allowed to beneighbours. For each of these new assemblies the radial distribution function andcoordination distribution function was calculated.In general outline, the effect of introducing holes into a random close-packedhard-sphere assembly is analogous to the effect of temperature on a real liquid. Thepeaks of the radial distribution functions become lower and broader (fig. 1). In thecoordination distribution functions the mean of the curve moves to lower values ofcoordination, and the coordination distribution function itself broadens and becomeslower. This is clearly seen by comparing the coordination distribution histogramobtained from the random model with no holes with that for the random model with35 % holes. The mode falls from 9 in the first case to 6 in the latter (fig. 2)

ISSN:0366-9033    

DOI:10.1039/DF9674300060

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

62 RANDOM CLOSE-PACKED HARD-SPHERE MODEL11. Geometry of Random Packing of Hard SpheresBY J. D. BERNAL AND J. L. FINNEYA preliminary Voronoi polyhedron analysis of a random model has enabled an accurate fixing ofthe density of random close-packing. A wide variation (-15 % of the density of closest regularpacking) in local density is found, correlated with various topological types of polyhedra. The sus-pected predominance of 5-sided faces is verified, while the average number of faces per polyhedron,and the average number of sides per face, are a few per cent greater than the theoretical predictions ofCoxeter. Much data remain still to be analyzed for both physical and computer models.The radial distribution function of a random close-packed array of hard spheresshows a striking similarity with simple liquids, as does the difference in densitybetween this packing and the corresponding crysta1.Z 3 Ball-bearing models havebeen studied by Bernal 1 and Scott,49 5 and the synthesis of such models attemptedusing computers.Complete descriptions in terms of the coordinates of the packedcentres have been obtained, providing data for a thorough examination of the arrays.Among other reasons, the random packing is of special interest because of itssimilarities with the liquid state, and the geometry of the arrangement itself. Amanageable overall description is necessary if the thermodynamics of such an array isto be computed, though some progress (e.g., in the heat of fusion) can be made usingcomputer computations directly on the coordinates, and further advances may bepossible by exploring Monte Carlo techniques.The development of some kind of“ statistical geometry ” appears to be one of the most promising lines of attack.From our first measured ball-bearing array, a large ball and spoke model wasconstructed. This showed many interesting features.1 For example, consideringvectors between neighbours and near neighbours, the whole model was seen to consistof only five different kinds of polyhedra (“ canonical holes ”) as the basic units,suggesting that perhaps these might be treated as the building bricks of a statisticalgeometry. Our present work, however, concentrates on describing the model as anetwork of Voronoi polyhedra. If we take an array of points, and perpendicularlybisect the vectors between them, we obtain a large number of intersecting planes.To define the Voronoi polyhedron of a point, we select the smallest polyhedron SOformed about that point, ensuring no further possible planes can cut this chosen set.This polyhedron contains all those points closer to the chosen centre than to any other :hence, a network of such polyhedra completely fills space.The construction of thetwo dimensional “Voronoi polygon” is shown in fig. 3. Now each polyhedroncontains information sufficient to describe completely the neighbourhood of a point J. D . BERNAL AND J . L. FINNEY 63we are trying to develop some statistical theory to describe an extended random modelin these terms.In a pilot experiment,l a number of plasticene spheres were compressed together,and the distribution of faces on each resulting polyhedron recorded.More recently,we have developed a computer programme to elucidate the polyhedra exactly usingthe coordinates of a measured hard-sphere array. Before we could use the data ofScott,5, 6 it was necessary to correct the coordinates of a few centres apparentlyseparated by much less than a sphere diameter. This was done by building a largescale ball-and-spoke model-which incidentally exhibited the same structural featuresof our original model ; no rigorous structural analysis on the basis of the five canonicalholes has, however, been attempted. All the errors thus found were simply ones oftranscription. This large model also proved a help in choosing those centres whoseVoronoi polyhedra would not intersect the outer boundary : in this way 407 of the1006 centres were chosen.FIG.3.-The construction of a Voronoi polygon.The Voronoi computer programme starts by calculatingpolyhedron. As every vertex is equidistant from four centres,the vertices of eachthis basically reducesto solving four siinultaneous equations. All possible combinations of four centreswithin a restricting combination radius rc are taken, and the resulting solution acceptedas a vertex if no other centre is closer than the initial four. The value of rc was fixedat 1.6 diameters by trial runs on a limited number of polyhedra : if this proved in-sufficient for any particular polyhedron (as it did in two cases where rc was con-sequently increased to 1.61d) an inconsistency would arise in the form of a missingvertex.The reference numbers of the four centres associated with each vertex arestored : this provides sufficient data for the second half of the programme to sort outthe topology. Each polyhedron takes about two seconds on the University ofLondon’s Atlas computer.The programme produces the following data for each polyhedron : (1) the centrenumber and its coordinates ; (2) the coordinates of the vertices, with their associatedcentres ; (3) for each face : (a) the number of edges, (b) the area of the face, (c) thevolume subtended by each face at the centre (d) the centre-face perpendicular distance,(e) the number of the centre with which the face is shared; (4) the total number offaces, total area, total volume, and percentage density ; (5) the number of faces with3, 4, .. . etc. edges. e.g., 03651000 denotes three quadrangular, six pentagonal, fivehexagonal, and one septagonal faces.Cumulative histograms are output of: (1) the number of M-sided faces ; (2) thenumber of N-faced polyhedra ; (3) polyhedron density64 RANDOM CLOSE-PACKED HARD-SPHERE MODELAll these data are stored on tape, and are easily accessible for further computationif required.SURVEY OF INITIAL RESULTSDENSITY OF RANDOM PACKINGScott 4 has measured the density directly, allowing for surface effects by extra-polation. We have previously attempted a measurement by calculating the actualvolume of material within spheres of increasing radius from the centre of the array :this also involved extrapolation as the oscillations had not quite died out when thesphere cut the mass boundary.In contrast, the Voronoi polyhedron provides, within the limits of the data, anexact measurement of average density, as well as the local variations.The resultsfor 407 polyhedra are shown in fig. 4, which is a histogram of polyhedron density(percentage occupied by the sphere) in intervals of 0.25 % : the large spread from57 to 70 % is obvious. The overall average is 63.42 %, which occurs at a minimum% densityFIG. 4.-Polyhedron densities.in the histogram. If we plot separate histograms for polyhedra with the same numberN of faces ( N = 12, 13, . . ., 19), we find, as expected, that the mean density p increasesas the number of faces N falls, with considerable overlap between them. The resultsof these separate histograms are plotted in fig.5. The lengths of the ordinates indicatestandard deviations from the mean (except for polyhedra with 18 and 19 faces),and the figure by each point is the number of such N-faced polyhedra. The averagefor the whole 407 centres is also plotted for comparison.Not only the average density, but the peaks also appear to have no obviousstructural significance (e.g., in terms of polyhedron type). The first implication is thatthe average density is a purely statistical figure. But this does not invalidate thepossibility of there being a statistical upper limit to the overall density of randompacking, as we cannot consider the array solely in terms of isdated polyhedra.Forexample, although some very high density polyhedra occur, we cannot necessarilyincrease the proportion of these at the expense of those of lower density, as thegeometry of the packing may not allow this. Thus, although there may be localregions of density greater than the average they may occur only in conjunction withbalancing polyhedra of lower density such as to keep the average density below J . D . BERNAL A N D J . L . FINNEY 65maximum value. A study of the packing of polyhedra around a central polyhedronis an obvious next step. We are also planning to examine the changing local densitypattern in packings of different average densities, and see how the distribution of+I201 i-Ii LI II 12 13 14 I\ Ik 17 I$ I91 1 I 1number of faces per polyhedronFIG.5.-Polyhedron density against number of faces.polyhedron types, as well as their sizes, changes. A detailed comparison of ourresults with the density curve of Kiang,7 obtained from a completely random fragmenta-tion of space, may also be useful, although our system is much more restricted.FACES AND EDGES OF POLYHEDRABoth the average number N of faces per polyhedron, and the average number @of edges per face may be important in the overall statistical picture. The distributionof M and N are shown in figs. 6 and 7, which are histograms of number of polyhedrawith a given number N of faces, and number of faces with a given number M ofedges respectively.The mean N of 14.28 is considerably higher than the 13.6 obtainedfrom the plasticene spheres 1 and the 13.56 from Coxeter’s theoretical model.8 Thehistogram of the number M of edges per face verifies the predominance of five-fold,and shows traces of eight-, nine-, and even ten-fold faces. is 5-160, higher thanCoxeter’s 5.1 15.8, 9These higher values immediately suggest two points. (a) In the plasticene model,very small faces are likely to be difficult to observe. Taking a sphere diameter about*”, and faces less than 1 mm2 in area to be unobserved, the Voronoi analysis of theScott model shows that about 130 faces come below this limit. for the plasticenemodel is thus increased to about 13-9. This factor does nothing to remove thediscrepancy with Coxeter’s theoretical model.(b) If the packing density were lowerthan the maximum limit, J7 would be greater than the Rmin expected at that limit.Moreover, as fig. 5 suggests, a decrease in the spread of the density histogram keepingthe same average would decrease N, and a change in polyhedron distribution of a more66 RANDOM CLOSE-PACKED HARD-SPHERE MODELgeneral nature could have a similar effect. However, with regard to changes indistribution, we must bear in mind the possible restrictions imposed by the geometrvof the packing, as discussed above.Analysis of different density packings will provide some idea of how NI (and hence a) vary with density : although it may be that several different values of NI may occur-I214012cI008 06 04020,I13 - 14number of facesFIG.6.-Histogram of number of faces per polyhedron.2 8 0 02400.2000.VJ1600.44lw t3 '200.d800.400.0' 3 b 1number of edgesFIG. 7.-Histogram showing distribution of number of edges per faceJ . D. BERNAL AND J . L . FINNEY 67for the same density, it may be possible to deduce a relationship between overalldensity and N. As the values of N and density for the Scott model fit reasonablywell on the N against density curve of fig. 5, an investigation of local regions ofdifferent densities may be of interest in this connection.TOPOLOGICAL CHARACTERISTICSIf the examination of some characteristics of a set of polyhedra showed theoccurrence of certain " types ", in the simplest case an array could be described interms of the frequency of occurrence M i of type Ai.This system, however, gives usno information about how different polyhedra are spatially arranged : an extensiontherefore is to formulate a matrix showing the frequency of occurrence nil of type Acnext to type Aj :A further extension adds a third dimension to the matrix, showing the type of faceBk by which Ai and Aj are joined. The general case occurs when all polyhedra are ofdifferent types, SO that nijk = 0 or 1.One possible characteristic by which similarity between polyhedra can be assessedis the topology. In this connection, it might be worth elucidating all the possibletopologically different types of polyhedra with N faces, and see whether any of thesetypes are specificially excluded or frequently occurring in a Voronoi network. Thistype elucidation has been done for N = 4 to 12 inclusive, after which the task of furtherextension for large N becomes unmanageable.If interesting results come out ofcomparing our 24 twelve-faced polyhedra with the possible 77 types, it may be worthconsidering further work along these lines.PRELIMINARY TOPOLOGICAL RESULTSA predominance of certain types is evident. Table 1 shows those polyhedrontypes with five or more examples present, together with the full plasticene sphereanalysis for comparison. 255 of the 407 polyhedra (63 %) are thus accounted for :the total number of plasticene polyhedra was 65. The number of exact correspond-ences between the two models, especially in the six and eight pentagonal face groups,is worth noting.And there are other near correspondences between polyhedra differ-ing only slightly, e.g., 04460 and 13451, where one quadrangular face has becometriangular, and one hexagonal has become septagonal ; similarly with 13532 and0454 1.What this representation fails to do is to allow for the differences in relativearrangements of faces; for example, there are five distinct polyhedra of type 03630(the 6 referring to the pentagonal faces). However, we now have a programme tosort out these sub-types. So far we have found that even where two polyhedra hav68 RANDOM CLOSE-PACKED HARD-SPHERE MODELdeintical arrangements of the same types of faces, they are not quantitatively similar.A case in point is polyhedra 1 and 322 in table 2, which have only four faces (1,9, 10, 13) corresponding vaguely in size and distance from the centre.Can two suchTABLE 1A B(Scott model) (plasticcne)1141 50( 1)4B04420(1)04440(1)04450(7) 0#50( 1)04460( 1 1) 04460(1)04470(5)1 343 1 (1)13441(2)1345 l(l0) 1345 l(3)13461 (2)6A B03620( 1)03630(?) 03630(5)03640(30) 03640(5)03650(17) 03650(3)03660( 17) 03660(5)12641 (5)04642( 1)8A B028 1 O(2)02820(6) 02820(3)02830( 13) 028 30( 3)02840( 15) 02840(4)02 8 50( 1 0)02860( 1)11831(8)10A B01”10”20(16) 01“10”20(4)01”10”30(11)01 ”1 0”40(2)A B30533 l(1)5A B13532(6)04541 (2)0453 1 (1 )12530( 1)7A B03 72 1 (3)0373 1 (1)0374 l(8)12722(7)12732(11)9A B10930( 1)02951 (1)12A B00” 1 2”00( 1 )00” 12”20( 1)01 ”1 2”00(5)polyhedra be considered to be similar for the purposes of our required geometricaldescription, and if not, how can the differences be represented? Much remains to bestudied merely with respect to this one analysis, which may throw further light onproblems such as these.FUTURE WORKThe polyhedron analysis of this and other models will be pursued, as variouslymentioned above, using in addition auxiliary tools such as computer programmes toproduce stereographic projections and stereoscopic pictures of the polyhedra. WJ .D . BERNAL AND J . L . FINNEY 69intend to measure a large (-3000 centres) random close-packed model on our co-ordinate measuring machine, which will give a much greater accuracy of measurementthan anything at present available : the increased size will give greater flexibility andstatistical reliability, and more scope for investigating local density fluctuations.Analyses of models of different densities may give some indication of how to deal withTABLE 2faceno.1234567891011121314no.ofedges75455664645555face area1 3225.84.60.34-33-65.04.90.95.00-53.74.24-20.55.72-72.73.31 *75.85-62-34-80.22.75.14.22.8distance betweencentre and face1 3221.581 -632.131-661.68.1-601.592.051 -592.171 -601.601-642.181.611.921 *891 *782.051.591.621 -921 6 02.281 -781-591.591-81temperature variations (such as by polyhedron dilatations or a change in the distribu-tion of types). The use of the Voronoi polyhedron as a basis for thermodynamiccalculations is to be investigated, and we are also intending an approach to the thermo-dynamics via Monte Carlo computer techniques. Along these lines we see thedevelopment of further consequences of a random-packed liquid, and thus hope toprovide more points for experimental comparison.1 J. D. Bernal, Proc. Roy. SOC. A , 1954, 280, 299.2 J. D. Bernal, Nature, 1960, 185, 68.3 J. D. Bernal and J. Mason, Nature, 1960, 188,910.4 G. D. Scott, Nature, 1960,188,908.5 G. D. Scott, Nature, 1962, 194,956.6 J. D. Bernal, J. Mason and K. R. Knight, Nature, 1962, 194,958.7T. Kiang, 2. Astrophysik, 1966, 64, 433.8 H. S. M. Coxeter, Introduction to Geometry, (John Wiley, 1961), p. 411.9 J. D. Bernal, in Liquids ; Structure, Properties, Solid Interactions, ed. T. J. Hughel. (ElsevierPub. Co. Amsterdam, 1965), p. 25

ISSN:0366-9033    

DOI:10.1039/DF9674300062

出版商:RSC    

年代:1967

数据来源: RSC

摘要:

Hidden Variables in the Critical Region*B Y MARSHALL FiXMAKChemistry Dept., Yale University, New Haven, Connecticut, U.S.A.Receiued 19th December, 1966A gas in the one-phase region near the critical point exhibits unusually large fluctuations ofdensity in volume elements much larger in linear dimension than are molecular diameters. Thesefluctuations in density give rise to forces and energy fluxes which may be treated by essentiallymacroscopic equations of thermodynamics and hydrodynamics, because the dominant wavelengthsare so large in a Fourier decomposition of the density fluctuation. The anomalous thermal con-ductivity and shear viscosity of a gas near the critical point have been calculated, and comparedwith experiments on CO,. The comparison of thermal conductivities is satisfactory, but thecomparison of shear viscosities is inconclusive, because of the smallness of the predicted effect andthe difficulty of its measurement.1.INTRODUCTIONIn a variety of relaxation problems, an unexpected dissipation of energy is explainedthrough the existence of an uncontrolled or " hidden " variable, such as an intra-molecular mode of energy, whose probability distribution is unable to remain inequilibrium with the varying hydrodynamic variables. The anomalous ultrasonicabsorption of gases near the critical point has been explained on this basis,l* withthe hidden variables taken to be the spectrum of long-wave-length density fluctuations.The density fluctuations appear in an anomalous contribution to the entropy, fromwhich a complex, dynamic heat capacity was derived and used to calculate the soundabsorption.In the previous work the anomalous entropy varied in time because ofa uniformly varying temperature, but here we consider a non-uniform variation oftemperature and from the variation in entropy derive a formula for the anomalousheat flux. The anomalous shear viscosity of the gas is also calculated, by techniquesmathematically similar to those used for the shear viscosity of the binary mixturein the critical region (although the physical bases are different). The thermalconductivity calculation may significantly be compared with experiment, but theanomalous gas viscosity is calculated to be so small, and the discrepancies betweendifferent experiments are so large, that no useful comparison can be made.Zwanzig and Mountain have also calculated the anomalous thermal conductivityand shear viscosity of a gas.Apparently, their results for the thermal conductivityare qualitatively similar to ours, but they predict no anomaly in the viscosity. Theirmethod of calculation differs from ours in two respects. First, they calculated thetransport coefficients from correlation function expressions (as did K a ~ a s a k i , ~ inhis calculations on binary mixtures). Secondly, they retained explicitly in theirequations of motion an intermolecular potential. The use of correlation functionsmakes unnecessary any search for a coupling between the motion of the fluid and thedriving force for shear flow or thermal conduction. We have thought it worthwhileto find this coupling, as something interesting in itself and also to provide an insight* supported in part by the Public Health Service, GM 13556, and the Alfred P.Sloan Foundation.7M . FIXMAN 71into non-linear transport processes. Regarding the intermolecular potential,Zwanzig and Mountain used a long-range potential to justify some of the manipula-tions. We use an alternative justification which consist in the restriction of densityvariations to those describable by long-wave-length Fourier components. Thepotential itself is presumed to be short range and to manifest itself only implicitlyin the various thermodynamic derivatives that are present. The quasi-thermodynamicforms of our equation facilitates a numerical comparison with experimental results.Only one parameter occurs which is microscopic in nature,6 the short-range correlationlength 1.In principle, 1 may be determined from light scattering, but in practice Iis here determined for COz from a previous comparison of heat capacity predictionswith experiment.The following remarks are essentially a summary of methods and results, thedetails of which must be published elsewhere for lack of space here.2. EQUATIONS OF MOTIONWe assume that at each point a specific energy and specific entropy are known,and that the specific entropy is a function of the local density and temperature, whilethe specific energy contains one ordinary part, a function of the specific entropy andmass density, and another part proportional to the square of the density gradient.ThusI-Je = eL+eI, (2.3)e, = a I Vp 1 ’.eL is the ordinary specific energy, while e, becomes significant in the critical region,not because a is supposed to become large, but because the second derivitive of eLat constant s becomes small near the critical point, and vanishes as the reciprocalcompressibility.For a given state of density fluctuation, the force acting on any unit element ofmass is presumed to be given bydu/dt = p - l V * a0+F, (2.5)whereq) = 2q*(Vu)”Yrn C2.6)and ro is a hypothetical shear viscosity which the fluid would have if density fluctua-tions were not allowed.The force which arises because of the density fluctuationsis designated F in eqn.(2.5).The force F is calculated on the basis of a mechanical variational principle. Theposition ro of any element of mass is taken to be a function of time and the positionw of the element in an initial reference configuration. The positions w are takento be stochastic variables, and an implicit average over their possible values is requiredwhen macroscopic quantities are calculated :7ro = r&(w,t72 HIDDEN VARIABLES I N THE CRITICAL REGIONThe mechanical variational principle is that6E = - F Grodm,dm = p(r,t)dr, swhere Gr, is the variation in the path of a mass element, the variation being at aconstant value of the entropy of the mass element. A straightforward calculationyieldspF = V * G,,, (2.9)(2.10) Gp = -p1+2a[p 1 vp I "p2v2p]1-2ap(vp)(vp).The last term on the right-hand side of eqn.(2.10) is of utmost significance attwo points in the analysis. First, with the consistent neglect of third-order fluctua-tions, and the consistent interpretation of p as an average over a small volume element(with dimensions much larger than the range of the intermolecular potential, butmuch smaller than the range of correlation in the critical region), an ensembleaverage of eqn. (2.5) gives- qo(e,eZ + eZeJ + 2ap lim d 2G(1)(r)/dxdz. (2.11)r-tOIn this equation, i: is at the rate of shear of the fluid in laminar flow, G(I)(r) is theperturbation in the radial correlation function evaluated to first order in i. Thelimiting process in eqn. (2.11) is not taken towards a mathematical zero, but to avalue of r much smaller than the range of the correlation.The second point wherethe last term in eqn. (2.10) plays a major role is in a thermal conduction experiment,where this term provides the coupling between the uniform temperature gradientand the density fluctuations. We derive below an expression for the anomalousthermal conductivity.3. HEAT CONDUCTIONFor the anomalous specific entropy, Botch and Fixman derived an expressionequivalent to(6s) = (*)s:uPi 4<(6P) 2> (3.1)where the subscripts on the right-hand side designate the second derivative of thespecific entropy with respect to the specific volume. Interpretation of p as anaverage over a small volume element led to a transformation of the right-hand sideto an integral over the long wave-length components of the radial correlation function.Then, a differential equation for the radial correlation function, the one introducedbelow, allowed a dynamic heat capacity to be calculated.However, there is anotherway in which the anomalous specific entropy can vary with time, and that is by heatconduction. In the absence of heat sources, an anomalous heat flux can be relatedto (6s) bypoToa(ss>/at = -v . q', (3.2)and eqn. (3.1) and (3.2) together giveq' = T ~ S , , ~ ; ~ Iim u(r),r-rO(3.3)where n(r) is the mean fluid velocity at a point r due to the combined effects of amolecule fixed in location at the origin, and the presence of a temperature gradientM. FIXMAN 73We have calculated u(r) for a fluid in the critical region, subject to a uniformtemperature gradient independent of time.This velocity must be determined fromthe modified Navier-Stokes equation, eqn. (2.5) together with eqn. (2.9) and (2.10)for the forces. The linearized equation of mass conservation shows that the diver-gence of u must vanish in the steady state; that is, u can be written as the curl of avector potential. Therefore any part of eqn. (2.5) which can be written as the gradientof a scalar potential will not contribute to the desired solution. The rejection of suchterms givesq0V2u = 2p&x(ap/aT),(VT)V2Go(r), (3 -4)G"(r) = a e-Kr/r, (3.5)where Go is the equilibrium radial correlation function and K is related to the thermo-dynamic parameters byK~ = -pv/(2ap:), (3.6)or, at the critical density, by Debye's relation,I C ~ = 6( T - Tc)/12Tc.Eqn.(3.4) still includes some scalar potential contributions to u, but these may berejected and a thermal conductivity determined from eqn. (3.3). The result for theanomalous conductivity A' isA' = -(3)Toaq, 1(ap/aT),(a2p/aTa In U ) K - l. (3-7)4. SHEAR VISCOSITYThe fluctuation formula for the radial correlation function,(MrIFjP(r2)) = P m l ~ r 2 ) (4.1)and eqn. (2.5) linearized in the density fluctuations, together with the ordinary thermalconduction equation gives a result, previously derived by Botch and Fixman,hdG(r)/dt = h(aG/at + [u(r) - u(O)] * VG) = V2[x2G - 2VG], (4.2)where ;lo is, like yo, a hypothetical quantity, the thermal conductivity in the absenceof long wavelength density fluctuations.Like yo, Ro must be determined by extra-polation from outside the critical region. The previous expression for h was too largeby a factor of two, due to the neglect of the fact that two time-dependent densityfluctuations enter the formula (4.1) for G. A first-order solution of eqn. (4.2) forthe perturbation induced by simple shear flow is straightforward, and gives onsubstitution into eqn. (2.1 1) a formula for the anomalous part of the shear viscosity :q' = aTop+/(40Ao~). (4.4)A virtually identical formula, with different interpretation of the coefficients, has beenderived for the shear viscosity of a binary liquid mixture by severalThe present method may also be used to reproduce the previous results.* *5.COMPARISON WITH EXPERIMENTFor COz in the critical region all of the thermodynamic and transport coefficientsrequired in eqn. (3.7) and (4.4) are available. The value of I was previously estimate74 HIDDEN VARIABLES I N THE CRITICAL REGIONto be 4.5& and this value has been used here. The numerical predictions are, atthe critical density :y’ = 38( T - Tc)-* micropoise,2.’ = 2-6 x 10-4(T- Tc)-* cal/deg. cin sec.(5.1)(5.2)The predicted y’ is some five times smaller than indicated by the data of Michelset a2.,l0 measured in a capillary viscometer, and some five times larger than indicatedby the data of Kestin et a2.,l1 measured in a rotating disc viscometer. The latterexperimental anomalies were increased by factors of about three in sample correctionsfor density gradients.Eqn. (5.2) has been compared with Sengers’ data,I2 the total conductivity beingtaken to be ;1 = Lo +A’, with A,, = cal/deg. cm sec. Over a 40°C intervalabove T,, the calculated and observed conductivities agree within 10 %, althoughthe observed results show a slightly sharper dependence on T than given by eqn. (5.2).1 W. Botch and M. Fixman, J. Chem. Physics, 1965,42, 196.2 W. Botch and M. Fixman, J. Chem. Physics, 1965, 42, 199.3 M. Fixman, Adv. Chem. Physics, 1963, 6, 175.4 R. Zwanzig and R. D. Mountain, private communication.5 K. Kawasaki, Physic. Rev,, 1966, 150, 291.6 In principle even this parameter may be determined from macroscopic observations on inhomo-7 The explicit computation of averages proceeds in terms of correlation functions which are the8 M. Fixman, J. Chem. Physics, 1962, 36, 310.9 J. M. Deutsch and R. Zwanzig, preprint.10 A. Michels, A. Botzen and W. Schuurman, Physica, 1957, 23, 95.11 J. Kestin, J. H. Whitelaw and T. F. Zien, Physica, 1964, 30, 161.12 J. V. Sengers, Thesis (University of Amsterdam, 1962).geneous systems.solutions of differential equations

ISSN:0366-9033    

DOI:10.1039/DF9674300070

出版商:RSC    

年代:1967

数据来源: RSC