Acta Chim. Slov. 1999, 46(3), pp. 375-388 THE CALCULATION OF THE THERMODYNAMICAL PROPERTIES IN THE LIQUID-GAS REGION Jurij Avsec, Milan Marcie Faculty of Mechanical Engineering, University of Maribor, Maribor, Slovenia (Received 25.8.1998) ABSTRACT The paper features the mathematical model of computing phase diagrams and thermodynamic functions of the state in the liquid, gas and two-phase domain with the help of statistical thermodynamics. The paper features all important contributions (translation, rotation, internal rotation, vibration, intermolecular potential energy and influence of electron and nuclei excitation). To calculate the thermodynamic properties of real gases we developed the cluster theory, which yields better results than the classical virial equation. For the realm of real liquids the Johnson-Zollweg-Gubbins model based on the modified Benedict-Webb-Rubin equation was applied. The Lennard - Jones intermolecular potential was used. The analytical results are compared with the thermodynamical data and models obtained by classical thermodynamics and show relatively good agreement. INTRODUCTION In engineering practice processes in the liquid-gas region are of vital importance. In order to design devices for this field of activity, it is necessary to be familiar with the thermodynamic properties of state in a one-phase and a two-phase environment. In most cases the thermodynamic tables, diagrams and different empirical functions obtained from measurements are used. The calculation of thermodynamic properties of state is possible by classical or statistical thermodynamics. Classical thermodynamics has no insight into the micro structure, but it allows the calculation of the thermodynamic function of state with 376 the assistance of macroscopic observation of phenomena. Unlike the classical thermodynamics, however, the statistical thermodynamics does enable the computation of the thermodynamic functions of the state by studying the molecular structure of matter. This paper presents a mathematical model for computing the speed of sound and other thermophysical properties using statistical thermodynamics. For real gases we developed a cluster theory, based on the principle of average clusters, which yields better results than the virial equation. For real liquids the Johnson-Zollweg-Gubbins model based on the modified Benedict-Webb-Rubin (BWR) equation of state and a great number of Monte Carlo and molecular dynamics simulations were applied. For the calculations in the two-phase region we applied the method of equilibrium conditions between two phases. The mathematical model enables the calculation in both sub- and supercritical region. To establish the critical point, the stabilization conditions known from the theory of fluctuation were applied. Using the mathematical model described above we were able to compute the thermophysical properties in the one- and two-phase region and draw phase diagrams for some technically significant substances. The results of the analysis are compared with the experimental data and show a relatively good agreement, especially for real gases. Somewhat larger deviations can however be found in the real liquid region due to the large influence of the attraction and repulsion forces, since the Lennard- Jones potential is an approximation of the actual real intermolecular potential. To calculate the thermodynamic functions of state we applied the canonical partition [3]. Utilising the semi-classical formulation for the purpose of the canonical ensemble for N indistinguishable molecules the parition function Z can be expressed as follow [2,3]: Z = N!h Nf-f"J'eXxp ~kT dr1dr2..drNdp1dp2..dpN (1) where f stands for the number of degrees of freedom of individual molecule, H designates the Hamiltonian molecule system, vectors r1,r2..rN .. describe the positions of N molecules and p1,p2...pN momenta, k is Boltzmann s constant and h is Planck`s constant. 377 The canonical ensemble of partition function for the system of N molecules can be expressed [3,5]: Z = Z0Ztrans ZvibZrotZir ZelZnucZconf . Thus the partition function Z is a product of terms of the ground state (0), the translation (trans), the vibration (vib), the rotation (rot), the internal rotation (ir), the influence of electrons excitation (el), the influence of nuclei excitation (nuc) and the influence of the intermolecular potential energy (conf). The computation of the individual terms of the partition function and their derivatives with the exception of configurational integral is dealt with in the works [3,4,5]. The velocity of sound The term velocity of sound refers to the velocity of the mechanical longitudinal pressure waves propagation through a medium. It is very important parameter in the study of compressible fluids flows [6,28] and in some applications of measurement (acoustic resonance level gauge) [24]. The propagation of sonic waves for real fluids is almost in all cases nearly isentropic [7,8,9]. Therefore, we can calculate the isentropic speed of sound for the real fluid c0: 21 dp ^ c0 = J—v I —^- I ) cp(dT\ p T I 3v I ¦v----------—-------j-—^—,( dv | cp | 3T | 3tJp T dp 1 /v where T is temperature, p is pressure, s is intropy, v is the specific volume, cp is the specific heat capacity at constant pressure. 378 CONFIGURATIONAL INTEGRAL The metho ds for solving the configurat io nal integral in the liquid - gas regio n are generally divided most crudely into solving the configurational integral for real gases and for real liquids. In a real gas the intermolecular forces are relatively weak. To solve the configurational integral we devised a method of clusters [1]. The real gas molecules move either individually or in small instantaneous and random clusters. The size and shape of the cluster alter due to the existence of attractive and repulsion forces between molecules. The number of molecules in clusters depends on the intermolecular distances. Molecules which are sufficiently close to one another are considered to be in the same cluster. More distant molecules, though close enough to each other, are said to be in another cluster. The method is based on the principle of average clusters. Fig. 1: Schematic outline of clusters. By dealing with a sufficiently large number N of molecules in the system an average cluster of N1 molecules can be determined. Figure 1 features the outline of the idea of the metho d o f clusters illust rating the activit y of int ermolecular fo rces in clusters as well as the a ctivi ty of in termolecular forces between clusters. The analysis of the motion of molecules by means of the Monte Carlo method [10] indicates that the number of molecules in the average cluster is relatively low. The principal idea of t he met hod of clust ers st ipulates t hat the intermolecular po ten tial energy of the system can be split up to intermolecular potential energy in clusters (Epot1) and 379 intercluster potential energy due to interactions between clusters (Epot2). Therefore the potential intermolecular energy can be written as the sum of both parts: E -F +E -^pot ~~ Epot1 T -Lpot2 (4) On the basis of the previously mentioned assumtions the configurational integral can be put as follows: Z =Z Z conf conf 1 conf 2 . (5) In the study of systems with the intermolecular potential uij acting between molecules i,j, it is more suitable to use the Mayer function fij: fij=exp(-uij/kT)-1. (6) For the intermolecular potential the Lennard-Jones potential [11,12] was applied: uij=4e r-1 v J J r-6 V J J (7) In equation (7) g and e are Lennard-Jones parameters and rij is the intermolecular distance. When computing the configurational integral Zconf1 the effects of mutual interactions of up to three molecules in the cluster were taken into account. The configurational integral Zconf1 that takes interactions in clusters into account can be written [5,11,12]: ^conf1 Vi J...Jdri..drNi + J...J Lfydri..drNi + J...jL(fijfkl)ir1..drNi +. (8) 1 N Equation (8) can also be expressed as follows [25]: ln Z conf1 N n7 •ln 1+ (N1-1)NI1 + ( N2 2V 1 6V2 (9) 380 Integrals I1 and I2 are defined: N, N I1 = vllf« dr1 dr2, I2 = ^TT" J/J C3fy fkl +fij fkl fmnldr1 dr2 dr3 (10) V Integral I1 for the case of the Lennard-Jones potential is solved by literature [12,13,14]. The integral I2 was transformed by A.Münster [11] as follows: I2=3I12+I21, (11) I = N1 21 V JJJfijfkrfmndr1dr2dr3 Integral I21 is resolved in the works of LE. Reichl [13] and Hirschfelder, Curtiss and Bird [12] for the case of the Lennard-Jones potential. The effect of the intercluster interaction can be presented by the configurational integral Zconf2 which includes the mutual interactions of two clusters in the system. U, Jconf2 -1NJ.... exp V N1 __ij_ kT dr1dr2...drN/N1 , (12) where Uij stands for the potential energy between i- and j- clusters. Applying Corner s method [26] in the case of diluted systems the inercluster potential energy can be expressed: U^N12uij, (13) where uij stands for the potential energy between molecule in the cluster i and molecule in the cluster j. Utilizing the above and taking only the mutual interaction of two clusters in the system we can write: Fij = exp ( U-- ^ * 1 rr N __ij_ kT (14) / where N2 and Fij are the number of clusters in the system, and Mayer function. Integral I1 can be solved in the same way as I1. lnZconf2= + N, 2V L1 ¦> (15) To determine the number of molecules in an average cluster of real gases we took advantage of the experimental results of thermodynamic functions of the state obtained 381 by the experimental results of J.B. Maxwell [15] and W.C. Edmister [16], as well as J.A. Barker [10] where the results of simulation by means of the Monte Carlo method are reported. After a thorough analysis the number of molecules in individual average clusters was established: 1 < Nj < 6 . For a real liquid we used the Johnson-Zollweg-Gubbins (JZG) model [17] based on molecular dynamics and Monte Carlo simulations and the modified Benedict-Webb-Rubin equation of state. The JZG model is based on Lennard-Jones intermolecular potential and contains 32 linear parameters and one nonlinear parameter. On this basis we can express reduced configurational free energy Aconf : Aconf *=I^ + IbiGi, (16) i=1 i i=1 where the coefficients ai, bi and G are presented in Table 1. The coefficients ai and bi are functions of reduced temperature T only, the coefficients Gi are function of the reduced density p and nonlinear adjustable parameter y. P =-V-, T = — , Aconf =-NT, F = exp(-yp2j ,7=3. (17) Equation (16) accurately correlates the thermophysical properties from the triple point to about 4 to 5 times the critical temperature. In equation (16) xj´s are the adjustable parameters in the equation of state. In equation (17) Aconf is conigurational free energy and V is volume of molecular system. To determine the equilibrium states between the liquid and the gaseous phases the conditions for equilibrium were applied: = T , p =p , n =\i, (18) where ´ in equation (18) means the liquid phase, " means the gaseous phase and \i constitutes the chemical potential. Due to the mathematical complexity of the equations in the model, the states on the coexistence curve are obtained numerically, with the help of the iteration procedure. By applying these states, the thermodynamic properties in the two-phase environment can be calculated. In the two-phase region the computation of thermodynamic functions of state is based on the mixing rule [7]. 382 Table 1: Coefficients for the JZG model i 1 2 3 4 5 6 7 8 ai bi Gi x1T* +x2vT* + x3 +x4 /T* +x5 /T*2 x20 /T*2 +x21 /T*3 (1- F)/(2g) x6T +x7 + x8 /T* +x9 /T*2 x22 /T*2 +x23 /T*4 -(Fr*2 - 2G1 )/(2g) x10T*+x11+x12/T* x24 /T*2 +x25 /T*3 -(Fr*4 - 4G2 )/(2g) x 13 x26 /T*2 +x27 /T*4 -(Fr*6 - 6G3 )/(2g) x14T*+x15/T*2 x28 /T*2 +x29 /T*3 -(Fr*8 - 8G4 )/(2g) x16/T* x30 /T*2 +x31 /T*3 +x32 /T*4 - (Fr*10 -10G5 )/(2g) x17 /T* +x18 /T*2 x19 /T*2 RESULTS AND COMPARISON WITH EXPERIMENTAL DATA The constants necessary for the computation, such as the characteristic rotation-, vibration-, electronic- etc. temperatures, are obtained from the experimental data [18,19,20]. The inertia moments are obtained analytically by applying the knowledge of the atomic structure of the molecule. Constants for the Lennard-Jones potential are obtained from the literature of [3, 4, 12]. We carried out calculations for methanol (CH4O), difluordichlormethane R-12 (CF2Cl2), ammonia (NH3), benzene (C6H6) and carbon dioxide (CO2) . The comparison of our calculations with thermodynamical data and experimental results (exp.) by Maxwell [15], Edmister [16], Eckert & Drake [21], Petrak [22] and Younglove, Frederick, and McCarthy [9], Borgnakke and Sonntag [23] are presented in Tables 3,4,5,6. Some important constants (vibration (qvib) and electron (qel) characteristic temperatures, internal rotation moment of inertia (Iir), symmetry number for internal rotation (sint)) need for the calculation of thermodynamic properties are presented in Table 2. 383 Table 2: Some important constants [18,19,20] Comp 2\ Iir (gern ) sir ---3 - 0el(K) qvib (K) CO2 - - 1890, 3360, 954, 954 NH3 - 66379 4780, 1360, 4880,4880, 2330, 2330 R-12 - - 374,654,623,955,460,1262,1322,1566, 1650 CH4O 1.027-10-40 - 5295,4279,4093,2127,2052,1936,1485, 1546,4279,2092,1769 C6H6 - - 581,581,869,869,969,1119,1215,1415,1419,144 1,1572,1572,1637,1658,1673,1873,1701,2116, 2352,2651,2552,2173,4500,4700,4900,5100, 5300, 5500, 5700 Fig. 2: Entropy and enthalpy of vaporisation for R-12 and C6H6 384 Table 3:The comparison between the analytical calculation and the experimental data [21,22] for difluordichlormethane in the region of superheated vapour T (K), H (kJ/kmol), S (kJ/kmolK), V (m3/kmol) T 10 bar 20 bar V H S V H S 373 JZG 2.81 29120 97.1 1.25 28440 91.0 ST 2.80 29060 98.5 1.25 28270 91.1 CT 2.80 29990 97.9 1.22 28280 89.2 Exp. 2.79 29936 97.9 1.21 28646 89.6 423 JZG 3.30 36900 106.1 1.55 35380 105.4 ST 3.30 33190 109.0 1.55 32560 102.2 CT 3.29 34880 109.3 1.56 33950 102.4 Exp. 3.29 34637 109.1 1.53 33358 101.8 473 JZG 1.82 40840 114.1 ST 1.80 36950 111.1 CT 1.82 38760 109.3 Exp. 1.79 38274 104.5 Table 4:The speed of sound c0 for carbon dioxide (CO2) in the region of saturated vapour T (K), V (m3/kmol), c0 (m/s) T V M O D E L S VDW RK PR JZG CT Exp. 223 2.432 222.1 221.6 219.7 226.7 224.1 223.6 233 1.667 223.9 222.7 220.9 227.5 224.2 224.2 243 1.144 221.8 222.3 219.8 226.9 222.6 222.9 253 0.843 220.9 222.2 219.6 226.2 221.1 220.8 263 0.613 216.1 220.1 218.4 223.5 217.9 217.5 273 0.448 209.0 214.2 217.0 219.3 213.9 212.7 283 0.325 202.4 205.5 214.7 212.9 209.9 206.3 293 0.228 188.0 181.9 204.8 203.7 208.6 197.9 301 0.155 180.2 92.9 139.1 194.1 219.3 189.8 385 Table 5:Pressure of saturation for ammonia (NH3) T (K), V (m3/kmol), p(bar) T V M O D E L S VDW PR JZG CT ST Exp. 223 44.76 0.43 0.40 0.41 0.41 0.41 0.40 243 16.37 1.29 1.21 1.22 1.20 1.22 1.19 263 7.072 3.53 3.42 3.03 2.98 3.03 2.91 283 2.941 7.61 7.41 7.65 6.20 7.60 6.16 303 1.853 12.63 12.18 12.77 11.73 12.76 11.6 323 1.067 21.91 21.74 22.90 20.61 22.85 20.33 343 0.643 35.43 34.30 38.78 34.36 38.66 33.07 363 0.395 57.18 52.54 63.71 54.40 63.3 51.02 383 0.238 83.70 76.90 106.5 83.20 104.0 75.60 We compared the models obtained by statistical thermodynamics (CT- cluster theory, ST- statistical theory with help of classical virial expansion [14], [12], JZG-Johnson-Zollweg-Gubbins model [17]), and the models obtained by classical thermodynamics (VDW- Van der Waals equation of state, RK- Redlich-Kwong equation of state, PR- Peng-Robinson equation of state) and experimental results (exp.)). In the procedure of computing the configurational integral with the help of statistical thermodynamics on the base of classical virial expansion (ST) the effects of mutual interactions of up to three molecules in the cluster were taken into account. Figure 2 illustrates the thermodynamic functions of the state computated by means of the mathematical model. Figure 2 shows the enthalpy of vaporisation and the entropy of vaporisation for benzene (C6H6) and difluordichlormethane (R-12). The comparison (Figure 2) of analytical calculations and thermodynamical data [7,21,22] for R-12 shows good agreement. In the real gas region we used the cluster theory (CT), and in the real liquid region we used the Johnson-Zollweg-Gubbins (JZG) model. The computed thermophysical functions for benzene show a slightly less good agreement of results. The reason lies in the complex structure of benzene molecules. Table 3 shows the deviation of the analytical computation of difluorodichlormethane for the real gas region (CT, ST, JZG) from the experimental values (exp.). The results obtained by the cluster theory show very good agreement. The cluster theory (CT) 386 yields better results, especially at high pressures. Tables 4 and 5 show the deviation of the analytical results obtained with statistical thermodynamics (CT, ST, JZG) for carbon dioxide and ammonia from the models obtained by classical thermodynamics (VDW, RK, PR) and experimental data (exp). The computed vapour pressure and velocity of sound conform well with the measured vapour pressure and velocity of sound. Somewhat larger deviations can be found in the region of critical conditions due to the large influence of fluctuation theory [27] and singular behaviour of some thermodynamic properties in the near-critical condition. The cluster theory (CT) yields better results, especially at high pressures, than the classical statistical virial form of the equation of state (ST). Somewhat larger deviations are observed in the real liquid region (Table 6) due to the large influence of attraction forces, since the Lennard-Jones potential is only an approximation of the actual real intermolecular potential. Table 6:The comparison between analytical calculations and experimental data for boiling liquid for methanol (CH3OH). T K P bar V m /kmol H kJ/kmol S kJ/kmolK 310 JZG 0.32 0.044 21894 110.6 Exp. 0.30 0.042 20214 128.4 350 JZG 1.81 0.047 25532 130.3 Exp. 1.60 0.044 23581 141.3 410 JZG 10.3 0.054 30347 166.7 Data 10.0 0.049 29366 157.4 450 JZG Exp. 26.5 25.3 0.063 0.055 33628 34214 178.3 167.4 387 CONCLUSION AND SUMMARY The paper presents the mathematical model for the computation of thermodynamical functions of the state in the liquid, gaseous and two-phase region. For the region of real gases we developed the method of clusters based on the average cluster with N1 molecules in the system with N molecules. The advantage of this method is in the computation of thermodynamical functions at high pressures. For the real liquid we used the Johnson-Zollweg-Gubbins [17] model based on molecular dynamics, Monte Carlo simulations and the modified Benedict-Webb-Rubin equation of state. The boiling curve and the saturation curve were determinated by means of equilibrium conditions. The analytical results obtained by statistical thermodynamics were compared with the thermodynamical data and the analytical calculations obtained by classical thermodynamics. The results show a relatively good agreement. The computed thermodynamical properties in the region of real gas conform well in comparison with the thermodynamical data. In the real gas region the comparison between the thermodynamical data and the analytical calculation shows a slightly less good agreement. The reason lies in the effect of molecular polarity. For our further research we intend to extend the presented mathematical model to the calculation of mixtures and to include the theory of fluctuation which is of great importance in the vicinity of the critical point. REFERENCES [I] M. Rigby, , E.B. Smith, W.A. Wakeham, G.C.Maithland, G.C., The Forces between Molecules, Clarendon Press,1986, Oxford [2] J.A. Barker, D. Henderson , Reviews of Modern Physics,48, 1976, 587-671 [3] K. Lucas, Applied Statistical Thermodynamics, Springer - Verlag, 1992, N. York [4] C.G. Gray, .K.E. Gubbins, Theory of Molecular Fluids, Clarendon Press, 1984 Oxford [5] B.J. McClelland, Statistical Thermodynamics, Chapman & Hall, 1980, London [6] B.E.Gammon, D.R. Douslin , J. Chem. Phys., 64, 1976,. 203-218 [7] F. Bosnjakovic, Nauka o toplini, Tehnicka Knjiga, 1976, Zagreb [8] Y.A. Cengel, M.A. Boles, Thermodynamics, McGraw Hill, 1994, N. York [9] B.A. Younglove, N.V. Frederick, R.D. McCarthy, Speed of sound Data and Related Models for Mixtures of Natural Gas Constituents, NIST, 1993, Washington [10] J.A. Barker, Lattice Theories of the Liquid State, Pergamon Press, 1963, Oxford [II] A. Münster, Statistical Thermodynamics, Springer - Verlag, 1974, New York [12] O.O. Hirschfelder, C.F. Curtiss, R.B. Bird,. Molecular Theory of Gases and Liquids, John Willey & Sons, 1954, London. 388 [13] L.E. Reichl, A Modern Course in Statistical Physics, University in Texas, 1990 [14] T. Kihara, Intermolecular Forces, University of Tokyo, John Wiley & Sons, 1976, Chichester [15] J.B. Maxwell, Data Book on Hydrocarbons, Van Nostrand Company, 1955, New York [16] W.C. Edmister, B.I. Lee, Hydrocarbon Thermodynamics"I,II, Gulf Publishing Company, 1964, London [17] J.K. Johnson, J.A. Zollweg, K.E. Gubbins, Mol. Phys., 1993,78, 591-618 [18] G. Herzberg, Electronic Spectra of Polyatomic Molecules, Van Nostrand Reinhold Company, 1966, London. [19] G. Herzberg, Infrared and Raman Spectra of Polyatomic molecules, Van Nostrand Reinhold Company, 1984, New York [20] L.J. Bellamy, The Infrared Spectra of Complex Molecules, Chapman & Hall, 1980, London [21] R. L. Eckert, R.M. Drake, Heat and Mass Transfer, McGraw-Hill, 1959, New York [22] J. Petrak, L. Ludek, Thermokinetic properties of refrigerants, Faculty of Mechanical Engineering, 1993, Praha [23] C. Borgnakke, R.E. Sonntag, Thermodynamic and transport properties, John Wiley & Sons, 1997, New York [24] D. Đonlagić, , M. Završnik , D. Đonlagić, Sens. Act.,1996, A57, 209-215. [25] N.A. Smirnova, Methods of the Statistical thermodynamics in the Physical Chemistry, University of Moscow, 1982, Moscow [26] J. Corner, J., Proc. Roy. Soc, 1948, A, 192 [27[ A. Pelt, J.V Sengers, J. Supercr. Fluids, 1995, 8, 81-99 [28] D.C. Giancoli, Physics, Prentice-Hall International Editions, 1991, London, Sydney, Toronto POVZETEK Članek obravnava matematični model izračuna termodinamičnih veličin stanja v kapljevitem, plinastem in dvofaznem območju s pomočjo statistične termodinamike. Članek obravnava vse pomembne vplive (translacijo, rotacijo, notranjo rotacijo, vibracijo, medmolekularno potencialno energijo kot tudi vpliv vzbujanja elektronov in jeder). Za izračun termodinamičnih veličin stanja realnih plinov smo razvili model gruč, s pomočjo katerega smo dobili boljše rezultate kot s klasično virialno teorijo realnih plinov. Za področje realne kapljevine smo uporabili Johnson-Zollweg-Gubbinsov model, ki je zasnovan na osnovi modificirane Benedict-Webb-Rubinove enačbe stanja. Za medmolekularni potencial smo uporabili Lennard-Jonesov potencial. Analitični izračuni so primerjani s termodinamičnimi podatki in modeli dobljeni s pomočjo klasične termodinamike in kažejo na relativno dobro ujemanje.