Molecular Simulation for Thermodynamic Properties and Process

Article pubs.acs.org/jced

Molecular Simulation for Thermodynamic Properties and Process Modeling of Refrigerants William R. Smith,*,†,‡ Susana Figueroa-Gerstenmaier,¶ and Magda Skvorova§ †

Faculty of Science, University of Ontario Institute of Technology, Oshawa, Ontario L1H 7K4, Canada Department of Mathematics and Statistics, University of Guelph, Guelph, Ontario N1G 2W1, Canada ¶ Division of Science and Engineering, University of Guanajuato-Leon Campus, Leon, Mexico § Faculty of Science, J. E. Purkinje University, 400 96 Ú stí nad Labem, Czech Republic ‡

ABSTRACT: The calculation of mixture dew and bubble points and the simulation of fluid processes occurring under specified enthalpy or entropy changes are important in the refrigeration and other industries. These calculations are typically carried out using an empirically based multiparameter equation of state approach, such as that incorporated in the industry-standard REFPROP software package. We consider the capabilities of more fundamentally based molecular simulation methodology for performing these calculations. This approach requires a much smaller parameter set, can reliably extrapolate beyond the thermodynamic conditions used to determine the parameters, and can be used for the prediction of multiple properties. We briefly review relevant algorithms and focus on a set of Monte Carlo molecular simulation methods for accomplishing the tasks. We demonstrate their applications to the calculation of isoenthalps for the two pure fluid alternative refrigerants R134a (CH2FCF3; 1,1,1,2-tetrafluoroethane) and R143a (CH3CF3; 1,1,1-trifluoroethane), and to the simulation of all stages of a vapor-compression refrigeration cycle involving a binary mixture of the refrigerant R32 (CH2F2: difluoromethane) and R134a of 30 mass % R32. The molecular simulation algorithms produce results for these problems in good agreement with those calculated by REFPROP.

1. INTRODUCTION Advances in computer technology and molecular simulation methodology have enabled the increasing use of such tools as valuable adjuncts to traditional macroscopic approaches for engineering process design in many areas; the approach is described in several recent reviews.1−5 The methodology is based on the use of both quantum mechanics and statistical mechanics, the latter of which is of primary interest in this paper. It involves the implementation of algorithms to calculate properties and simulate processes of interest using Molecular Dynamics (MD) and/or Monte Carlo (MC) techniques, employing a set of mathematical models (force fields) for the interactions of the constituent molecules. Given the use of appropriate algorithms, the accuracy of the calculated results in comparison with experimental data is a reflection of the quality of the underlying force fields. Thermodynamic property prediction and process simulation for fluids is traditionally carried out using empirical equations fitted to experimental data. Such equations for many refrigerants have been incorporated in the software program REFPROP,6 which for each fluid employs a Helmholtz free energy equation of state (EOS), typically involving more than 40 parameters fitted to many hundreds of experimental data points. The extension to mixtures requires the use of empirical rules to calculate the cross-species parameters. In contrast, the development of force fields currently requires the order of five parameters to be fitted to a relatively small experimental data set for each fluid, and the predicted results for many fluid © 2014 American Chemical Society

properties have been shown to be of quantitative accuracy. (Although force fields that incorporate models of intramolecular forces contain additional parameters, these can be determined by quantum mechanical approaches not requiring additional experimental data.) Force fields have become available in recent years for many fluids, particularly those relevant to the petrochemical industry.1−5 Although ideal gas thermodynamic property prediction, required to calculate ideal gas contributions to thermal properties, is relatively well in hand by means of quantum mechanical approaches, progress in first-principles force field development for newly proposed fluids such as those that arise in refrigeration is still in a state of development. Although the vast bulk of force fields available in the literature have been determined by fitting their parameters to experimental results such as vapor−liquid equilibrium (VLE) data, with mixture force fields obtained by various methods based on those of the pure fluids, it has recently been shown that it is possible to develop transferable force fields (the force field parameters for a given interaction site are transferable between different molecules) for families of complex fluids such Special Issue: Modeling and Simulation of Real Systems Received: March 17, 2014 Accepted: July 10, 2014 Published: July 29, 2014 3258

dx.doi.org/10.1021/je500260d | J. Chem. Eng. Data 2014, 59, 3258−3271

Journal of Chemical & Engineering Data

Article

as fluoropropene-based refrigerants using a minimum of experimental data.7,8 In addition to force fields, the application of molecular simulation methodology requires algorithms to implement property predictions and process simulations. The majority of such applications have been in the areas of the prediction of volumetric and certain vapor−liquid equilibrium properties, for which relatively mature simulation methodologies are available. 9,10 Although Joule-Thomson inversion curves for refrigerants, which also depend solely on volumetric properties, have also been of interest,11−17 studies relevant to realistic industrial fluids for the prediction of caloric fluid properties, such as heat capacities11,18−21 and Joule-Thomson coefficients11,18−22 have been more recent. The calculation of isoenthalps and the simulation of processes at specified values of changes in the total system enthalpy, H,11,22−24 internal energy, U,22−24 or entropy, S,25−27 have also been the object of recent studies, as has the calculation of mixture bubble points28−44 and dew points.32−35,41−44 Calculations of these properties and for related processes are of intrinsic scientific interest and are also of practical importance in the refrigeration and natural gas processing industries. The purpose of this paper is to discuss molecular simulation methods for their calculation, and to demonstrate their use in examples from the refrigeration industry. The paper is organized as follows. In the next section, we describe the vapor-compression refrigeration cycle and thermodynamic conventions used. In the subsequent section, we discuss the computational methods employed. The next section presents and discusses new results for example applications to calculations of isoenthalps for the alternative refrigerants R134a (CH2FCF3; 1,1,1,2-tetrafluoroethane) and R143a (CH3CF3; 1,1,1-trifluoroethane), and to the simulation of a VCRC involving a refrigerant mixture of R32 (CH2F2; difluoromethane) and R134a. The paper concludes with a Conclusions section.

In Figure 2b, point 1′ represents the saturated vapor in the evaporator corresponding to the schematic of Figure 1 at its dew point, at (Plow, T1′). The vapor is then superheated by ΔTsh to point 1 at (Plow, T1 = T1′ + ΔTsh) before entering the compressor stage. The compressor of isentropic efficiency ω then compresses the vapor to point 2 at (Phigh, T2); Point 2s is the result of the ideal isentropic compression from point 1. The condenser operates at Phigh and condenses the vapor to its bubble point at point 3′ at (Phigh,T3′), and subsequently subcools it by ΔTsc to point 3 at (Phigh,T3 = T3′ − ΔTsc). A Joule-Thomson expansion (at constant enthalpy) then takes the subcooled liquid to a vapor−liquid mixture at point 4 in the evaporator at P = Plow, where it evaporates to its dew point at point 1′ to complete the cycle. The processes are represented by straight lines in Figure 2, and the heat transfer in the condenser and the evaporator are respectively equal to the lengths of the corresponding horizontal lines, h2−h3 and h1−h4. The efficiency of the VCRC is measured by its coefficient of performance (COP), the ratio of the evaporator heat transfer to the compressor work, given by

COP =

h1 − h4 h2 − h1

(1)

Typical VCRC design parameters are (Tlow ≡ T1′, Thigh ≡ T3′, ΔTsh, ΔTsc, ω). Under these specifications, simulation of the VCRC of Figure 2b and calculation of its COP requires the following property calculations and thermodynamic process simulations: 1. calculation of the dew point pressure Plow at T1′; 2. calculation of the properties at point 1, specified by Plow and T1 = T1′ + ΔTsh; 3. calculation of T2s and other properties at point 2s for the isentropic process 1−2s occurring from the initial state point 1 to the pressure Phigh 4. calculation of T2 and other properties at point 2 for the process 1−2 from the initial state point 1 to the pressure Phigh, occurring with the enthalpy change given by

2. THERMODYNAMIC BACKGROUND 2.1. The Vapor-Compression Refrigeration Cycle. The single-stage vapor-compression refrigeration cycle (VCRC) is illustrated schematically in Figure 1, and by means of a

h2 − h1 =

h2s − h1 ω

(2)

5. calculation of the bubble point pressure Phigh at T3′; 6. calculation of the properties at point 3, specified by Phigh and T3 = T3′ − ΔTsc; The property calculations and process simulations required to simulate the VCRC of Figure 2b thus include one dew point pressure calculation, one bubble point pressure calculation, two calculations of properties for specified P and T, one simulation of an ideal isentropic compression process for a specified P change from an initial state point, and one simulation of a process undergoing a specified enthalpy change from an initial state point. Molecular simulation methodologies for accomplishing these tasks are described in section 3. 2.2. Conventions for Enthalpy and Entropy. We express the molar enthalpy as a sum of ideal-gas and residual contributions, via

