Absorption Refrigeration Cycles with Ammonia-Ionic Liquid Working Pairs Studied by Molecular Simulation.

For absorption refrigeration, it has been shown that ionic liquids have the potential to replace conventional working pairs. Due to the huge number of possibilities, conducting lab experiments to find the optimal ionic liquid is infeasible. Here, we provide a proof-of-principle study of an alternative computational approach. The required thermodynamic properties, i.e., solubility, heat capacity, and heat of absorption, are determined via molecular simulations. These properties are used in a model of the absorption refrigeration cycle to estimate the circulation ratio and the coefficient of performance. We selected two ionic liquids as absorbents: [emim][Tf2N], and [emim][SCN]. As refrigerant NH3 was chosen due to its favorable operating range. The results are compared to the traditional approach in which parameters of a thermodynamic model are fitted to reproduce experimental data. The work shows that simulations can be used to predict the required thermodynamic properties to estimate the performance of absorption refrigeration cycles. However, high-quality force fields are required to accurately predict the cycle performance.

Introduction

Temperature control devices are essential in our society, and they are an integral part of numerous industrial processes.[1,2] The area of application is huge, ranging from electronic devices,[3] food preservation,[4] and the energy management of buildings,[5,6] to industrial cooling and heating.[7] Hence, refrigeration plays a significant role in improving living conditions.[8] In this context, the absorption refrigeration cycle has gained increasing interest, because it enables the efficient utilization of produced low-grade heat, for example, as a byproduct of industrial processes or by solar thermal collectors.[2,6,9−14] The most widely used fluids in absorption cycles are aqueous solutions of LiBr or NH3.[14,15] The H2O/LiBr working pair is often applied in room air conditioning.[15] NH3/H2O is deployed for subzero cooling and free of air infiltration.[14] However, both conventional working pairs are characterized by some inherent drawbacks. The H2O/LiBr system suffers from problems such as corrosion and crystallization, while the NH3/H2O system requires an expensive rectifier to separate the two components.[16−19] Ionic liquids (ILs) have been proposed to overcome some of the problems associated with conventional working fluids.[15,18,20] Room temperature ILs are salts with very low melting points.[18,21] Their main advantages are nonvolatility, chemical and thermal stability, and the possibility of tuning their thermodynamic properties by combining varying anions and cations.[15,22,23] By selecting particular combinations of cations and anions, ILs potentially provide suitable operating temperatures for waste heat applications and might offer wider operating ranges leading to better process flexibility.[15,24] To perform a thermodynamic analysis of the absorption refrigeration cycle, the miscibility, the heat capacity, and the enthalpy of absorption of the working fluids are required. As an initial step in the development of new absorption refrigeration cycles, these thermodynamic properties of various ILs and their mixtures with refrigerants (here NH3) need to be evaluated to find the most promising ones.[25] NH3 is chosen because it is a promising refrigerant due to its low costs, large latent heat,[26] and low impact on the environment (zero ozone depletion and global warming potential).[27] In addition, NH3 is particularly interesting since only a limited amount of experimental measurements exists for NH3/IL systems.[15] For the vast majority of ILs the necessary thermodynamic properties (e.g., densities, vapor pressures, heat capacities, refrigerant solubility) have not been investigated or are experimentally very difficult to access.[22,23,28−30] Moreover, the number of possible anion–cation combinations is huge.[23,28,31] Conducting lab experiments for a large set of ILs is a tedious and, due to the currently high price of ILs, very expensive task. Computer simulations provide a convenient and cheap alternative to predict relevant thermodynamic properties of ILs and to gain a better understanding of the underlying physical behavior.[32−37] Based on the simulations, a preselection of particularly promising ILs for the design of new absorption refrigeration cycles can be made. It is important to note that the used molecular model to describe the NH3/IL mixtures (i.e., the force field) and statistical uncertainties related to the simulation method influence the accuracy of the computer simulations. It is still unknown if the degree of accuracy is sufficient for reasonable performance predictions of the absorption refrigeration cycle. Moreover, molecular simulations of ionic liquids can be time-consuming, which makes the screening of thousands of ILs still impossible. Hence, we provide a proof-of-principle study.

