Thermal Energy Storage of R1234yf, R1234ze, R134a and R32/MOF-74 Nanofluids: A Molecular Simulation Study

Thermal energy storage can be carried out by working fluid adsorbing and desorbing in porous materials. In this paper, the energy storage properties of four refrigerants, R1234yf, R1234ze, R134a and R32, with M-metal organic framework (MOF)-74 (M = Zn, Ni, Mg, Co) nanoparticles are investigated using molecular dynamics simulations and grand canonical Monte Carlo simulations. The results show that M-MOF-74 can adsorb more R32 and R134a than R1234yf and R1234ze, as the molecular structures of R32 and R134a are smaller than those of R1234yf and R1234ze. Mg-MOF-74 owns a higher adsorbability than the other MOFs. The energy storage properties of the studied refrigerants can be enhanced when the sum of thermodynamic energy change of MOF particles and the desorption heat of fluid in MOFs is larger than the enthalpy change of pure organic fluid. The R1234yf/M-MOF-74 (M = Co, Mg, Ni) nanofluid can store more energy than other refrigerants/M-MOF-74 (M = Co, Mg, Ni) nanofluid. The energy storage enhancement ratios of R1234yf, R1234ze and R134a with Mg-MOF-74 nanoparticles are higher than those of other M-MOF-74 (M = Co, Ni, Zn) materials.


Introduction
The energy crisis and environmental pollution are serious issues in human society, as a consequence of the devastating consumption of fossil fuels such as coal, oil and natural gas. In China, the waste heat generated in industries is equivalent to the combustion energy of 340 million tons of standard coal each year. Finding an efficient way to recover the low-grade energy, such as industrial waste heat, geothermal and solar energy, is, hence, essential to combat these issues [1].
The organic Rankin cycle (ORC) uses an organic refrigerant as a working fluid to recover the heat from low-grade energy resources and turns it into electricity, providing a potential way to use the low-grade energy [2]. Further, the efficiency of ORC in low-grade energy recovery could be remarkably enhanced by increasing the thermal properties of the working fluid. McGrail et al. [3] observed reversible adsorption of R123 molecules in metal organic framework (MOF) nanoparticles. They showed that adding MOF nanoparticles into refrigerant (also known as metal-organic heat carrier nanofluids (MOHCs)) could significantly improve the heat capacity of the working fluid and subsequently enhance the efficiency of ORC in recovering the low-grade energy. During the evaporation process in ORC, excess heat is absorbed by MOHCs from the heat resource for the desorption of adsorbate molecules from MOF particles. Further, excess heat is released by the MOHCs during the condensation process as a consequence of the exothermic adsorption of refrigerant molecules into MOF particles. Additionally, the thermal conductivity of the working fluid is enhanced by adding nanoparticles [4]. The reported research works [5,6] indicated that the enhancement of the thermal properties of nanofluid are mainly due to the effective medium theory. Besides, the Brownian motion, layering phenomenon, clustering, ballistic phonon motion, thermal boundary resistance and mass difference scattering will also impact the thermal properties of nanofluid [7,8].
MOFs are highly crystallized by metal ions or metal clusters and organic ligands through the coordination bond or intermolecular interactions, possessing good thermal stability, high specific surface area and strong adsorption affinity for adsorbates. The adsorptive performance of MOFs is found superior to the traditional nanoporous materials such as zeolites and activated carbons [9], making them a promising material for energy storage and gas separation [10,11]. Among a wide range of MOF structures, MOF-74 owns unique secondary building units that have a strong fault tolerance of different metal particles [12]. Nevertheless, its adsorptive performance for organic refrigerates is little understood [13]. Thus, M-MOF-74 (M = Zn, Ni, Mg, Co) is employed as the adsorbent added in the organic refrigerants in this work.
Hitherto, the development of refrigerants has evolved for four generations, with the third and fourth generation of refrigerants being most widely used in practice. In resemblance to the thermodynamic properties of R410a, R32 (difluoromethane, also known as HFC-32), a halogenated hydrocarbon, is treated as an ideal substitute for the second generation ozone-depleting substances [14]. In addition, R134a (tetrafluoroethane, also known as HFC-134a) with zero ODP (ozone depletion potential) is widely used in industry for its excellent cooling performance [15]. Unlike the third generation refrigerants, R32 and R134a, the fourth generation refrigerants are sustainable and environmental friendly, which are proposed to reduce the emission of fluorinated greenhouse gases [16]. Among the hydrofluoroolefins (HFOs) refrigerants, R1234ze (1,3,3,3-tetrafluoropropene, also known as HFO-1234ze) and R1234yf (2,3,3,3-tetrafluoropropene, also known as HFO-1234yf) have demonstrated to be the best working substances to replace the third generation refrigerant, R134a, for their zero ODP, low GWP (global warming potential) and short residence time in the atmosphere [17]. Meanwhile, it has been widely recognized that both R1234yf and R1234ze are promising working fluids in the refrigeration cycle, heat pump and ORC [18,19].
For comparative purposes, the adsorption and desorption heats of the third and fourth generation refrigerant molecules, R32, R134a, R1234yf and R1234ze, in MOF-74 were measured by using grand canonical Monte Carlo (GCMC) simulations. In this connection, the maximized heat capacity among the four different nanofluidic fluids comprising the organic refrigerant (R32, R1234yf, R1234ze and R134a) and the MOF-74 particles is unveiled with the aid of molecular dynamics (MD) simulations [20][21][22].