Figure 1. Schematic diagram of the equipment employed in a singlestage vapor-compression refrigeration cycle.

pressure-enthalpy (P−h) diagram (where P is the pressure and h is the molar enthalpy) in Figure 2; in Figure 2a, for a pure fluid the vapor−liquid equilibrium isotherms (represented by the dashed lines at Tlow and Thigh) are horizontal, whereas in the typical case of a zeotropic mixture in Figure 2b, the corresponding isotherms for T1′ and T3′ are at a negative angular displacement.

h(T , P) = hIG(T ) + hr(T , P)

(3)

IG denotes the ideal-gas value, R is the universal gas constant, superscript r denotes a residual property, and 3259



Article

Figure 2. P−h diagram for a VCRC for a pure refrigerant working fluid (a) and for a zeotropic mixture (b). The numbers indicate the state points of cycle of Figure 1, and primes indicate points on the vapor−liquid saturation curve of the refrigerant: the 1−2 process is compression, 2−3 is condensation, 3−4 is expansion, and 4−1 is evaporation (1′−1 is the superheating of vapor and 3′−3 the subcooling of liquid). Point 2s, connected to point 1 by a dashed line, denotes the final state of an ideal isentropic compression from point 1.

hIG(T , P) ≡ hIG(T ) = hIG(Tref , Pref ) +

∫T

T

P calculations are required for the simulation of the VCRC of section 2.1.

cpIG(T ) dT

ref

(4)

Table 1. Important Types of Phase Equilibrium Problems for Fluid Mixtures Involving T,P and Phase Composition Specificationsa

cp is the molar isobaric heat capacity, defined by cp =

⎛ ∂h ⎞ ⎜ ⎟ ⎝ ∂T ⎠ P

(5)

Here, we use Tref = 298.15 K and Pref = 1 bar, and set hIG(Tref, Pref) = ΔHf(298.15), where ΔHf(298.15) is the standard enthalpy of formation at 298.15 K. Values of ΔHf(298.15), hIG(T) − hIG(Tref, Pref) and cIG P (T) can be calculated essentially exactly using quantum mechanical approaches or obtained from experimental measurements, and values are also available in standard thermochemical compilations.45−47 A curve of constant enthalpy in the plane of two other variables, such as (P,T), defines an isoenthalpic curve, HC. Such curves are useful for graphically determining the final state of an adiabatic process (occurring at Δh = 0) from an initial (P,T) state. Similarly to eq 3, the molar entropy, s, may be expressed as s(T , P) = s IG(T ) + s r(T , P)

∫T

T

(6)

ref

cpIG(T ) T

dT

problem type

T,P,w h,P,w T,x P,x T,y P,y

x,y T,x,y P,y T,y P,x T,x

(T,P) flash isenthalpic (h,P) flash bubble-P bubble-T dew-P dew-T

3. COMPUTATIONAL METHODS 3.1. Simulations for Properties at Specified (P,T) and for VLE. We employ the conventional NPT MC algorithm for calculating single-phase properties,9,10 which uses a fixed number of particles, N, and specified values of P and T. The total system enthalpy H is calculated in such a simulation from

(7)

s is related to h by

1 s = (h − μ) T

calculated parameters

a x = (x1, x2, ..., xc), where c is the number of components, denotes the liquid phase mole fraction vector. y = (y1, y2, ..., yc), where c is the number of components, denotes the vapor phase mole fraction vector. w = (w1, w2, ..., wc), where c is the number of components, denotes the overall mole fraction vector.

where s IG(T ) = s IG(Tref , Pref ) +

specified parameters

H = HIG + < ext + PV − NkT

(9)

ext

where < is the intermolecular configurational energy and k is the Boltzmann’s constant. We calculate the total entropy S for binary mixtures from eq 6, using 1 S = (H − N1μ1 − N2μ2 ) (10) T where subscript i denotes a species, and

(8)

where μ is the molar chemical potential, which may also be separated into its IG and residual parts, similarly as in eqs 3 and 6. 2.3. Dew and Bubble Points. Calculations involving specifications of the composition of the liquid or vapor phase of a vapor−liquid equilibrium (VLE) state of a fluid mixture are important in refrigeration and other applications. Six common types of VLE problems are listed in Table 1. Bubble-P and dew-

⎞ ⎛ NP i μi = μi IG + kT ln⎜ ⎟ + μie ⎝ N1 + N2 ⎠ 3260

(11)


Journal of Chemical & Engineering Data μi IG = hiIG − TsiIG

Article

also be used to generate an HC by selecting a set of P values and using the algorithm to calculate the resulting T values satisfying ΔH = 0 from one of the points. We use this approach in the examples of section 4.1. Another means of calculating the T resulting from a finite P change for a process occurring from an initial state (P1, T1) to a pressure Pf for a specified enthalpy change ΔH is by numerically solving for Tf in the equation

(12)

and μei is the excess chemical potential of species i with respect to an ideal gas mixture at specified (P,T), which is calculated by the conventional Widom test-particle-insertion method. Although the ideal-gas quantities can be taken from any thermochemical data compilation, 45−47 since we make comparisons with REFPROP6 calculations, we use the same equations for the relevant ideal-gas quantities as those used by REFPROP. We employ the Gibbs Ensemble Monte Carlo (GEMC) algorithm9,10 for pure-fluid VLE calculations. T is specified and the vapor pressure P(T) and the densities of the coexisting phases are calculated. For a binary mixture, there are two typical extensions of the method: NPT GEMC, for which T and P are specified; and NVT GEMC, for which T and the total system volume V are specified. For a binary mixture, we use instead the NPT Reaction Gibbs Ensemble Monte Carlo (RGEMC) algorithm,48,49 a variant of GEMC that improves mixture VLE results by incorporating methodology to ensure exact agreement with known experimentally measured pure component vapor pressure data, similar in spirit to the approach taken by traditional thermodynamic models for calculating mixture VLE properties.50,51 3.2. Simulating a Process Occurring for a Specified Enthalpy or Entropy Change and Calculating HC. We remark in passing that NPH MD simulation algorithms in the special case when H consists only of kinetic and intermolecular energy terms have been developed by Ray and Graben52 and by Kioupis and Maginn.53 Since such an enthalpy excludes additional ideal-gas contributions arising from the internal partition function for a real fluid, these algorithms are not applicable to the real fluids considered here. Escobedo and Chen11 calculated HC using NPT MC simulation by means of an integration scheme involving a first-order Taylor expansion of h in P and β = 1/(kT), incorporating IG contributions as required. This integration scheme can also in principle be used for simulating a process occurring for a finite pressure change, although it would be time-consuming. The NPH MC algorithm of Smith and Lsal22−24 includes all terms in H and allows the direct calculation of the final T due to a finite pressure change from an initial state within a single simulation run. The methodology has been described in our previous papers,22,23 and only a brief summary is given here. The algorithm generates a Markov chain to simulate the properties; for a single-phase system this consists of three types of state transitions (moves), involving particle positions and orientations, system volume, and inverse temperature β = 1/kT. The transition probabilities for particle position, orientation, and volume change moves use standard approaches,9,10 and the β change transition probability is given in the original references quoted above. If the final state point lies in a vapor−liquid region, simulations are performed in two boxes that correspond to the vapor and liquid phases, and the algorithm can thereby be considered to be determining the temperature for an isenthalpic flash calculation at the specified pressure (see Table 1). The Markov chain in such a two-phase region consists of additional moves involving particle transfers between the two boxes. Since an HC can be viewed as a curve of (P,T) points for which the total enthalpy change ΔH is set to 0 for a process occurring between any two points, the NPH MC algorithm can

H(Pf , Tf ) = H(P1 , T1) + ΔH

(13)

by an iterative method that performs a sequence of NPT MC simulations at states (Pf,Tn) that converges to the solution. We have developed25 an NPS MC method for simulating a process occurring for a specified entropy change, which combines the use of the NPH MC algorithm with an iterative scheme that incorporates calculation of the chemical potential in eq 8; see the indicated reference for details. For the isentropic process in step 3 of the VCRC simulation procedure in section 2.1, to reduce the computer time, we instead use a numerical method analogous to that for solving eq 13, for the equation S(Pf , Tf ) = S(P1 , T1)

(14)