We choose a multiscale modeling approach to evaluate the accuracy of the prediction on the circulation ratio (f) and the coefficient of performance (COP) of the absorption refrigeration cycle using two ionic liquids with NH3 as absorbent. With only minor adjustments, the same approach is applicable to absorption heat pump cycles. On the smallest scale, quantum mechanical (QM) calculations are used to compute the ideal gas heat capacities of ionic liquids. Monte Carlo (MC) simulations are conducted to determine the residual heat capacity, solubilities of the refrigerant in the absorbent, and the enthalpy of absorption. These properties are determined by intermolecular interactions and can be described well by MC simulations.[38] The computed thermodynamic properties are then used in a developed equilibrium-based cycle model of a single-effect absorption refrigeration cycle. Similar computational based approaches have been developed for vapor-compression cycles.[39,40] However, simpler molecular models that have been previously fitted to experimental vapor–liquid equilibrium data were applied in these studies. Moreover, different simulation techniques have been performed.

Here, we want to introduce the concept for absorption refrigeration cycles with NH3/IL working pairs, to explain how the various simulation techniques can be combined to create an initial cycle design, and to assess uncertainties in the performance prediction of the cycle related to the different methods. As exemplary cases, we selected two working pairs: NH3/[emim][Tf2N] and NH3/[emim][SCN]. The IL [emim][Tf2N] is one of the most studied ILs, and a large amount of experimental data is available.[29,41−43] This is crucial to validate our approach. The IL [emim][SCN] has been proposed as a particularly promising candidate for absorption refrigeration cycles with NH3.[44,45] Other promising ILs have been suggested, especially ILs containing functional OH groups. Unfortunately, we could not find any reliable force field for this type of IL. In the traditional approach, experimental measurements are conducted to determine thermodynamic properties of ILs.[25] Afterward, these experimental measurements can be directly used in the cycle design. Frequently, the available experimental measurements do not cover the complete range of conditions necessary in the cycle model. Hence, equations of states (EoSs) or activity coefficient models are fitted to reproduce experimental measurements and are used to extrapolate thermodynamic properties for conditions for which the properties have not been determined directly.[45,46] This extrapolation may introduce undesired errors. Alternative methods to predict thermodynamic properties of mixtures are COSMO-based.[23] However, previous studies suggest that the accuracy of these models is often insufficient for mixtures containing ILs.[47,48] Especially for NH3/IL working pairs the uncertainty can be large.[49] For the screening of large quantities of ILs, less accurate COSMO-based[49] methods seem more suitable due to the advantages in computation time.

This paper is organized as follows. Section presents a brief introduction of the absorption refrigeration cycle along with the assumptions and the equations used in the cycle description. Subsequently, in section we provide the background and the methodology to predict the required thermodynamic properties. Section provides the simulation details for computing the required thermodynamic properties. Section presents the results for the thermodynamic properties and the thermodynamic performance of the cycle. In section , our findings are summarized.

Description of the Absorption Refrigeration Cycle

A schematic diagram of a single-effect absorption refrigeration cycle is shown in Figure a. The basic units of the cycle are a generator (GEN), a condenser (CON), an absorber (ABS), and an evaporator (EVA). To illustrate the changes between the different state points, Figure b shows the relation between vapor pressure (ln(P)) and temperature (−1/T). In the EVA, the pure liquid refrigerant (NH3) is vaporized and provides the cooling load QEVA. The low-pressure refrigerant at state 1 enters the ABS and is exothermically absorbed by the weak solution at state 5 (weak in NH3) coming from the GEN. Thereby, the weak NH3 solution is enriched with NH3 and it becomes the rich solution at state 2. Here, the redundant heat load QABS is released to the surroundings. Before entering the GEN, the rich solution is pumped from state 2 to a higher pressure level and the temperature is increased in the solution heat exchanger (SHX). In the GEN, the heat load QGEN is supplied and the refrigerant vapor is released from the rich solution at state 8. Hence, the concentration of refrigerant in the stream leaving the GEN at state 7 is reduced. The weak NH3 solution is cooled in the SHX. Before entering the absorber, the pressure is lowered in an expansion valve. The refrigerant vapor exiting the GEN at state 8 passes the condenser, where it is condensed while releasing the condensation heat QCON to the surroundings. Subsequently, the pressure of the liquid refrigerant at state 9 is reduced in an expansion valve before the cycle is completed by the liquid refrigerant reentering the EVA. To predict the thermodynamic performance of the single-effect absorption refrigeration cycle with NH3/IL as working pair, the following assumptions are made to facilitate the calculations:

Figure 1

(a) Schematic diagram of a single-effect absorption refrigeration cycle. (b) ln(P)–1/T diagram of the same absorption refrigeration cycle. QGEN, QCON, QABS, and QEVA are respectively the transferred heat at the generator, the condenser, the absorber, and the evaporator. TGEN, TCON, TABS, and TEVA are the corresponding temperatures, and PCON and PEVA the pressures. ṁs and ṁr are the mass flow rates of the strong NH3 solution and of the refrigerant, respectively.