Materials and Methods
The thermal energy stored in MOHCs can be calculated following [3,23], or be written as, where the enthalpy gain of MOHCs (∆h MOHCs ) consists of the enthalpy change of pure organic fluid (∆h Fluid ), the energy change of the MOF nanoparticles (( T 1 T 0 C p dT) MOFs ) and the enthalpy of desorption (∆h desorption ) of the refrigerant molecules for the temperature difference of T 1 − T 0 . Here, T 1 and T 0 are the temperatures of the heating and cooling resources. C p is the heat capacity of MOF particles, and x is the mass fraction of MOF particles in MOHCs. It can be concluded that the MOHCs are able to store more thermal energy than the pure fluid when the sum of thermodynamic energy change of MOF particles and desorption heat of fluid in MOFs are larger than the enthalpy change of pure organic fluid. In other words, as long as the second term in the right side of Equation (2) is positive, the energy storage in MOHCs is enhanced compared to that in the pure organic fluid.
Since the enthalpy of pure species, h Fluid (R32, R1234yf, R1234ze and R134a), is available in NIST [24], the first term on the right side of Equation (2), ∆h Fluid , can be determined as ∆h Fluid = h Fluid.T − h Fluid.T0 , where T 0 and T are the reference and interested temperatures of the working fluid. Meanwhile, the heat capacity of MOF, C p , was measured in our equilibrium molecular dynamic (EMD) simulations. The last part, ∆h desorption , was calculated as the isosteric heat of adsorption in our GCMC simulations. Subsequently, ∆h MOHCs can be readily calculated.
The EMD and GCMC simulations were performed by using Materials Studio [25]. The COMPASS force field [26] was adopted to describe the intra-and inter-molecular interactions. The Ewald summation method [27] was employed to correct the long-range Coulomb interactions.