We solve eqs 13 and 14 simultaneously to determine T2s and T2 in steps 3 and 4 of the VCRC of section 2.1 as follows. We first perform (in parallel) NPT simulations at a series of (Pf, Tn) states over an interval including initial approximations to the values of T2s and T2, at which s(Pf, Tn) and h(Pf, Tn) are calculated. Interpolating equations h̃(T) and s̃(T) are then fitted to the resulting (sn, Tn) and (hn, Tn) values. s̃(T) is then used to find T2s satisfying eq 14, and h̃(T2s) is used to calculate h2 from eq 2. Finally, h̃(T) is used to calculate T2 from eq 13. By performing the NPT simulations in parallel, this procedure allows the rapid determination of T2. 3.3. Dew and Bubble Point Calculations. Four important types of dew and bubble point problem are listed in the final four rows of Table 1. Several molecular simulation approaches have been developed for bubble-P problems. Calculations for binary and ternary mixtures were first considered by Vrabec and Fischer,28 who developed the NPT + test particle method, and the related Grand Equilibrium method was later proposed by Vrabec and Hasse. 29 Escobedo33−35 proposed an algorithm related to a method of Kofke for flash calculations59 to solve bubble-P problems. Ungerer et al.36,37 solved bubble-P problems by performing “virtual” GEMC flash calculations in which the liquid phase composition is held fixed. This approach was applied to several systems by Ferrando et al.39,40 In addition to flash problems, Escobedo’s approach can in principle also solve all four types of dew and bubble point problems. We have recently developed an algorithm41−44 to solve these problems by implementing a Newton−Raphson (NR) method to solve the equation

z(P , T ) = z *

(15)

where z is a species mole fraction in the vapor (z ≡ y for a dew point) or liquid (z ≡ x for a bubble point) phase, and superscript asterisk (∗) denotes its specified value. Either P or T is also specified, and to obtain the remaining thermodynamic variable, Y (T or P) and the composition of the other (liquid or vapor) phase, we use the NR method on eq 15. We start from an initial estimate, (Y0, T) or (P, Y0), corresponding to a state at phase equilibrium; each iteration corresponds to a point on the 3261



Article

For mixtures, the standard Lorentz−Berthelot rules were used for the cross LJ parameters. We have previously shown that the RGEMC algorithm,48,49 used within the algorithm for calculating dew and bubble points described, produces excellent VLE results for mixtures using these rules (the RGEMC methodology can be used with any combining rule). Here, we used REFPROP pure-component vapor pressure data in the RGEMC implementation. For the NPH MC simulations, for pure fluids we used N = 256 molecules for single-phase simulations and N = 364 for two-phase simulations. For single-phase NPT mixture simulations we used N = 864 and for RGEMC VLE calculations we generally used N = 2304. In all simulations, we used cubic boxes with the minimum image convention and periodic boundary conditions. For the LJ interactions, long-range corrections for the configurational energy and the pressure were included,9 assuming that the radial distribution functions are unity beyond a cutoff radius of half the box length; the dipolar long-ranged corrections were treated by the reaction field method, with the dielectric constant set to infinity for the liquid state and unity for the gaseous state. To prevent overlap of the dipole sites a hard sphere cutoff diameter of 0.5L was employed. The simulations were organized in cycles consisting of three (single-phase simulations) or four (two-phase simulations) types of attempted state transitions (moves): nD particle displacement moves (translation and rotation, chosen with equal probability); nV volume change moves; (where relevant) nT particle transfer moves; and (where relevant) nβ inverse temperature change moves. The types of moves were selected at random with fixed probabilities, chosen so that an appropriate ratio of moves was obtained. For single-phase NPH MC simulations, for each cycle, we set (nD, nV, nβ) = (N,1,1); for two-phase NPH MC simulations, (nD, nV, nT, nβ) = (N,2,10N,1). In the case of single-phase NPT simulations, we set (nD, nV) = (N,1), and for two-phase NPT RGEMC simulations, we set (nD, nV, nT) = (N,1,10N). The maximum step-size parameters governing the acceptance ratios for translational and rotational moves and for volume changes were adjusted to yield an acceptance ratio of approximately 30 %. The corresponding parameter governing the acceptance ratio for inverse temperature change moves was adjusted to yield an acceptance ratio of approximately 50 %. Each simulation was performed from an initial molecular configuration (typically an FCC crystalline lattice) and an initial equilibration period was followed by a production run during which statistics were collected over blocks of cycles. For all (single-phase and two-phase) NPH MC simulations, following an initial equilibration period of 60 000 cycles we generated 100 000 cycles in a production run and used 1000 cycles per block. In the case of mixture simulations, for single-phase NPT MC simulation, the equilibration period was 3·104 cycles and the production run was 1.5·105 cycles, with 30 000 cycles per block. For RGEMC VLE simulations, the equilibration period was 3·105 cycles and the production run was 1.5·106 cycles; 6000 cycles per block were used to collect statistics. In addition, we monitored the convergence profiles of the thermodynamic quantities, in order to keep the development of the system under careful control.62 Uncertainties, u(x) in the production run results were calculated using

coexistence curve, which is calculated by means of a RGEMC (flash) simulation. The NR iterations are given by Y (n + 1) = Y (n) −

z(n) − z*

( ∂∂Yz )σ (n)

(16)

where n denotes the iteration number, (z(n),Y(n)) are states on the coexistence curve (denoted by the subscript σ), and the partial derivatives (∂z/∂Y)σ are calculated along this curve as follows during each RGEMC simulation by means of fluctuation formulas.54−58 In the case of a bubble-T (z ≡ x) or dew-T (z ≡ y) calculation, Y ≡ T, and ⎛ ∂⟨z⟩ ⎞ ⎛ ∂z ⎞ ⎜ ⎟ ≡ ⎜ ⎟ = ⟨z⟩⟨H ⟩ − ⟨zH ⟩ = ⟨z⟩⟨Ĥ ⟩ − ⟨zĤ ⟩ ⎝ ∂T ⎠σ ⎝ ∂T ⎠ σ (17)

where Ĥ = < ext + PV is the intermolecular configurational enthalpy, and the brackets denote averages obtained during the course of the RGEMC simulation at the current iteration. For a bubble-P (z ≡ x) or dew-P calculation (z ≡ y), Y ≡ P, and ⎛ ∂⟨z⟩ ⎞ ⎛ ∂z ⎞ ⎜ ⎟ ≡ ⎜ ⎟ = β[⟨z⟩⟨V ⟩ − ⟨zV ⟩] ⎝ ∂P ⎠σ ⎝ ∂P ⎠ σ

(18)

3.4. Molecular Models and Simulation Details. We consider the case of R32/R134a mixtures, for which the molecules are modeled using force fields based on a rigid 2center homonuclear Lennard-Jones dipolar (2CLJD) potential model, composed of two identical Lennard-Jones (LJ) sites at distance L apart, and an axial point dipole of moment D embedded in the geometric center of the molecule and oriented along the molecular axis: u 2CLJD(r12, ω 1, ω 2) = u 2CLJ(r12, ω 1, ω 2) + uD(r12, ω 1, ω 2)

(19)

where r12 is the center-to-center vector between molecules 1 and 2 and the vectors ω1 and ω2 represent the orientation of molecules 1 and 2. u2CLJ denotes the contribution of the four LJ sites, ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ u 2CLJ(r12, ω 1, ω 2) = ∑ ∑ 4ϵ⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎢⎣⎝ rab ⎠ ⎝ rab ⎠ ⎥⎦ a=1 b=1 2

2

(20)

where rab is the distance between LJ sites a and b on molecules 1 and 2, respectively, and σ and ϵ are the LJ size and energy parameters, respectively. uD is the dipolar contribution: uD(r12, ω 1, ω 2) =

3(D1r12)(D2 r12) ⎤ 1 ⎡ ⎢D1D2 − ⎥ 3 2 4π ϵ0r12 ⎣ r12 ⎦ (21)

where ϵ0 is the permittivity of free space, r12 = |r12, and D1 and D2 are the dipoles of molecules 1 and 2. There are four force field parameters for each pure fluid, for which we used the values of Lisal60 and of Stoll et al.61 in section 4.1. The parameters in the latter case were obtained by fitting to the entire experimental VLE curves. For the former case, the parameters were fitted to a smaller data set for each fluid, consisting of the experimental critical temperature, Tc, the saturated liquid density at T = 0.75Tc, and the steepness of the vapor pressure curve. 3262



Article

Table 2. Comparison of Simulation Results for R134a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −900 kJ mol−1a P/MPa 0.01 0.03 0.05 0.07 0.10 0.20 0.30 0.50 1.00 2.00 3.00 5.00 7.00 10.00 15.00 20.00 30.00 40.00 50.00 60.00 70.00 100.00 110.00 130.00 150.00 300.00 400.00 500.00