1. The fluid streams are in a steady state.

2. The operating pressures of the EVA and the ABS and, likewise, the operating pressures of the GEN and the CON are equal.

3. The NH3 stream leaving the CON and the EVA is saturated liquid and saturated vapor, respectively.

4. The solution leaving the GEN is at equilibrium.

5. The refrigerant enters the CON as superheated vapor at the end generation temperature.

6. The solution leaving the ABS is subcooled by 5 K.

7. The pinch temperature of the SHX is set to 5 K.

8. Heat and pressure losses are neglected.

9. Throttling is assumed to be an isenthalpic process.

10. The power required for pumping is insignificant, and therefore neglected.

Considering these assumptions, the conditions of all state points can be specified. The temperatures TEVA, TCON, TABS, and TGEN are determined by the heat sources, the surroundings, and the cooling application. The pressures PEVA and PCON are determined by the vapor pressure of pure NH3 assuming that the IL is nonvolatile. Hence, the conditions of pure NH3 in states 1 and 9 are set. The pressures PGEN and PABS follow directly from the assumption of equal pressures at GEN and CON, and EVA and ABS, respectively. Based on the vapor–liquid equilibrium, the mass fractions of each component at the outlet of the GEN and the ABS (subcooled by 5 K) can be determined for TGEN and TABS and the corresponding pressure levels PGEN and PABS. Thereby, states 2 and 7 are specified. Between states 2 and 4, and 7 and 5, the mass fractions of the absorbent do not change (w2 = w3 = w4, w5 = w6 = w7). Thus, all remaining states can be determined from energy balances around the SHX and the isenthalpic throttling process. The cycle conditions in this work are TCON = 35 °C, TABS = 30 °C, TEVA = 10 °C, TGEN = 74–120 °C, PEVA = 6.15 bar, and PCON = 13.5 bar. The enthalpies of pure NH3 are directly obtained from Refprop[50] with the Helmholtz energy EoS.[51] For NH3/IL mixtures, the enthalpy of the solution h̅sol can be calculated asin which h̅NH and h̅IL are the enthalpies of pure NH3 and IL at the specified conditions, wNH is the mass fraction of NH3, and Δh̅abs is the enthalpy of absorption of the mixture. The influence of the pressure on the heat capacity and the enthalpy of the liquid IL c̅IL can be neglected. Hence, the enthalpy of the pure IL can be determined:where h̅o is the reference enthalpy at an arbitrary reference state for eq with To = 250.15 K and Po = 1 MPa. Subsequently, it is straightforward to calculate the exchanged heat with the surroundings (see Figure ):where QGEN, QCON, QABS, and QEVA are respectively the heats, the generator, the condenser, the absorber, and the evaporator exchange with the surroundings. ṁs and ṁr are the mass flow rates of the strong NH3 solution and of the refrigerant, respectively, and h̅ is the specific enthalpy of state i. The mass balance of the absorbent in the solution yieldsw2 and w5 are the mass fractions of the absorbent of the corresponding states in the cycle. Finally, the performance parameters, i.e., the circulation ratio f and the coefficient of performance COP, can be determined as

Predictions of Thermodynamic Properties

Force Field

MC simulations are performed to calculate the solubility of NH3 in [emim][Tf2N] and [emim][SCN], the residual heat capacity of these ILs, and their enthalpies of absorption at various temperatures and pressures. The classical force field developed by Liu et al.[22,30,32,38,52] is used to describe the ILs. It includes intramolecular contributions such as bond stretching, angle bending, and torsions. Intermolecular contributions are described with a Lennard-Jones potential, and electrostatic interactions are considered via the Ewald summation technique with a relative precision of 10–5.[53] The Lennard-Jones interactions are truncated and shifted at 12 Å. No tail corrections are applied. The Lorentz–Berthelot mixing rules are used to calculate the interactions between unlike atoms.[54] The anion and the alkyl part of the cation of the IL molecules are considered flexible, whereas the ring of the cation is rigid. NH3 is described via the TraPPE force field.[55] The TraPPE force field was developed to reproduce the vapor–liquid equilibrium for pure ammonia and an ammonia–methane mixture.[55] These equilibria can be reproduced very well. Besides, the TraPPE force field has the potential to model hydrogen bonds,[55] which is important for [emim]/NH3 mixtures.[32] The force field parameters are not adjusted and were directly taken from Liu et al.,[28] Tenney et al.,[30] and Zhang and Siepmann.[55] The force field parameters of the ILs are summarized in the Supporting Information.

Calculation of the NH3 Solubility in ILs