MD Computational Details
M-MOF-74 is formed by divalent metal ions (Zn, Mg, Co, Ni, etc.) connected with ligand 2, 5-dihydroxy terephthalic acid. It forms a 2-dimensional six-party channel and a 3-dimensional space network structure like a honeycomb, and metal atoms are octahedrally coordinated (dominated by five oxygen atoms and a water molecule) [28]. The molecular structure of MOF-74 (M = Ni) employed in this work is depicted in Figure 1. The MOF-74 sample comprised of 1 × 1 × 4 unit cells, with 648 atoms (including 288 carbon atoms, 72 oxygen atoms, 216 hydrogen atoms and 72 Ni atoms) was used in our GCMC and EMD simulations. The lattice parameters for the MOF-74 structure are a = b = 25.7856 Å, c = 27.0804 Å, α = β = 90 • and γ = 120 • .
The molecular configurations of refrigerants R1234yf, R1234ze, R134a and R32 are presented in Figure 2. In addition, the distribution of the adsorption sites of R1234ze in Co-MOF-74 is shown in Figure 3, which was obtained at 353 K and 5 MPa. store more thermal energy than the pure fluid when the sum of thermodynamic energy change of MOF particles and desorption heat of fluid in MOFs are larger than the enthalpy change of pure organic fluid. In other words, as long as the second term in the right side of Equation (2) is positive, the energy storage in MOHCs is enhanced compared to that in the pure organic fluid.
Since the enthalpy of pure species, hFluid (R32, R1234yf, R1234ze and R134a), is available in NIST [24], the first term on the right side of Equation (2), ∆hFluid, can be determined as ∆hFluid = hFluid.T − hFluid.T0, where T0 and T are the reference and interested temperatures of the working fluid. Meanwhile, the heat capacity of MOF, Cp, was measured in our equilibrium molecular dynamic (EMD) simulations. The last part, ∆hdesorption, was calculated as the isosteric heat of adsorption in our GCMC simulations. Subsequently, ∆hMOHCs can be readily calculated.
The EMD and GCMC simulations were performed by using Materials Studio [25]. The COMPASS force field [26] was adopted to describe the intra-and inter-molecular interactions. The Ewald summation method [27] was employed to correct the long-range Coulomb interactions.

MD Computational Details
M-MOF-74 is formed by divalent metal ions (Zn, Mg, Co, Ni, etc.) connected with ligand 2, 5dihydroxy terephthalic acid. It forms a 2-dimensional six-party channel and a 3-dimensional space network structure like a honeycomb, and metal atoms are octahedrally coordinated (dominated by five oxygen atoms and a water molecule) [28]. The molecular structure of MOF-74 (M = Ni) employed in this work is depicted in Figure 1. The MOF-74 sample comprised of 1 × 1 × 4 unit cells, with 648 atoms (including 288 carbon atoms, 72 oxygen atoms, 216 hydrogen atoms and 72 Ni atoms) was used in our GCMC and EMD simulations. The lattice parameters for the MOF-74 structure are a = b = 25.7856 Å , c = 27.0804 Å , α = β = 90° and γ = 120°.
The molecular configurations of refrigerants R1234yf, R1234ze, R134a and R32 are presented in Figure 2. In addition, the distribution of the adsorption sites of R1234ze in Co-MOF-74 is shown in Figure 3, which was obtained at 353 K and 5 MPa.      In EMD simulations, the thermodynamic energies of M-MOF-74 nanoparticles at different temperatures (293 K, 313 K, 333 K, 353 K, 373 K and 393 K) were computed in canonical (NVT) ensemble via the Forcite module in Materials Studio. The time step was set as 1 fs, with a 200-ps run being used to equilibrate the MOF structures. A 500-ps run after equilibration was further conducted to obtain the thermodynamic energy. The assignment of charges in MOF structures is based on the COMPASS force field, keeping the structures electronically neutral. The Berendsen thermostat was employed to adjust the temperature of the system. The cut-off distance is 12.5 Å . Periodic boundary conditions were applied in the X, Y and Z directions in all of our GCMC and EMD simulations.

GCMC Computational Details
The GCMC simulations are able to calculate the adsorption isotherms (293 K, 313 K, 333 K, 353 K, 373 K and 393 K) of refrigerants (R1234yf, R1234ze, R134a and R32) in the M-MOF-74 nanoparticles via the Sorption module in Materials Studio. The pressure range of 1-10,000 kPa is investigated in this work. The fugacity is calculated by the Peng-Robinson equation. For each point of the adsorption isotherm, the equilibration time was 10,000 cycles with another 100,000 cycles for statistics. Here, the simulation runs 5,500,000 cycles to obtain the adsorption isotherm at each temperature (293 K, 313 K, 333 K, 353 K, 373 K and 393 K).