ρ/mol m−3

T/K RFP

SFF

LFF

RFP

SFF

LFF

375.55 375.75 375.95 376.15 376.44 377.44 378.42 380.40 385.31 395.00 404.47 422.06 436.67 452.05 NA NA NA NA NA NA NA NA NA NA NA NA NA NA

376.10 376.40 376.49 376.72 377.12 378.04 379.33 381.56 387.02 397.46 407.06 424.63 438.95 454.37 467.37 474.15 479.89 481.33 481.62 480.61 478.76 471.80 469.32 462.92 456.81 405.64 370.01 332.85

376.10 376.30 376.33 376.62 376.96 377.80 378.53 380.46 385.06 394.02 402.88 419.50 434.15 449.70 463.32 469.74 474.87 476.28 475.97 474.12 472.37 464.23 461.28 455.10 448.21 392.76 353.66 313.85

3.21 9.63 16.07 22.52 32.23 64.84 97.85 165.05 340.27 723.19 1151.2 2126.3 3167.1 4557.1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA

3.19 9.58 15.98 22.40 32.04 64.54 97.32 164.39 340.45 726.38 1160.4 2149.7 3189.5 4581.1 6110.3 7091.7 8357.9 9185.8 9811.8 10318. 10743. 11781. 12064. 12559. 12996. 15339. 16472. 17412.

3.19 9.58 15.97 22.38 32.02 64.45 97.24 164.03 337.76 718.48 1142.4 2136.2 3240.5 4736.0 6399.1 7427.8 8730.7 9598.5 10230. 10709. 11153. 12182. 12464. 12967. 13407. 15716. 16817. 17747.

All points are in the gaseous phase. NA denotes “not available”, due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in the densities. a

⎛ Nb (x − x )2 ⎞1/2 i ̅ ⎟ u(x) = 1.96⎜⎜∑ ⎟ N ( N − ⎝ i = 1 b b 1) ⎠

For comparison purposes, when using the REFPROP software we set the reference state temperature to 298.15 K, the reference state enthalpy to the appropriate ΔHf value, and the reference state pressure to 10−7 MPa. We calculated HC for each refrigerant at representative values of h by performing NPH MC simulations at a sequence of P values along the isoenthalp to calculate corresponding T and density values. The numerical results are given in Tables 2 to 8, and graphs are displayed in Figures 3 and 4. We first note that simulation results are available in the tables and figures at states outside the range of the REFPROP data predictions, illustrating a general advantage of the fundamentally based molecular simulation approach over the empirical macroscopically based approach of REFPROP. The results indicate that the predictions of both force fields are in good agreement with those of REFPROP, particularly in single-phase regions, with the StollFF results being generally superior to those of the LisalFF. The most important quantity is the predicted temperature, since this has an important effect on the density predictions. Table 8 indicates that the StollFF-predicted temperatures in single-phase regions deviate from those of REFPROP by a maximum of 0.7 % for R134a and 2.4 % for R143a; the corresponding LisalFF maximum deviations are 6.9 % and 4.6 %, respectively. In two-phase regions, the StollFF maximum T deviations are 2.6 % for R134a and 4.2 % for R143a, and the corresponding LisalFF maximum deviations are 7.1 % and 7.5 %. The maximum single-phase density deviations

(22)

where Nb is the number of blocks, xi is the average value of x in block i, and x̅ is the average value of x over all blocks. This yields an approximate 95 % uncertainty value for the calculated value of x ≡ x.̅ For the NPH MC calculations we used in-house code running on a Mac Pro 12-core (two Intel Xeon 6-core 3.06 GHz processors) using the Intel Fortran Compiler. For the VCRC simulation, we used in-house code running on an Intel Xeon X5365 3 GHz processor using the IntelClusterStudio IFORT v.11.0 2008 1105 compiler.

4. RESULTS AND DISCUSSION 4.1. Isoenthalps for R134a and R143a Using Two Different Force Fields. We consider the calculation of HC for R134a (1,1,1,2-tetrafluoroethane) and R143a (1,1,1-trifluoroethane) using the NPH MC algorithm described in section 3.2. We compare results using the two molecular force fields described in section 3.4, one due to Lisal et al.60 (referred to as the LisalFF) and the other due to Stoll et al.61 (referred to as the StollFF). We compare the simulation results with pseudoexperimental results generated using the REFPROP software package.6 We used ΔHf(298.15) = −907.1 kJ mol−1 for R134a45 and ΔHf(298.15) = −748.7 kJ mol−1 for R143a.47 3263



Article

Table 3. Comparison of Simulation Results for R134a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −910 kJ mol−1a P/MPa 0.01 0.03 0.05 0.07 0.10 0.20 0.30 0.50 1.00 1.30 1.50 1.80 2.00 2.50 5.00 7.00 10.00 15.00 20.00 30.00 40.00 50.00 60.00 70.00 100.00 110.00 130.00 150.00 300.00 400.00 500.00

ρvap/mol m−3

T/K

ρliq/mol m−3

RFP

SFF

LFF

RFP

SFF

LFF

263.02 263.68 264.35 265.01 266.00 269.31 273.82 288.88 312.54 322.61 328.38 336.05 340.63 350.73 384.29 395.20 402.44 407.96 410.59 412.44 412.19 410.89 408.95 406.58 NA NA NA NA NA NA NA

263.34 264.18 265.01 265.73 267.05 270.92 274.51 281.41 305.24 NC 329.56 337.32 341.67 350.38 386.08 395.97 402.75 407.97 411.19 412.90 412.59 411.34 409.52 407.35 398.71 395.95 389.40 382.81 327.98 290.39 253.30

263.31 263.92 264.42 264.88 265.93 268.62 271.45 276.82 294.10 314.92 320.92 328.53 333.11 342.83 381.55 390.62 397.14 401.30 403.69 404.67 404.12 402.15 400.22 397.27 387.63 384.46 377.00 369.61 309.72 267.49 224.38

4.59 13.82 23.11 32.48 46.68 95.20 144.76 238.33 482.42 639.80 750.70 927.96 1054.80 1412.20

4.57 13.74 22.98 32.28 46.41 94.63 144.89 252.66 524.74 NC 762.01 941.39 1075.85 1484.41

4.57 13.72 22.93 32.25 46.30 94.19 143.88 294.02 534.66 654.03 760.48 942.16 1071.3 1448.2

RFP

SFF

LFF

12668. 12161. 11264. 10828. 10557. 10167. 9912.2 9273.1 5473.4 6778.3 7743.0 8646.3 9238.1 10051. 10631. 11089. 11473. 11806. NA NA NA NA NA NA NA

11439. NC 9783.7 9600.3 9421.5 8993.2 5706.9 6850.6 7731.5 8614.7 9237.4 10076. 10682. 11172. 11600. 11956. 12842. 13100. 13552. 13961. 16196. 17315. 18293.

12803. 10411. 10410. 10377. 10139. 9690.0 6193.1 7260.9 8219.5 9067.4 9671.6 10502. 11123. 11609. 12043. 12396. 13270. 13522. 13972. 14374. 16587. 17683. 18652.

Italics denote two-phase points. NC denotes “not calculated”; NA denotes “not available” due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in single-phase densities, and 2 % in two-phase densities.

a

4.2. Molecular Simulation of a VCRC. We consider the simulation of a VCRC involving a 30 mass % R32/R134a mixture for the cycle described in section 2.1, according to the procedures described in section 3. The VCRC specifications are a condenser outlet temperature of Thigh = 328.15 K, an evaporator inlet temperature of Tlow = 278.15 K, condenser subcooling and evaporator superheating of 5 K, and a compressor of 80 % isentropic efficiency. We used the simulation protocols and the force fields of Stoll et al.61 for the mixture components, described in section 3.4. We used pure-component vapor pressure data from REFPROP within the RGEMC algorithm, used for the calculation of the dew and bubble points. To facilitate comparison with REFPROP results, we used the REFPROP equations for ideal-gas properties, and employed the default reference state used by REFPROP. The results for selected properties are given at the points of the cycle in Tables 9 and 10, and the P−h cycle diagram is shown in Figure 5. The dew-P and bubble-P simulation results of Table 9 were achieved in three or four iterations of the procedure described in section 3.3. They agree very well with those of REFPROP, the largest discrepancy of about 3 % occurring for ρvap at the dew point 1′. We implemented the bubble-P method of Ungerer et al.36,37 for point 3′ by