The osmotic ensemble is well-suited to calculating the solubilities of gases in nonvolatile liquids such as ILs.[30,56] In this ensemble, the temperature, the pressure, the fugacity of the solute, and the number of solvent molecules are fixed. The fugacity and the hydrostatic pressure are coupled via the Peng–Robinson EOS applied for NH3. The volume of the system and the number of solute molecules fluctuate. In equilibrium, the number of solute molecules in the system determines the solubility. A representation of a simulated system is shown in Figure . An inherent difficulty of a fluctuating number of solute molecules (NH3) is the insertion of additional molecules into already dense systems. To overcome this issue, Shi and Maginn[57] proposed the continuous fractional component MC method (CFCMC). In this method, the system consists of whole molecules and a fractional one. Insertion and deletion MC moves of the solute are achieved by scaling the intermolecular interactions (Lennard-Jones and Coulombic) of the fractional molecule. The scaling parameter λ can range between 0 and 1. λ = 0 signifies that the fractional molecule does not interact with the surroundings, whereas λ = 1 implies that the fractional molecule has full interactions with the surrounding molecules. In the MC algorithm, besides trial moves for thermalization, additional trial moves to change the value of λ are included. If λ is changed to a value larger than 1, the molecule is completely inserted and a new fractional molecule is added to the system with λ = λ – 1. If λ is changed to a value smaller than 0, the molecule is removed and a random molecule is chosen as the new fractional one with λ = 1 + λ. This method significantly facilitates the insertion of molecules in dense systems.[58] Subsequently, the fractional molecule and its surroundings can slowly adjust while λ is changed and thus the interactions are gradually switched on. Thereby, the probability of successful insertions increases. To ensure an uniform sampling of the scaling parameter λ, the Wang–Landau scheme[59] is applied during equilibration to determine a biasing function W(λ). The advantage of the Wang–Landau scheme is that it does not require any former knowledge of the biasing function. For more details on computing solubilities with the CFCMC method, the reader is referred to Shi and Maginn[32,57,60] and Ramdin et al.[56,61−63]

Figure 2

Representation of the simulated system consisting of NH3, [emim]+, and [Tf2N]−. Exemplary, the molecules are marked by dashed lines. White, red, gray, purple, yellow, and green spheres represent hydrogen, oxygen, carbon, nitrogen, sulfur, and fluorine atoms, respectively.

In this paper, we compare solubilities computed from MC simulations with the NRTL model which was fitted to reproduce the experimental data of Yokozeki and Shiflett.[29,44] In the range of experimental data, the NRTL model describes the experiments very well. Consequently, the calculated solubilities are compared to the NRTL model at the exact cycle conditions. Details of the NRTL model and the fitted parameters can be taken from Wang and Infante Ferreira.[45]

Calculation of the IL Heat Capacity

The heat capacity describes how much the temperature changes if a certain amount of heat is added. The heat capacity at constant pressure C is defined aswhere T, P, and ⟨H⟩ are the temperature, the pressure, and the average enthalpy computed in the NPT ensemble, respectively. The enthalpy is defined asUintra and Uinter are the intramolecular and the intermolecular contributions to the potential energy, respectively, K is the kinetic energy, and V is the volume of the system. The average enthalpy can be expressed as the sum of two separate contributions, i.e., the ideal gas enthalpy and the residual enthalpy:[52]where ⟨Hres⟩ is the ensemble average of the residual enthalpy, N is the number of molecules the system comprises, and kB is the Boltzmann constant. The ensemble average of the ideal gas enthalpy ⟨Hig⟩ is defined as the sum of the intramolecular energy Uintra, the kinetic energy K, and the PV term which can be replaced with NkBT according to the ideal gas law. Consequently, the total heat capacity can be split into an ideal gas and a residual part by differentiating the two enthalpy contributions individually:

It is assumed that the ensemble average of the ideal gas enthalpy is independent of the residual enthalpy.[52] The splitting of the heat capacity is convenient, because intramolecular potentials used in classical force field based molecular simulations may result in large deviations for the ideal part of the heat capacity.[52] The causes of these deviations are the overestimate of the vibrational energy due to the harmonic oscillator approximation and the negligence of off-diagonal coupling terms.[52] Likewise, it is unfavorable to use QM calculations to predict the residual contribution of the heat capacity due to rapidly increasing computational costs for systems containing more than a single molecule or ion. Therefore, a separate calculation of the ideal and residual terms is advantageous. The derivative of the residual enthalpy with respect to temperature can be determined from fluctuations of thermodynamic variables throughout the simulations. Here, the procedure of Lagache et al.[64] is applied:

The configurational enthalpydiffers from the enthalpy H by not including the kinetic energy K. The ideal contribution to the heat capacity can be reliably predicted from ab initio electronic structure calculations.[52] Therefore, QM calculations of isolated ions are performed where only intramolecular interactions are considered. The combination of both computational methods has been successfully applied to predict the heat capacity of ILs and other molecules.[22,28,30,38,52,64−66]

The results of this simulation based approach are compared to experimental data of Paulechka et al.[67] and Navarro et al.[68] for [emim][Tf2N] and [emim][SCN], respectively. These experimental results are fitted to a polynomial function to describe the temperature dependency of CIL.

Calculation of the Enthalpy of Absorption

The enthalpy of absorption Δh̅abs can be calculated from MC simulations following the procedure of Shi and Maginn.[32]Equation can be rearranged toThree separate MC simulations at the same conditions are necessary to compute Δh̅abs: one for the enthalpies of the solution (h̅sol), one for NH3 (h̅NH), and one for the IL (h̅IL). Here, the enthalpies of pure NH3 and the ILs are computed from simulations in the NPT ensemble whereas the enthalpy of the solution is computed from simulations in the osmotic ensemble, simultaneously with the solubility of NH3. The specific enthalpies in both ensembles are calculated via h̅ = u̅ + Pv̅, where u̅ is the specific energy and v̅ is the specific volume. The enthalpies of pure NH3 and the ILs are computed at T and P corresponding to the condition of the mixture. Note that, at the conditions chosen for the cycle design, pure NH3 is in the gaseous state.

There is no experimental data for the enthalpy of absorption available for NH3/IL mixtures. Therefore, the mixing enthalpy Δh̅mix derived from EoS or activity coefficient models has to be used in the traditional cycle design. To compare the enthalpy of absorption computed from MC simulations with the mixing enthalpy obtained from EoSs or activity coefficient models with the liquid state as reference, the latent heat Δh̅lat of pure NH3 at the same T and P has to be subtracted:The latent heat of NH3 is taken from Refprop.[50,51] Here, the calculation results are compared with the predictions of the Redlich–Kwong EoS (RK-EoS) and mixing rules based on experimental vapor–liquid equilibrium data usingwhere h̅NHres, h̅ILres, and h̅solres are the residual enthalpies for NH3, the IL, and the solution at liquid state, respectively. The detailed procedure and the necessary critical parameters and binary interaction parameters are explained in the work of Yokozeki and Shiflett.[29,44] In a publication by Wang et al.,[69] it was shown that the RK-EoS performed well in comparison to other EoS models, activity coefficient models, and the Clausius–Clapeyron equation. It is possible to predict mixing enthalpies from the NRTL model. However, the NRTL model is not recommended for determining mixing enthalpies due to a purely empirical temperature dependency of the fitting parameters.[44,69]

Simulation Details

The simulations are performed with the RASPA software package.[58,70] In RASPA, the number of MC steps conducted during every cycle is equal to the total number of molecules considered in the simulation with a minimum of 20 MC steps per cycle. The total number of simulation cycles is divided into five blocks. The provided statistical uncertainty in the computed properties is the standard error calculated from the standard deviation of the block averages.

The solubility of NH3 in [emim][Tf2N] and [emim][SCN] is computed by conducting MC simulations in the osmotic ensemble. NH3 molecules are inserted or removed via the CFCMC technique.[57,58] A production run of 5 million cycles is performed. Simulations are performed with 52 and 70 IL molecules for [emim][Tf2N] and [emim][SCN], respectively. Three different MC trial moves with an equal probability are considered for the ILs: translational, rotational, and configurational bias to account for different configurations of the molecules. For NH3, translational, rotational, and λ moves are considered with an equal probability. The probability of a volume change move is 1%.

The calculation of the heat capacity is divided in two parts. QM calculations with the Gaussian software package[71] are performed to determine the ideal gas part of the heat capacity. In these calculations the isolated ions [emim], [Tf2N], and [SCN] are considered. Possible conformers are analyzed using molecular mechanics with the Merck molecular force field,[72] and the results are compared to the literature.[73] Paulechka et al.[73] concluded that the B3LYP functional[74−77] with a 6-31+G(2df,p) basis set is adequate for calculation of the frequencies of [emim] and [Tf2N] which is also chosen here. First, the geometry is optimized, and subsequently, a frequency analysis is performed. A scaling factor of 0.965 is applied to scale the calculated vibrational frequencies consistent with the NIST database.[78] In the NIST database,[78] these kind of well-established scaling factors are tabulated to improve the agreement between QM calculations and experiments for different basis sets and theories. To compute the residual part of the heat capacity of the pure ILs, MC simulations in the NPT ensemble are conducted. For [emim][Tf2N] and [emim][SCN], simulations are performed with 55 and 70 molecules, respectively. After reaching equilibrium, between 22 and 24 million MC cycles are conducted to compute the residual heat capacity according to eq . Again, translational, rotational, and configurational-bias MC trial moves with equal probability are considered for the ILs and the probability of a volume change move is 1%.