Thermodynamic Energy of MOF-74
Our previous work [29] indicated that the present method is an efficient approach to obtain the thermodynamic parameters of MOFs. The changes of the thermodynamic energies of M-MOF-74 structures at different temperatures relative to that at the reference temperature of 293 K are plotted in . These values are larger than ρCp-MOF-5 = 1.24 (J/cm 3 ·K), which was calculated in our previous work [29]. The difference of heat capacities between MOF-74 and other MOFs was expected to be a result of the difference in their structures and components [3,30]. In EMD simulations, the thermodynamic energies of M-MOF-74 nanoparticles at different temperatures (293 K, 313 K, 333 K, 353 K, 373 K and 393 K) were computed in canonical (NVT) ensemble via the Forcite module in Materials Studio. The time step was set as 1 fs, with a 200-ps run being used to equilibrate the MOF structures. A 500-ps run after equilibration was further conducted to obtain the thermodynamic energy. The assignment of charges in MOF structures is based on the COMPASS force field, keeping the structures electronically neutral. The Berendsen thermostat was employed to adjust the temperature of the system. The cut-off distance is 12.5 Å. Periodic boundary conditions were applied in the X, Y and Z directions in all of our GCMC and EMD simulations.

GCMC Computational Details
The GCMC simulations are able to calculate the adsorption isotherms (293 K, 313 K, 333 K, 353 K, 373 K and 393 K) of refrigerants (R1234yf, R1234ze, R134a and R32) in the M-MOF-74 nanoparticles via the Sorption module in Materials Studio. The pressure range of 1-10,000 kPa is investigated in this work. The fugacity is calculated by the Peng-Robinson equation. For each point of the adsorption isotherm, the equilibration time was 10,000 cycles with another 100,000 cycles for statistics. Here, the simulation runs 5,500,000 cycles to obtain the adsorption isotherm at each temperature (293 K, 313 K, 333 K, 353 K, 373 K and 393 K).