of the StollFF results from the REFPROP values are 4.3 % for R134a and 2.3 % for R143a, and the corresponding LisalFF deviations are 13.1 % and 3.9 %. For two-phase states, the StollFF maximum density deviations are 7.3 % for R134a and 6.9 % for R143a, and the corresponding LisalFF maximum deviations are 14.4 % and 12.4 %. The maximum density deviations of the simulation results from those of REFPROP generally occur at pressures in Tables 2 to 7 for which the maximum temperature deviations occur. The primary reason for the relative superiority of the StollFF over the LisalFF results is readily ascertained by examining their VLE predictions, which are displayed in Figures 6 and 7 of Stoll et al.61 These figures show that the StollFF predicts the VLE properties of the two refrigerants much more accurately than does the LisalFF. This is not unexpected, in view of the fact that the StollFF parameters were obtained by fitting to experimental data over a range of VLE results, whereas the LisalFF used only the critical temperature, Tc, the saturated liquid density at T = 0.75Tc, and the steepness of the vapor pressure curve for each fluid. Our results demonstrate that the advantages of requiring a smaller experimental data set to determine the force field parameters is offset by a decline in the accuracy of the resulting predictions. 3264



Article

Table 4. Comparison of Simulation Results for R134a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −925 kJ mol−1a P/MPa 0.01 0.03 0.05 0.07 0.10 0.12 0.15 0.18 0.20 1 2 3 5 7 10 15 20 30 40 50 60 70 100 110 130 150 300

ρvap/mol m−3

T/K

ρliq/mol m−3

RFP

SFF

LFF

RFP

SFF

LFF

RFP

SFF

LFF

206.29 223.47 232.70 239.29 246.79 250.84 256.02 260.44 263.07 309.46 309.52 309.56 309.57 309.49 309.26 308.61 307.72 305.46 302.77 299.77 296.56 293.19 NA NA NA NA NA

207.07 222.39 230.34 236.28 244.05 247.30 252.59 256.81 258.39 307.37 307.63 307.85 307.66 307.38 307.14 306.72 305.96 303.70 301.21 298.36 295.27 291.99 281.31 277.78 270.23 262.66 202.53

190.10 206.34 216.54 222.24 231.44 234.94 241.11 246.62 249.05 293.27 293.20 293.08 293.11 292.69 292.55 291.31 290.12 287.11 284.09 280.41 276.62 272.83 260.58 256.36 247.76 239.13 170.40

5.88 16.45 26.54 36.39 50.90 60.45 74.66 88.76 98.13

5.90 16.90 27.25 37.46 52.36 62.34 76.85 91.77 101.34

6.50 18.06 28.77 39.40 54.93 65.16 79.95 94.56 104.28

14635. 14166. 13908. 13719. 13501. 13381. 13226. 13091. 13009. 11396. 11466. 11533. 11659. 11775. 11936. 12175. 12387. 12755. 13072. 13352. 13606. 13838. NA NA NA NA NA

14734. 14326. 14118. 13944. 13717. 13618. 13459. 13302. 13284. 11442. 11524. 11593. 11725. 11841. 12008. 12275. 12513. 12920. 13268. 13589. 13880. 14138. 14821. 15034. 15403. 15757. 17846.

15203. 14760. 14454. 14294. 13995. 13891. 13693. 13512. 13435. 12017. 12077. 12149. 12285. 12393. 12576. 12806. 13027. 13432. 13782. 14083. 14365. 14618. 15304. 15504. 15881. 16235. 18348.

Italics denote two-phase points. NA denotes “not available”, due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in single-phase densities, and 2 % in two-phase densities.

a

Table 5. Comparison of Simulation Results for R143a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −745 kJ mol−1a P/MPa 0.001 0.003 0.01 0.03 0.1 0.3 0.7 1 2 5 10 25 40 100 200 400

ρ/mol m−3

T/K RFP

SFF

LFF

RFP

SFF

LFF

343.33 343.36 343.44 343.67 344.50 346.82 351.42 354.81 365.81 395.08 425.00 447.29 450.03 435.85 NA NA

343.72 343.89 343.71 343.83 344.75 347.40 351.95 355.72 366.49 395.26 425.20 447.90 451.50 439.15 404.37 325.79

343.88 343.69 343.80 344.09 344.77 346.80 351.10 354.45 364.34 391.84 421.66 444.74 447.45 433.44 396.14 312.14

0.35 1.05 3.51 10.53 35.27 107.23 257.00 374.58 800.28 2355.5 4899.7 8231.3 9664.3 12263. NA NA

0.35 1.05 3.48 10.47 35.01 106.42 255.62 372.27 789.43 2314.1 4771.8 8004.2 9430.5 12111. 14379. 17085.

0.35 1.05 3.48 10.45 35.03 106.41 254.88 371.23 788.48 2326.8 4903.9 8345.7 9792.6 12525. 14787. 17471.

a All points are in the gaseous phase. NA denotes “not available”, due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in single-phase densities.

performing 105 cycles of a virtual NPT RGEMC simulation keeping the liquid phase mole fraction fixed at xliq = 0.4567, followed by 1.5·105 cycles of an NVT RGEMC simulation from the final configuration. The results agree reasonably well with

those of our algorithm, nearly falling within their mutual uncertainties. The results indicated as 1′(NVT) and 3′(NVT) are supplementary “verification” NVT RGEMC simulations that were conducted at the points 1′ and 3′. We initiated the 3265



Article

Table 6. Comparison of Simulation Results for R143a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −755 kJ mol−1a P/MPa 0.001 0.003 0.01 0.03 0.3 0.7 1 2 5 10 25 40 100 200 400

ρvap/mol m−3

T/K

ρliq/mol m−3

RFP

SFF

LFF

RFP

SFF

LFF

208.13 208.31 208.91 210.50 251.80 277.13 289.48 316.90 354.92 366.46 372.80 371.58 353.01 NA NA

208.34 208.34 208.96 210.52 246.21 277.12 290.58 314.95 355.01 367.64 372.72 374.02 357.63 319.36 236.65

208.28 208.47 208.79 209.97 233.02 269.68 283.52 309.27 351.12 362.69 368.06 366.69 348.00 306.35 214.77

0.58 1.74 5.80 17.48 159.15 367.65 532.92 1183.3

0.58 1.73 5.76 17.36 164.63 367.96 523.4 1198.4

0.58 1.72 5.76 17.35 175.77 378.89 536.47 1194.53

RFP

SFF

LFF

12998. 12025. 11485. 9996.4 6967.1 8639.3 10455. 11416. 13474. NA NA

13107. 11560. 10989. 9928.3 6819.2 8445.1 10300. 11206. 13419. 15490. 18165.

14615. 12097. 11602. 10356. 7239.9 8855.5 10643. 11626. 13848. 15912. 18589.

Italics denote two-phase points. NC denotes “not calculated”; NA denotes “not available”, due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in single-phase densities, and 2 % in two-phase densities.

a

Table 7. Comparison of Simulation Results for R143a Using Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 with Those of REFPROP (RFP)6 at h = −770 kJ mol−1a P/MPa 0.01 0.03 0.05 0.07 0.10 0.12 0.15 0.7 1 2 5 10 25 40 100 200 400

ρvap/mol m−3

T/K

ρliq/mol m−3

RFP

SFF

LFF

RFP

SFF

LFF

RFP

SFF

LFF

186.77 203.20 212.06 218.41 225.63 229.54 234.55 254.48 254.43 254.24 253.59 252.32 247.60 242.05 NA NA NA

185.29 200.19 210.59 215.86 222.33 227.34 231.42 258.35 258.69 258.31 257.96 257.11 252.80 247.78 224.59 180.64 111.14

175.27 193.04 202.51 208.36 215.77 220.17 225.10 244.38 244.39 244.23 243.46 242.07 236.87 230.91 203.64 154.50 NC

6.51 18.17 29.28 40.11 56.04 66.53 82.13

6.58 18.51 29.50 40.95 57.13 67.58 83.52

6.96 19.16 30.72 42.26 58.67 69.27 86.19

15090. 14595. 14320. 14120. 13888. 13760. 13594. 12919. 12936. 12989. 13143. 13376. 13962. 14443. NA NA NA

15281. 14816. 14478. 14323. 14074. 13900. 13765. 12725. 12739. 12799. 12968. 13216. 13838. 14359. 15895. 17699. 19061.

15499. 14937. 14608. 14426. 14191. 14015. 13870. 13286. 13312. 13362. 13533. 13774. 14374. 14884. 16416. 18245. NC

Italics denote two-phase points. NC denotes “not calculated”; NA denotes “not available”, due to REFPROP limitations; simulation uncertainties (confidence level 0.95) are 1 K in T, 1 % in single-phase densities, and 2 % in two-phase densities.

a