For the enthalpy of absorption, the energies for NH3/IL mixtures are computed in the osmotic ensemble simultaneously with the NH3 solubility. For pure ILs, energies are computed from simulations in the NPT ensemble together with the residual heat capacity. Additional, NPT simulations at the same conditions are conducted for pure NH3. These simulations are performed with 160 molecules. First the systems are equilibrated, and then continued for 105 production cycles. Translational and rotational MC trial moves with an equal probability are considered for NH3. Volume change moves have a probability of 1%.

Results and Discussion

NH3 Solubility in ILs

In accordance with the conditions occurring in the absorption refrigeration cycle, the solubility of NH3 is computed in [emim][Tf2N] and [emim][SCN] at 308.15, 347.15, 373.15, and 393.15 K for pressures between 4 and 19 bar. The computed results for [emim][Tf2N] and [emim][SCN] are shown in Figure , parts a and b, respectively. The simulation results (colored) are compared to the NRTL model (black). Results of experimental measurements are shown exemplary as open symbols. The predictions of the MC simulations are in qualitative agreement with the NRTL model. The comparison between experimental and simulation results shows that the high solubility of NH3 in the two ILs and the influence of temperature and pressure are reproduced. The average relative deviation between the NRTL model and MC simulations is 17–28% for [emim][Tf2N], and 20–29% for [emim][SCN], whereby our MC simulations predict a higher solubility of NH3 in the ILs than the NRTL model. The relative deviations decrease at higher pressures, indicating that the force field performs better at higher loadings of NH3. At higher loadings of NH3, NH3–NH3 interactions become more important. The TraPPE force field of NH3 is designed for vapor–liquid equilibrium data and describes these interactions well. The magnitude of the deviations observed between computed values and experimental results is comparable with the work of Shi and Maginn[32] and Urukova et al.,[79] who also computationally investigated highly soluble gases in ILs. Shi and Maginn[32] also computed the absorption of NH3 in [emim][Tf2N]. For comparison, one of the absorption isotherms computed by these authors is shown in Figure a (dashed line). The force field of [emim][Tf2N] used in this work is almost identical to the one used by Shi and Maginn.[32] The only difference is that Shi and Maginn[32] consider the IL completely flexible. However, their force field describes NH3 with smaller charges and a smaller Lennard-Jones energy parameter. These differences in the description of NH3 result in a poorer prediction of the vapor–liquid equilibrium of pure NH3. In contrast to our results, these authors underpredict the solubility of NH3 in [emim][Tf2N]. This suggests that NH3/IL interactions play a major role. Deviations in the predicted solubility can be ascribed to the applied force fields. Hence, improvements in the force field parameters are required. Nevertheless, MC simulations can be used to predict thermodynamic properties at high temperatures and pressures which are otherwise difficult to determine experimentally.[40] A recently published COSMO-RS study reports an average underprediction of 34% for vapor pressures of a large set of NH3/IL mixtures. However, detailed conditions are not provided, which complicates a direct comparison. Experimental solubilities of NH3 in [emim][Tf2N] have been measured up to 347.6 K. Using the NRTL model to extrapolate the solubility results in a change of the curvature of the absorption isotherm for 393.15 K (see Figure a). This change in curvature is not expected since the nature of the physical interactions between the IL and NH3 should remain similar. Therefore, the change in curvature might be an artifact of the extrapolation using the NRTL model. Hence, we recommend caution when using the NRTL model to extrapolate experimental data. To extrapolate experimental data, the qualitative behavior predicted by MC simulations may be more reliable.[40]

Figure 3

Computed NH3 solubilities (blue/green/cyan/magenta) in (a) [emim][Tf2N] and (b) [emim][SCN], compared to solubilities calculated with the NRTL model (black), experimental data (○), and simulation results of Shi and Maginn[32] (red) at 308.15 (▼), 347.15 (●), 373.15 (■), and 393.15 K (⧫). The determined standard error is smaller than the size of the symbols.

IL Heat Capacity