Thermodynamic Energy of MOF-74
Our previous work [29] indicated that the present method is an efficient approach to obtain the thermodynamic parameters of MOFs. The changes of the thermodynamic energies of M-MOF-74 structures at different temperatures relative to that at the reference temperature of 293 K are plotted in Figure 4. The thermodynamic energy of M-MOF-74 increases linearly with temperature. The increment of thermodynamic energy per unit temperature is the heat capacity of MOFs, C p , i.e., ρC p-Co-MOF-74 = 3.34 (J/cm 3 ·K), ρC p-Mg-MOF-74 = 3.74 (J/cm 3 ·K), ρC p-Ni-MOF-74 = 2.14 (J/cm 3 ·K), ρC p-Zn-MOF-74 = 2.62 (J/cm 3 ·K). These values are larger than ρC p-MOF-5 = 1.24 (J/cm 3 ·K), which was calculated in our previous work [29]. The difference of heat capacities between MOF-74 and other MOFs was expected to be a result of the difference in their structures and components [3,30].   Figure 5 depicts the adsorption isotherms of refrigerants including R1234yf, R1234ze, R134a and R32 in M-MOF-74 at different temperatures, in which symbols represent the results measured in our GCMC simulations. Among the four refrigerants investigated, R32 was demonstrated to have the highest adsorption capacity in the M-MOF-74 within the temperature range considered, evident in Figure 5c,g,k,o, followed by R134a. This was due to the molecule structure of refrigerants. As is shown in Figure 2, the R32 molecule has a structure with a smaller effective diameter. R1234yf and R1234ze own structures with large effective diameters. At the high pressures, when the adsorption approached saturation, the adsorbate molecules experienced a strong repulsive force from the tightly-packed neighbouring molecules and the host structure [31,32]. This repulsive interaction should be further enhanced for the molecules having a larger size, leading to reduction of the adsorption volume available in MOF-74 for R1234yf, R1234ze and R134a compared to R32. Consequently, the highest adsorption capacity was achieved for R32 in the pores of MOF-74. On the other hand, the Mg-MOF-74 adsorbed more refrigerants than other MOF-74s. This is because the ionic radius of Mg 2+ is smaller than those of Co 2+ , Ni 2+ and Zn 2+ , which results in a larger adsorption volume in the Mg-MOF-74. More importantly, the interactions between Mg 2+ and adsorbate molecules were stronger than the counterpart interactions for Co 2+ , Ni 2+ and Zn 2+ ions [33]. The enthalpy of desorption calculated by GCMC simulations is shown in Appendix A, confirming this statement.    Figure 5c,g,k,o, followed by R134a. This was due to the molecule structure of refrigerants. As is shown in Figure 2, the R32 molecule has a structure with a smaller effective diameter. R1234yf and R1234ze own structures with large effective diameters. At the high pressures, when the adsorption approached saturation, the adsorbate molecules experienced a strong repulsive force from the tightly-packed neighbouring molecules and the host structure [31,32]. This repulsive interaction should be further enhanced for the molecules having a larger size, leading to reduction of the adsorption volume available in MOF-74 for R1234yf, R1234ze and R134a compared to R32. Consequently, the highest adsorption capacity was achieved for R32 in the pores of MOF-74. On the other hand, the Mg-MOF-74 adsorbed more refrigerants than other MOF-74s. This is because the ionic radius of Mg 2+ is smaller than those of Co 2+ , Ni 2+ and Zn 2+ , which results in a larger adsorption volume in the Mg-MOF-74. More importantly, the interactions between Mg 2+ and adsorbate molecules were stronger than the counterpart interactions for Co 2+ , Ni 2+ and Zn 2+ ions [33]. The enthalpy of desorption calculated by GCMC simulations is shown in Appendix A, confirming this statement.   Figure 5 depicts the adsorption isotherms of refrigerants including R1234yf, R1234ze, R134a and R32 in M-MOF-74 at different temperatures, in which symbols represent the results measured in our GCMC simulations. Among the four refrigerants investigated, R32 was demonstrated to have the highest adsorption capacity in the M-MOF-74 within the temperature range considered, evident in Figure 5c,g,k,o, followed by R134a. This was due to the molecule structure of refrigerants. As is shown in Figure 2, the R32 molecule has a structure with a smaller effective diameter. R1234yf and R1234ze own structures with large effective diameters. At the high pressures, when the adsorption approached saturation, the adsorbate molecules experienced a strong repulsive force from the tightly-packed neighbouring molecules and the host structure [31,32]. This repulsive interaction should be further enhanced for the molecules having a larger size, leading to reduction of the adsorption volume available in MOF-74 for R1234yf, R1234ze and R134a compared to R32. Consequently, the highest adsorption capacity was achieved for R32 in the pores of MOF-74. On the other hand, the Mg-MOF-74 adsorbed more refrigerants than other MOF-74s. This is because the ionic radius of Mg 2+ is smaller than those of Co 2+ , Ni 2+ and Zn 2+ , which results in a larger adsorption volume in the Mg-MOF-74. More importantly, the interactions between Mg 2+ and adsorbate molecules were stronger than the counterpart interactions for Co 2+ , Ni 2+ and Zn 2+ ions [33]. The enthalpy of desorption calculated by GCMC simulations is shown in Appendix A, confirming this statement.

Thermal Energy Storage Properties of MOHCs
As previously mentioned, the thermal energy storage in MOHCs can be calculated according to Equation (1) or Equation (2). The relationships between thermal energy storage and temperature difference (the temperature of the cold source is 293 K) of the adsorbates (R1234yf, R1234ze, R134a, R32)/MOF-74 mixture at 3,500 and 5,000 kPa are shown in Figures 6 and 7, respectively. The thermal energy storage capacity of adsorbates (R1234yf, R1234ze, R134a, R32) increased with the mass ratio of MOF-74 nanoparticles in the working fluid. The percentage enhancement of the thermal energy was calculated as [∆hMOHCs(∆T) − ∆hFluid(∆T)]/∆hFluid(∆T), where ∆T is the temperature difference relative to the reference temperature of 293 K.
In the M-MOF-74 (M = Co, Mg, Ni) structures, the energy storage enhancement ratio of R1234yf was higher than that of the other three adsorbates. However, in Zn-MOF-74, R134a achieved the highest energy storage enhancement ratio among the four refrigerants considered. The energy