simulations in each case from the final configurations of our dew-P and bubble-P simulations, from which production runs of 1.5·105 cycles were carried out. The resulting 1′(NVT) and 3′(NVT) results for all properties are seen to agree reasonably well with our dew-P and bubble-P results. We note that the P uncertainties in our dew-P and bubble-P results were calculated from the discrepancies in exactly matching the respective vapor or liquid mole fractions, which are quite small. The uncertainties of the NVT simulations are larger, because of the pressure fluctuations inherent in the algorithm. Table 10 compares our simulation results with those of REFPROP at the single-phase points of the VCRC. The results for the properties at points 1 and 3 were calculated by means of standard NPT simulations using the methodology described in section 3. For point 1, the simulated value of h agrees with the

REFPROP value to within 0.7 %, s to within 4.0 %, and the (vapor) density to within 2.5 %. For point 3, the simulated value of h agrees with the REFPROP value to within 0.5 % and the (liquid) density to within 3.8 %. The simulation results at points 2s and 2 were obtained by solving for the temperatures resulting from the isentropic compression 1−2s, and from the specified enthalpy process 1− 2 according to the procedure of section 3.2. We first performed NPT simulations in parallel at the calculated bubble point pressure P2 ≡ Phigh with y = 0.4567 and nine temperatures, increasing from 335.15 K in 5 K intervals. We then fitted the resulting h and s values to T using weighted linear regression (using the squared inverses of the simulation uncertainties as weights), yielding 3266



Article

Table 8. : Maximum Deviations of Simulation Results in Tables 2 to 7 from Those of REFPROP6 Using the Force Fields of Stoll et al. (SFF)61 and of Lisal et al. (LFF)60 fluid h/kJ mol R134a -900 R134a -910 R134a -925 R143a -745 R143a -755 R143a -770

T/K

region −1

one-phase two-phase one-phase two-phase one-phase two-phase one-phase two-phase one-phase two-phase one-phase two-phase

ρvap/mol m−3

ρliq/mol m−3 SFF

LFF

4.3 7.3 2.2 2.1

13.1 14.4 5.6 4.2

2.3 6.9 1.5 1.5

3.9 12.4 2.8 2.7

SFF

LFF

SFF

LFF

0.6

0.6

1.1

3.9

0.6 2.6 0.7 1.8 0.8

2.3 5.9 6.9 7.1 0.8

1.9 8.8

1.1 10.8

3.3 2.8

10.5 2.1

1.3 4.2 2.4 1.5

1.4 7.5 4.6 6.2

0.6 3.4

0.9 10.4

2.1

6.9

Figure 4. Temperature, T, as a function of density, ρ, along representative isoenthalps for the refrigerants R134a and R143a. The enthalpy values are indicated, using the reference state convention discussed in the text. The points are simulation results, with squares representing those obtained using the force field of Stoll et al.61 and triangles corresponding to the force field of Lisal et al.60 The curves are the results from REFPROP.6 Single-phase points are open symbols and two-phase points are filled. The orange colors on the curve and the filled points indicate two-phase state points. The complete REFPROP vapor−liquid phase envelope is shown for reference, with the critical point indicated by a star symbol. The statistical uncertainties of the simulation results lie within the symbols.

h ̃ = −5469.370 + 120.842T

(23)

s ̃ = 33.546 + 0.33357T

(24)

where h̃ is in J mol−1 and s̃ is in J mol−1 K−1. Setting s̃ = s1 = 154.0 in eq 24 yields T̃ 2s = 361.1 K, for which eq 23 yields h̃2s = 38167. Substituting in eq 2 yields h̃2 = 39031., and substituting in eq 24 yields T̃ 2 = 368.25 K. Verification NPT simulations at Phigh, T̃ 2s and (Phigh, T̃ 2s) yielded s2s = 154.6(1), h2 = 39110. (10), which were deemed to be satisfactory, and the simulations resulted in the property values at points 2s and 2 given in Table 10. These indicate that the respective disagreements with the REFPROP results for ρ, h, s at point 2s are 5.4 %, 0.1 %, 4.0 %, and 6.2 %, 0.2 %, 3.7 % at point 2. Overall, the simulation results of Tables 9 and 10 are seen to be in good agreement with those of REFPROP. It is important to note that the quantities h1, h2, h4, which are involved in the calculation of the COP via eq 1, agree with those of REFPROP to better than 1 %. This is reflected in Figure 5, which shows that the VCRC simulation results are in excellent agreement

Figure 3. Temperature, T, as a function of pressure, P, along representative isoenthalps for the refrigerants R134a and R143a. The enthalpy values are indicated, using the reference state convention discussed in the text. The points are simulation results, with squares representing those obtained using the force field of Stoll et al.61 and triangles corresponding to the force field of Lisal et al.60 The curves are the results from REFPROP.6 Single-phase points are open symbols and two-phase points are filled. The orange colors on the curve and the filled points indicate two-phase state points. The complete REFPROP vapor−liquid phase envelope is shown for reference, with the critical point indicated by a star symbol. The statistical uncertainties of the simulation results lie within the symbols. 3267



Article

Table 9. Comparison of Simulation (sim) and REFPROP (RFP)6 Thermodynamic Properties at the Dew (1′) and Bubble (3′) Points of Figure 5 with Those Obtained at Verification Simulations and with the Bubble-P Method of Ungerer36 at 3′a T

P/MPa

xliq

state point

K

sim

RFP

sim

1′ (dew-P) 1′(NVT) 3′ (bubble-P) 3′(NVT) 3′ (Ungerer36)

278.15* 278.15* 328.15* 328.15* 328.15*

0.501(1) 0.506(2) 2.412(2) 2.398(9) 2.368(13)

0.502

0.257(1) 0.238(1) 0.456(1) 0.456(2) 0.457(2)

2.335

ρliq/mol m−3

yvap sim

RFP

sim

RFP

sim

RFP

Nvap

0.456(1) 0.458(1) 0.599(1) 0.583(1) 0.590(1)

0.4567*

13719.(8) 13596.(8) 12292.(11) 12286.(13) 12241.(18)

13841.

252.6(1) 261.2(6) 1247.(1) 1317.(2) 1371.(13)

245.8

1613 1679 439 385 433

RFP 0.259 0.4567*

ρvap/mol m−3

0.585

12373.

1261.4

a Nvap is the average number of particles in the vapor phase simulation box. All simulations used a total of 2304 particles, except for the 3′ (Ungerer) state point, which used 2501 particles. The asterisk (∗) denotes a pre-assigned parameter value. Simulation uncertainties, calculated according to eq 22, are indicated in parentheses as the standard error in the final digits. For example, 0.463(6) indicates an uncertainty of 0.006; 12220.(9) indicates an uncertainty of 9.

Table 10. Comparison of Simulation (sim) and REFPROP (RFP)6 Thermodynamic Properties at the VCRC State Points of Figure 5a T/K state point

sim

1 2s 2 3 4

283.15* 361.1(3) 368.3(3) 323.15* NC

ρ/mol m−3

P/MPa RFP 355.77 364.29

h/kJ/mol

s/J/mol/K

sim

RFP

sim

RFP

sim

RFP

sim

RFP

0.501* 2.412* 2.412* 2.412* 0.501*

0.502* 2.335* 2.335* 2.335* 0.502*

244.9(2) 1080.(2) 1030.(2) 12256.(135) NC

238.92 1024.70 970.28 12746.00

34.71(2) 38.38(2) 39.11(2) 21.98(15) 21.98(15)*

34.958 38.357 39.207* 22.084 22.084*

154.0(1) 154.6(1) 156.6(1) NC NC

148.7 148.7* 151.0

a All NPT MC simulations used 2304 particles. The asterisk (∗) denotes a pre-assigned parameter value. NC denotes the the value was not calculated. Simulation uncertainties, calculated according to eq 22, are indicated in parentheses as the standard error in the final digits. For example, 21.98(7) indicates an uncertainty of 0.07.

the REFPROP value. These results correspondingly indicate that the force field of Stoll et al. used in our simulations provides a good description of the mixture for the calculations involved in the cycle COP simulation. 4.3. General Discussion. The methodologies employed in this paper produce results comparable to those obtained using complex multiparameter equations of state such as those embodied in REFPROP.6 Furthermore, unlike the macroscopic approach, the fundamental basis of the molecular simulation approach should allow extrapolation to thermodynamic conditions beyond those used in the parameter estimation. It also has the advantage of requiring a relatively small number of parameters (typically much fewer than 10) in comparison with the roughly 40 parameters typically required for quantitative accuracy in the macroscopic thermodynamic approach; this results in a greatly reduced need for experimental data. Even for relatively complex force fields entailing molecular internal degrees of freedom, the parameters describing these portions of the force field may be determined without the need for experimental data, by means of ab initio quantum mechanical methods (e.g., Raabe and Maginn7). Consequently, the number of parameters determined by experimental data remains relatively small. We note in passing that the NPH MC algorithm can be used to implement an isenthalpic flash calculation to calculate the VLE T for a specified P; indeed, each two-phase point in our HC calculation was determined in this manner. This has also been previously demonstrated previously demonstrated for R32.22 Although the molecular simulation methodologies demonstrated in this paper produce quantitatively accurate results, the computation times required are currently typically orders of magnitude greater than those required by the macroscopically based approach. The total calculation time for the VCRC