The total heat capacity for [emim][Tf2N] and [emim][SCN] is obtained by adding the ideal gas part and the residual part (eq ). In Figure the computed and experimental heat capacities[41,42,67,68,80] are compared as functions of temperature. The computed values for both ILs are in agreement with experimental data. Average deviations between the experimental and the computed heat capacities are around 4 and 2.5% for [emim][Tf2N] and [emim][SCN], respectively. For temperatures from 303.15 to 333.15 K, the residual heat capacity is computed at 6.1505 bar (evaporation pressure), while for temperatures from 343.15 to 393.15 K, it is computed at 13.508 bar (condensation pressure). The experimental uncertainty is large, which is also depicted by the high degree of scatter between the different experimental data sets. For usage in the cycle model, the computed heat capacities are fitted to a quadratic polynomial in temperature. The resulting parameters are shown in Table . Subsequently, the enthalpies of the pure IL at different cycle conditions can be calculated from eq .

Figure 4

Comparison between computed total heat capacities (blue), computational results of Tenney et al.[30] (black), and experimental measurements of Ge et al.[42] (green), Paulechka et al.[67] (cyan), Ferreira et al.[41] (magenta), Navarro et al.[68] (purple), and Ficke et al.[80] (orange) for (a) [emim][Tf2N] and (b) [emim][SCN].

Table 1

Parameters Used To Fit the Polynomial Describing the Temperature Dependency of the Heat Capacity (c̅IL = a + bT + cT2) to Our Simulation Results

IL	a [J kg–1 K–1]	b [J kg–1 K–2]	c [J kg–1 K–3]
[emim][Tf2N]	–429.51	7.338	–6.429 × 10–3
[emim][SCN]	–552.20	10.947	–1.167 × 10–2

Enthalpy of Absorption

The computed enthalpies of absorption from MC simulations are reported in Table and compared to the results obtained from the RK-EoS.

Table 2

Enthalpies of Absorption at Different Cycle Conditions Computed from MC Simulations and the RK-EoS (Eqs and 22)

			MC simulations		RK-EoS
IL	T [K]	P [bar]	wNH3 [kg kg–1]	Δh̅abs [kJ kg–1]	wNH3 [kg kg–1]	Δh̅abs [kJ kg–1]
	308.15	6.15	0.1086	–144.2 ± 9	0.0688	–63.27
[emim]	347.15	13.15	0.0766	–94.6 ± 12	0.0576	–42.87
[Tf2N]	373.15	13.15	0.0404	–51.5 ± 17	0.0290	–15.98
	393.15	13.15	0.0277	–28.7 ± 7	0.0123	–4.61
	308.15	6.15	0.2235	–293.8 ± 3	0.1159	–94.42
[emim]	347.15	13.15	0.1480	–184.9 ± 5	0.1045	–58.67
[SCN]	373.15	13.15	0.0822	–99.4 ± 6	0.0580	–0.78
	393.15	13.15	0.0544	–54.7 ± 6	0.0357	22.71

The enthalpies of absorption computed from MC simulations show a consistent trend and are negative for all computed conditions. The determined absolute values increase as the temperature decreases and as the concentration of NH3 increases. This behavior signifies attraction between the ILs and NH3, which is consistent with the results reported by Shi and Maginn.[32] Enthalpies of absorption from MC simulations are consistently larger than the ones from the RK-EoS. In addition, the EoS predicts positive heats of absorption for [emim][SCN] and temperatures higher than 373 K. Without experimental data, it is difficult to comment on the accuracy of the obtained data and further experiments are necessary. However, it is known that it is problematic to describe the phase behavior of polar nonvolatile compounds such as ILs with a cubic EoS.[81] Hence, the enthalpies of absorptions predicted from molecular simulations might be more reliable.

Circulation Ratio

The circulation ratio f is defined as the ratio between the mass flow rate of the strong NH3 solution leaving the absorber and the mass flow rate of refrigerant (see eq ). The value of f depends solely on the solubility of NH3 in the absorbent. It is an important performance parameter as it is directly related to the size and cost of the equipment.[45] The calculated values for f as a function of TGEN following from MC simulations and the NRTL model for both ILs are compared in Figure . Deviations between the circulation ratio predicted from the NRTL model and from MC simulations can be observed. The average deviation for NH3/[emim][Tf2N] is 50%, and it is 67% for NH3/[emim][SCN]. The reason for these deviations is the discrepancy in the calculated NH3 solubility (see eq ). This deviation is particularly important for low-end generation temperatures. At these temperatures, the solubility of NH3 in the ILs is very high (above 60 mol %). Hence, the mass flow rate of the IL in comparison to NH3 is relatively low. This results in a sensitivity of f toward changes in the solubility of NH3. For both ILs, the simulations overestimate the solubility of NH3 and therefore predict lower mass fractions of ILs. As TGEN increases, the mass fraction of the ILs in the weak solution increases and f decreases to a nearly constant value at high-end generation temperatures. A high circulation ratio raises the generation heat input according to eq . Therefore, it is not recommended to operate an absorption refrigeration cycle at low TGEN.