Thermal Energy Storage Properties of MOHCs
As previously mentioned, the thermal energy storage in MOHCs can be calculated according to Equation (1) or Equation (2). The relationships between thermal energy storage and temperature difference (the temperature of the cold source is 293 K) of the adsorbates (R1234yf, R1234ze, R134a, R32)/MOF-74 mixture at 3500 and 5000 kPa are shown in Figures 6 and 7, respectively. The thermal energy storage capacity of adsorbates (R1234yf, R1234ze, R134a, R32) increased with the mass ratio of MOF-74 nanoparticles in the working fluid. The percentage enhancement of the thermal energy was calculated as [∆h MOHCs (∆T) − ∆h Fluid (∆T)]/∆h Fluid (∆T), where ∆T is the temperature difference relative to the reference temperature of 293 K.
In the M-MOF-74 (M = Co, Mg, Ni) structures, the energy storage enhancement ratio of R1234yf was higher than that of the other three adsorbates. However, in Zn-MOF-74, R134a achieved the highest energy storage enhancement ratio among the four refrigerants considered. The energy storage enhancement ratio of R32 in Zn-MOF-74 was higher than that of the other three M-MOF-74 (M = Co, Mg, Ni) structures, while the energy storage enhancement ratios of R1234yf, R1234ze and R134a in Mg-MOF-74 were higher than those of the other three M-MOF-74 (M = Co, Ni, Zn). It is worth noting that the energy storage enhancement ratio at 5000 kPa is higher than that of 3500 kPa. However, the negative enhancement of thermal energy storage properties is found in a certain temperature range in Figure 6c,g,i,j,k,l,o and Figure 7c,k,l,o. This is because the enthalpy difference of the refrigerant was larger than the sum of the thermodynamic energy change of M-MOF-74 and desorption heat at the corresponding temperature range. . It is worth noting that the energy storage enhancement ratio at 5,000 kPa is higher than that of 3,500 kPa. However, the negative enhancement of thermal energy storage properties is found in a certain temperature range in Figures 6c,g,i,j,k,l,o and 7c,k,l,o. This is because the enthalpy difference of the refrigerant was larger than the sum of the thermodynamic energy change of M-MOF-74 and desorption heat at the corresponding temperature range.  Tempreture difference (K)  Tempreture difference (K)

Conclusions
In this work, molecular simulations including MD and GCMC simulations were employed to investigate the energy storage by the adsorption of adsorbates (R1234yf, R1234ze, R134a, R32) in MOF-74. From the results, the following conclusions could be drawn.
Among the four different refrigerants, R32 has the largest adsorption capacity in M-MOF-74 materials, followed by R134a. This is due to the molecular structures of R32 and R134a being smaller than those of R1234yf and R1234ze. Mg-MOF-74 adsorbs more refrigerant than the other three MOFs, as the adsorption interaction between Mg-MOF-74 and refrigerants is stronger than the other MOFs. The

Conclusions
In this work, molecular simulations including MD and GCMC simulations were employed to investigate the energy storage by the adsorption of adsorbates (R1234yf, R1234ze, R134a, R32) in MOF-74. From the results, the following conclusions could be drawn.
Among the four different refrigerants, R32 has the largest adsorption capacity in M-MOF-74 materials, followed by R134a. This is due to the molecular structures of R32 and R134a being smaller than those of R1234yf and R1234ze. Mg-MOF-74 adsorbs more refrigerant than the other three MOFs, as the adsorption interaction between Mg-MOF-74 and refrigerants is stronger than the other MOFs.
The energy storage enhancement ratio of R1234yf with M-MOF-74 (M = Co, Mg, Ni) nanoparticles was higher than other refrigerants with M-MOF-74 (M = Co, Mg, Ni) nanoparticles. The energy storage enhancement ratio of R134a with Zn-MOF-74 nanoparticles was enhanced compared to other refrigerants with Zn-MOF-74 nanoparticles; while the energy storage enhancement ratio of R32 in Zn-MOF-74 was higher than those of the other three M-MOF-74 (M = Co, Mg, Ni) structures. The energy storage enhancement ratios of R1234yf, R1234ze and R134a in Mg-MOF-74 were higher than those of the other three M-MOF-74 (M = Co, Ni, Zn). Besides, the negative enhancement of thermal energy storage properties is found under certain conditions of MOHCs.