Figure 5. P−h diagram for the example vapor-compression refrigeration cycle discussed in the text. The numbers indicate the state points of cycle of Figure 1, and primes indicate points on the vapor−liquid saturation curve of the refrigerant: the 1−2 process is compression, 2−3 is condensation, 3−4 is expansion, and 4−1 is evaporation (1′−1 is the superheating of vapor and 3′−3 the subcooling of liquid). Point 2s denotes the final state of an ideal isentropic compression from point 1. Solid lines and points indicate the simulation results, and dashed lines indicate the REFPROP results. Fluid properties at the points of the cycle are given in Tables 9 and 10. The statistical uncertainties of the simulation results lie within the symbols.

with those of REFPROP for each stage of the cycle in the P−h diagram, with the simulated COP value being within 4.3 % of 3268



Article

parameters of the empirical equations have been determined, and the determination of force field parameters typically requires much less experimental data than is the case for the macroscopic equations of state. Although not demonstrated here, the molecular-based approach has the additional advantage that, once a force field is available, it can also be used to calculate many fluid properties without the need to determine additional parameters.2,4,5,7,63 This is unlike the macroscopic approach, which typically requires different empirical equations and corresponding experimental data to fit their parameters. Using the described hardware and software, some of the simulation algorithms demonstrated in this work take considerable amounts of computer time. This will surely improve in the future, and it can be argued that the simulation approach still has economic advantages over the multiparameter equation of state approach for new fluids, due to the reduced requirement for the expensive determination of experimental data. When improved theoretical and computational technology becomes available to determine force fields from first principles, it will be possible for molecular simulation methodology to be used for the prediction of the thermophysical properties of refrigerants and for the simulation of refrigeration cycles with minimal need for experimental data.

simulation depends on the computer hardware configuration and on the efficiency of the computer code. We first note that we have made no significant efforts to optimize the calculation times for the algorithms presented in this paper, nor did we seek to run them on the fastest currently available computers; the details of the hardware and software used in this paper are described in section 3, along with the simulation methodology. Typical computing times are shown in Table 11 for the Table 11. Typical Computer Clock Times Required for the Calculations of This Papera problem type single-phase pure-fluid NPH two-phase pure-fluid NPH binary mixture NPT binary mixture VLE (NPT or NVT) dew or bubble point

time per simulation cycle 0.04 0.27 0.11 0.42 0.42

s s s s s

total time 2h 12 h 5h 200 h 600 h

a

The hardware and software employed and the simulation methodology used are described in section 5, as is the definition of a simulation cycle.

different simulation problems considered in this paper. They range from 2 h for an NPH calculation, to about a week for a VLE calculation, to about a month for a dew- or bubble-point calculation. The VCRC simulation time is slightly longer than that required for a single dew- or bubble-point calculation, since each of these can be calculated simultaneously, and the remaining calculations are much less time-consuming. The NPT simulations at points 1′ and 2′ take only 5 h and can be calculated in parallel, as can the NPT calculations required at points 2s and 2. It may be possible to perform the dew- and bubble-point calculations similarly to that used for the latter points, which would reduce the required number of VLE calculations involved. We remark that the systems considered here involve rigid force fields. The implementation of the algorithms for flexible force fields will entail an increase in the required computing time. It can be argued that, for new candidate fluids for which experimental data is lacking, the computer time requirement is outweighed by the need of the macroscopically based approach for much more copious amounts of experimental data to determine the multiple parameters.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada (Grant OGP-1041), by the University of Guanajuato and the PROMEP program (grant support for new faculty members), and by the Grant Agency of the Academy of Sciences of Czech Republic under project No. P208/11/P392. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation in the Oil and Gas IndustryMonte-Carlo Methods; IFP Publications: Paris, 2005. (2) Ungerer, P.; Nieto-Draghi, C.; Rousseau, B.; Ahunbay, G.; Lachet, V. Molecular simulation of the thermophysical properties of fluids: From understanding toward quantitative predictions. Mol. Simul. 2007, 134, 71−89. (3) Theodorou, D. N. Progress and outlook in Monte Carlo simulations. Ind. Eng. Chem. Res. 2010, 49, 3047−3058. (4) Maginn, E. J.; Elliott, J. R. Historical perspective and current outlook for Molecular Dynamics as a chemical engineering tool. Ind. Eng. Chem. Res. 2010, 49, 3059−3078. (5) Meunier, M., Ed. Industrial Applications of Molecular Simulations; CRC Press: Boca Raton, 2011. (6) Lemmon, E. W.; Huber, M. L.; McLinden, M. O. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, version 9.1; National Institute of Standards and Technology, Standard Reference Data Program: Gaithersburg, USA, 2013. (7) Raabe, G.; Maginn, E. J. A force field for 3,3,3-fluoro-1-propenes, including HFO-1234yf. J. Phys. Chem. B 2010, 114, 10133−10142. (8) Raabe, G. Molecular modeling of fluoropropene refrigerants. J. Phys. Chem. B 2012, 116, 5744−5751. (9) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, UK, 1997.

5. CONCLUSIONS We have demonstrated the use of molecular simulation methodology for the calculation of isoenthalps for the alternative refrigerants R134a and R143a, and for the simulation of a vapor-compression refrigeration cycle involving a mixture of the refrigerants R32 and R134a. The algorithms of this paper can be readily applied for other fluids for which force fields are available. Although we have demonstrated the use of our dew- and bubble-point algorithms to calculate pressures, they can also be used to calculate the dew- or bubble-point temperature for a specified pressure; we have recently used it for this purpose elsewhere.44 The simulation results of this paper display generally excellent agreement with results obtained using the macroscopic empirically based multiparameter equation of state approach typified by the industry-standard software package REFPROP.6 Because of its fundamental nature, the molecular simulation approach can also perform calculations at thermodynamic conditions outside the range at which the 3269



Article