Figure 5

Comparison between f values calculated with NH3 solubilities from MC simulations (colors) and from the NRTL model (black) for [emim][Tf2N] (●) and [emim][SCN] (■). The cycle conditions in this work are TCON = 35 °C, TABS = 30 °C, TEVA = 10 °C, TGEN = 74–120 °C, PEVA = 6.15 bar, and PCON = 13.5 bar.

Coefficient of Performance

Figure shows the comparison between the COP calculated from simulations (MC simulations and QM calculations) and the traditional approach (NRTL/EoS model) for both working pairs as a function of TGEN. A different trend of the COP can be observed for both approaches. The COP predicted from simulations is almost constant over the whole range of considered end generation temperatures. In contrast, the COP based on the NRTL/EoS model increases rapidly with TGEN for low temperatures. This behavior results from a strong temperature dependency of the circulation ratio for low TGEN (see Figure ). For our simulations, a smaller value of f than for the NRTL model is predicted for both ILs. A smaller f decreases the required heat input for the generator (see eq ) and therefore has a favorable impact on the COP, which explains the larger COP for the simulation approach for low TGEN. However, for high TGEN, the larger enthalpy difference predicted from simulations between the solutions entering (state 4) and leaving (state 7) the generator affects the COP even more (see Table ) and results in smaller COPs for the computational approach. The difference in the heat capacity predicted from simulations and measured experimentally is rather small and is therefore not causing significant changes in the COP prediction. The average deviations between the two approaches for the COP are 32 and 38%, for NH3/[emim][Tf2N] and NH3/[emim][SCN], respectively. Overall, the results show that, for the investigated cycle model, NH3/[emim][Tf2N] performs better than NH3/[emim][SCN] for a preset TGEN. This statement is true for both the NRTL/EoS model and the simulations.

Figure 6

Comparison between the COP values calculated from simulations (colors) and from the NRTL/EoS model (black) for [emim][Tf2N] (●) and [emim][SCN] (■). The cycle conditions in this work are TCON = 35 °C, TABS = 30 °C, TEVA = 10 °C, TGEN = 74–120 °C, PEVA = 6.15 bar, and PCON = 13.5 bar.

Conclusions

It has been shown that the computational prediction of the thermodynamic properties from MC simulations along with QM calculations can be used to predict the performance of a single-effect absorption refrigeration cycle. The approach has been demonstrated for NH3/[emim][Tf2N] and NH3/[emim][SCN]. For these working pairs solubility, heat capacity, and enthalpy of absorption are computed at various cycle conditions. Subsequently, the circulation ratio f and the coefficient of performance COP are determined with a developed model of the cycle. The only prerequisite for the simulation based approach is a force field describing the interactions between the absorbent and the ions of the ionic liquids, as well as the intramolecular interactions of the ions. The determined thermodynamic properties and performance parameters are compared to a NRTL/EoS model which is fitted to experimental data. Average deviations between the traditional and the simulation based approach for the COP of 32 and 38%, for NH3/[emim][Tf2N] and NH3/[emim][SCN], respectively, were observed. For the circulation ratio f, the average deviations between the traditional and the simulation based approach are around 50% for NH3/[emim][Tf2N] and 67% for NH3/[emim][SCN]. The observed deviations are mainly due to discrepancies of the enthalpies of absorption predicted with the RK-EoS and with simulations, and the overprediction of the solubilities of NH3 in the ILs by MC simulations. These results show that accurate force fields for the investigated working pairs are crucial. To enable the systematic computational screening of working pairs for absorption refrigeration cycles, this issue needs to be further addressed. Two major problems are worth mentioning concerning the traditional approach. First, the extrapolation with the NRTL model of the NH3 solubility in the investigated ionic liquids to temperatures higher than experimentally measured is questionable. Deviations in the solubility can significantly influence f and thereby the COP. Second, the absence of experimental data for the enthalpy of absorption introduces a considerable uncertainty for the predicted enthalpy of absorption. Therefore, an error of unknown magnitude exists for this thermodynamic property. This error can also have a significant impact on the COP. Simulation techniques show potential to mitigate these issues. They can play an important role in the prediction of thermodynamic properties for conditions under which experiments are challenging to perform (such as high temperatures and high pressures) and the prediction of mixture properties (such as the enthalpy of absorption). We feel that the presented computational approach seems to be the best choice in the complete absence of experimental data.