(10) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, USA, 2002. (11) Escobedo, F. A.; Chen, Z. Simulation of isoenthalps and JouleThomson inversion curves of pure fluids and mixtures. Mol. Simul. 2001, 26, 395−416. (12) Vrabec, J.; Kedia, G. K.; Hasse, H. Prediction of Joule-Thomson inversion curves for pure fluids and one mixture by molecular simulation. Cryogenics 2005, 45, 253−258. (13) Vrabec, J.; Kumar, A.; Hasse, H. Joule-Thomson inversion curves of mixtures by molecular simulation in comparison to advanced equations of state: Natural gas as an example. Fluid Phase Equilib. 2007, 258, 34−40. (14) Colina, C. M.; Müller, E. A. Joule-Thomson inversion curves by molecular simulation. Mol. Simul. 1997, 19, 237−246. (15) Colina, C. M.; Müller, E. A. Molecular simulation of JouleThomson inversion curves. Int. J. Thermophys. 1999, 20, 229−235. (16) Chacin, A.; Vázquez, J. M.; Müller, E. A. Molecular simulation of the Joule-Thomson inversion curve of carbon dioxide. Fluid Phase Equilib. 1999, 165, 147−155. (17) Colina, C. M.; Lsal, M.; Siperstein, F. R.; Gubbins, K. E. Accurate CO2 Joule-Thomson inversion curve by molecular simulations. Fluid Phase Equilib. 2002, 202, 253−262. (18) Lagache, M.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Prediction of thermodynamic derivative properties of fluids by Monte Carlo simulation. Phys. Chem. Chem. Phys. 2001, 3, 4333−4339. (19) Lagache, M. H.; Ungerer, P.; Boutin, A. Prediction of thermodynamic derivative properties of natural condensate gases at high pressure by Monte Carlo simulation. Fluid Phase Equilib. 2004, 220, 211−223. (20) Lustig, R. Direct molecular NVT simulation of the isobaric heat capacity, speed of sound and Joule-Thomson coefficient. Mol. Simul. 2011, 37, 457−465. (21) Lustig, R. Statistical analogues for fundamental equation of state derivatives. Mol. Phys. 2012, 110, 3041−3053. (22) Lsal, M.; Smith, W. R.; Aim, K. Direct molecular-level Monte Carlo simulation of Joule-Thomson processes. Mol. Phys. 2003, 101, 2875−2884. (23) Smith, W. R.; Lsal, M. Direct Monte Carlo simulation methods for nonreacting and reacting systems at fixed total internal energy or enthalpy. Phys. Rev. E 2002, 66, 0111041. (24) Lsal, M.; Bendova, M.; Smith, W. R. Monte Carlo adiabatic simulation of equilibrium reacting systems: The ammonia synthesis reaction. Fluid Phase Equilib. 2005, 235, 50−57. (25) Smith, W. R.; Lísal, M.; Nezbeda, I. Molecular-level Monte Carlo simulation at fixed entropy. Chem. Phys. Lett. 2006, 426, 436− 440. (26) Maillet, J.-B.; Bourasseau, E.; Soulard, L.; Clérouin, J. Constant entropy sampling and release waves of shock compressions. Phys. Rev. E 2009, 80, 021135. (27) Fazelabdolabadi, B.; Bahramian, A. Prediction of sound velocity and compressibility via molecular simulation at fixed entropy. Fluid Phase Equilib. 2010, 293, 262−264. (28) Vrabec, J.; Fischer, J. Vapour liquid equilibria of mixtures from the NPT plus test particle method. Mol. Phys. 1995, 85, 781−792. (29) Vrabec, J.; Hasse, H. Grand equilibrium: vapour−liquid equilibria by a new molecular simulation method. Mol. Phys. 2002, 100, 3375−3383. (30) Miyano, Y. Vapor−liquid equilibria from molecular simulations using the algorithm in equation of state calculations. Fluid Phase Equilib. 1998, 144, 137−144. (31) Spyriouni, T.; Economou, I. G.; Theodoru, D. N. Phase equilibria of mixtures containing chain molecules predicted through a novel simulation scheme. Phys. Rev. Lett. 1998, 80, 4466−4469. (32) Escobedo, F. A. Novel pseudo-ensembles for simulation of multicomponent phase equilibria. J. Chem. Phys. 1998, 108, 8761− 8772. (33) Escobedo, F. A. Tracing coexistence lines in multicomponent fluid mixtures by molecular simulation. J. Chem. Phys. 1999, 110, 11999−12010.

(34) Escobedo, F. A. Simulation and extrapolation of coexistence properties with single-phase and two-phase ensembles. J. Chem. Phys. 2000, 113, 8444−8456. (35) Escobedo, F. A. Molecular and macroscopic modeling of phase separation. Am. Inst. Chem. Eng. J. 2000, 46, 2086−2096. (36) Ungerer, P.; Boutin, A.; Fuchs, A. H. Direct calculation of bubble points by Monte Carlo simulation. Mol. Phys. 1999, 97, 523− 539. (37) Ungerer, P.; Boutin, A.; Fuchs, A. H. Direct calculation of bubble points for alkane mixtures by molecular simulation. Mol. Phys. 2001, 99, 1423−1434. (38) Yazaydun, A. Ö .; Martin, M. G. Bubble point pressure estimates from Gibbs ensemble simulations. Fluid Phase Equilib. 2007, 260, 195−198. (39) Ferrando, N.; Lachet, V.; Boutine, A. Monte Carlo simulations of mixtures involving ketones and aldehydes by a direct bubble pressure calculation. J. Phys. Chem. B 2010, 114, 8680−8688. (40) Ferrando, N.; Defiolle, D.; Lachet, V.; Boutin, A. Ethanoled gasoline bubble pressure deetermination and experimental and Monte Carlo modeling. Fluid Phase Equilib. 2010, 299, 132−140. (41) Smith, W. R.; Francová, M.; Kowalski, M.; Nezbeda, I. Refrigeration cycle design by molecular simulation. Proceedings of the 3rd IIR Conference on Thermophysical Properties and Transfer Processes, Boulder, CO, 2009, Paper No. 157. (42) Smith, W. R.; Francová, M.; Kowalski, M.; Nezbeda, I. Refrigeration cycle design for refrigerant mixtures by molecular simulation. Collect. Czech Chem. Commun. 2010, 75, 383−391. (43) Skvorová, M.; Nezbeda, I., Smith, W. R. Molecular-level simulation of dew-points of fluid mixtures and application to refrigerant cycle design, Proceedings of the 4th IIR Conference on Thermophysical Properties and Transfer Processes, Delft, The Netherlands 2013, Paper No. TP-045. (44) Skvorová, M.; Smith, W. R. Molecular-level simulation of bubble and dew points of fluid mixtures and application to refrigerant cycle design. Int. J. Refrig. 2014, dx.doi.org/10.1016/j.ijrefrig.2014.02.007. (45) Frenkel, M.; Marsh, K. N.; Kabo, G. J.; Wilhoit, R. C., Roganov, G. N. Thermodynamics of Organic Compounds in the Gas State; Thermodynamic Research Center: College Station, TX, 1994. (46) Chase, M., Jr. NIST-JANAF Thermochemical Tables; J. Phys. and Chem. Reference Data Monograph No. 9, Am. Chem. Society, Am. Inst. Physics; National Institute of Standards and Technology: Gaithersburg, MD, 1998; NIST Chemistry WebBook http://webbook. nist.gov/chemistry/. (47) Linstrom, P. J.; Mallard W. G. Eds. NIST Chemistry WebBook; NIST Standard Reference Database Number 69; National Institute of Standards and Technology, : Gaithersburg, MD, 2010, http:// webbook.nist.gov. (48) Lsal, M.; Smith, W. R.; Nezbeda, I. Accurate simulation of phase equilibrium for complex fluid mixtures. Application to binaries involving isobutene, methanol, methyl tert-butyl ether, and n-butane. J. Phys. Chem. B 1999, 103, 10496−10505. (49) Lsal, M.; Smith, W. R.; Nezbeda, I. Accurate vapour−liquid equilibrium calculations for complex systems using the Reaction Gibbs Ensemble Monte Carlo simulation method. Fluid Phase Equilib. 2001, 181, 127−146. (50) Smith, J. M.; Van Ness, H. C.; Abbott, M. M. Introduction to Chemical Engineering Thermodynamics, 6th ed.; McGraw-Hill: New York, 2001. (51) Elliott, R. J.; Lira, C. T. Introductory Chemical Engineering Thermodynamics; Prentice-Hall: Upper Saddle River, NJ, 1999; pp 187−189. (52) Ray, J. R. Fundamental treatment of the isobaric-isoenthalpic ensemble. Phys. Rev. A 1986, 34, 2517−2519. (53) Kioupis, L. I.; Maginn, E. J. Pressure-enthalpy driven molecular dynamics for thermodynamic property calculation I. Methodology. Fluid Phase Equilib. 2002, 200, 75−92. (54) Tolman, R. C. The Principles of Statistical Mechanics; Oxford Univ. Press; London, UK, 1938; Section 141. 3270



Article

(55) Hill, T. L. Statistical Mechanics: Principles and Selected Applications; McGraw-Hill: New York, USA, 1956; Chapter 4. (56) Hill, T. L. An Introduction to Statistical Thermodynamics; Addison-Wesley: Reading, USA, 1960; Chapter 2. (57) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids, Vol. I: Fundamentals; Clarendon: Oxford, UK, 1984; Section 3.4.1. (58) Panagiotopoulos, A. Z.; Reid, R. C. On the relationship between pairwise fluctuations and thermodynamic derivatives. J. Chem. Phys. 1986, 85, 4650−4653. (59) Kofke, D. A. Coexistence diagrams of mixtures by molecular simulations. Chem. Eng. Sci. 1994, 49, 2633−2645. (60) Lísal, M.; Budinsky, R.; Vacek, V.; Aim, K. Vapor−liquid equilibria of alternative refrigerants by Molecular Dynamics simulations. Int. J. Thermophys. 1999, 20, 163−174. (61) Stoll, J.; Vrabec, J.; Hasse, H. A set of molecular models for carbon monoxide and halogenated hydrocarbons. J. Chem. Phys. 2003, 119, 11396−11407. (62) Nezbeda, I.; Kolafa, J. The use of control quantities in computer simulation experiments: Application to the exp-6 potential fluid. Mol. Simul. 1995, 14, 153−163. (63) Raabe, G. Molecular simulation studies on the thermophysical properties of the refrigerant blend R-445A. J. Phys. Chem. B 2013, 58, 3470−3476.

3271


Molecular Simulation for Thermodynamic Properties and Process

Recommend Documents