Molecular Dynamics Simulation of the Interaction between R1336mzz(Z) and POE Lubricants

1  Introduction

The Organic Rankine Cycle (ORC) systems [1,2] employ the organic working fluid as an energy carrier, facilitating the transformation of low-grade thermal energy into high-grade electrical energy. The ORC system presents numerous advantages, such as a straightforward design, ease of maintenance, high safety, and environmental sustainability. ORC technology [3] is based on the principles of the Rankine cycle, utilizing organic substances with low boiling points as the working fluid. The lower boiling points of these organic compounds enable their heating via low-temperature heat sources, facilitating the generation of high-pressure steam for the expansion machine to perform work, ultimately leading to power output. Different from the traditional water vapor Rankine cycle power generation technology, ORC technology can mainly be applied to medium and low temperature heat source conditions. In the ORC, the working fluid functions as a thermal carrier and is pivotal in influencing the efficacy of the the variable temperature heat source [4]. The choice of working fluid in ORC systems is one of the critical and significant issues investigated by researchers and scholars.

Hydrocarbons (HC), hydrofluorocarbons (HFC), hydrochlorofluorocarbons (HCFC), and chlorofluorocarbons (CFC) are potential working fluids suitable for ORC systems [5]. In the selection of an appropriate working fluid, it is crucial to prioritize not only the achievement of high cycle efficiency within the ORC system but to address considerations related to environmental protection and safety. Recent research has extensively examined the physical property parameters of ORC systems. However, in response to increasingly stringent environmental protection standards, the international community has designated R113 and R11 refrigerants previously employed in the initial phases of ORC development as phased out working substances. This decision is attributed to their high global warming potential (GWP) and ozone depletion potential (ODP), despite their favorable energy conversion efficiency [6,7]. In alignment with the amendments stipulated in the Montreal Protocol, developed nations have pledged to phase out refrigerants possessing high global warming potential within the period from 2019 to 2036 [8]. As a result, a new generation of working fluids, known as the fourth generation, has been developed. This includes compounds such as R1233zd(E) and R1336mzz(Z) [9,10], which are intended to replace hydrofluorocarbon (HFC) fluids.

Traditionally, experimental data on the thermal physical properties of working fluids are recorded or plotted in thermodynamic properties diagrams or tables for query. With the development of computers, experimental data on the physical properties of working fluid are fitted using computers according to the developed physical properties model formula. Where, some of the fluid physical property models have been developed specifically for one or some fluids, such as the IAPWS-95 [11] and IAPWS-IF97 [12] formulas for water and water vapor, and the GERG-2004 [13] and GERG-2008 [14] formula for natural gas mixtures. Some are developed for multiple or one class of fluids, such as Helmholtz energy equation of state and Modified Benedict-Webb-Rubin (MBWR) state equation [15,16], etc. These thermophysical property models have been developed into computer software, among which the more popular ones include NIST REFPROP [17,18], CoolProp [19] and Helmholtzmedia [20], etc. The use of working mass physical properties calculation and query software to obtain working mass physical properties is the most widely used method in the industry. Take the Helmholtz energy equation of state as an example, which is established using two binary functions, whose the two independent variables are independent state parameters: temperature and density; two dependent variables are non-dimension ideal gas construction Helmholtz energy and non-dimension residual contribution to Helmholtz energy. All thermodynamic property parameters can be obtained by the analytic derivatives of these two dependent variables. In practice, temperature and density are selected as their non-dimensional relative values to the critical temperature and critical density. This method fits the Helmholtz energy function entirely on the basis of experimental data, focusing on the macroscopic experimental data fitting, without the study of microscopic particles, and can not predict the thermal physical property parameters of the new working medium. Accordingly, the Helmholtz energy equation of state for R1336mzz(Z) [21] has been developed and applied to the NIST REFPROP. This is the most mature and accurate model for calculating the thermal physical properties of R1336mzz(Z) in industry. Furthermore, molecular dynamics (MD) simulations are also frequently employed to investigate the thermophysical properties of refrigerants. Although the quantum thermodynamic microscopic nature of temperature is not well defined, molecular dynamics gives a classical thermodynamic definition of temperature at the microscopic scale. The 26th General Conference on Weights and Measures (CGPM, Conférence Générale des Poids et Mesures) held by the Bureau International des Poids et Mesures (BIPM) in 2018 revised the International System of Units (SI, Système International d’Unités) [22], in which the Boltzmann constant was defined as a fixed value, establishing a fixed relationship between temperature and energy, and the dimension of temperature was derived from energy. Therefore, the relationship between the dynamic energy of microscopic molecules and the macroscopic temperature is fixed and clear. This definition clarified the microscopic nature of temperature in classical thermodynamics and improved the theoretical basis of molecular dynamics. In classical physics, molecular dynamics can link the microscopic motion of molecules with the macroscopic thermodynamic properties, and the study of the microscopic motion of molecules is helpful to predict the macroscopic thermophysical properties of new organic working fluids. Khan et al. [23] use dynamics simulation to study the macroscopic thermodynamic properties of R1336mzz(Z), and the simulation results were compared with NIST REFPROP. The results were almost identical, which also proved that molecular dynamics simulation is a quite accurate and mature method to study R1336mzz(Z) fluid.

R1336mzz(Z) [24] exhibits an Ozone Depletion Potential (ODP) of 0, a Global Warming Potential (GWP) of merely 2, is non-combustible, and demonstrates excellent compatibility with commonly utilized lubricating oils. It is considered one of the most promising working fluids for ORC applications, garnering significant scholarly attention [9,24] Molés et al. [25] conducted comparative experiments between R245fa and the novel working fluid R1336mzz(Z) within an ORC system. The findings indicated that the system employing R1336mzz(Z) achieved a higher thermal efficiency than the system using R245fa, with efficiency improvements ranging from −0.3% to 17%. Li et al. [26] and Tanaka et al. [27] conducted in-depth studies on the saturated vapor pressure and density of R1336mzz(Z), respectively. Other physical parameters, such as critical temperature [28], viscosity [29], and thermal conductivity [30], have been investigated experimentally. In a study by Hu et al. [31], MD simulations were utilized to examine the energy storage characteristics of refrigerants, specifically R1234yf, in conjunction with metal-organic framework (MOF) nanoparticles. This research underscores the utility of molecular dynamics in analyzing both physical and thermodynamic properties. The homogeneous condensation, thermal properties, and structural characteristics of R450A, R513A, and R515A were examined utilizing the MD method [32,33]. In contrast, the majority of research concerning the physical properties of R1336mzz(Z) has primarily concentrated on experimental measurements, with a significant paucity of studies analyzing its physical property parameters at the nanoscale.

The thermal physical properties of the working fluid are recognized as critical prerequisites for its processing. In ORC systems, the lubricating oil inevitably interacts with the refrigerant, resulting in notable alterations to the thermal physical properties of refrigerant upon mixing, which subsequently impacts system performance [34,35]. Consequently, it is imperative to investigate the interactions between HFC refrigerants, such as R1336mzz(Z) and lubricating oil. Polyol Ester Oil (POE) [36] is one of the most extensively utilized lubricants in refrigeration and heat pump systems, primarily due to its remarkable thermal stability and lubricative properties. POE demonstrates compatibility with a variety of refrigerants characterized by low GWP. POE lubricants have high temperature stability, high lubricity, economic friendliness, and is widely used in engine and low temperature heat engine mechanical structure lubrication. Lubricants are used in heat engines to lubricate mechanical structures. The working fluid in the ORC expands and does work in the turbine, and is condensed by the pump pressure back to the high pressure state, while the lubricant is the lubrication of the mechanical structure, and the working fluid is inevitably mixed with the lubricating oil in every link. The fundamental composition of POE comprises polyols and ester compounds synthesized from organic acids, with polyols functioning as the principal constituents of POE oil. Pentaerythritol esters (PEC) are the most frequently used polyols, extensively applied in lubricants for refrigeration and air conditioning systems due to their superior lubrication properties and thermal stability [37].

This study aims to bridge the existing research gap concerning the interaction mechanisms between the fourth-generation refrigerant R1336mzz(Z) and lubricating oil at the nanoscale by utilizing molecular dynamic simulation to explore the interactions between R1336mzz(Z) and PEC. Additionally, it examines the variations in physical properties within pure R1336mzz(Z) and the mixed system of R1336mzz(Z)/PEC. The findings provide theoretical support for the application and optimal design of R1336mzz(Z) in practical ORC systems.

2  Simulation Detail

All molecular simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package [38,39]. Figs. 1 and 2 depict the molecular structures of R1336mzz(Z) and PEC, with PEC serving as the primary component of the POE. ORC systems working condition necessitates that OPE exhibits both high thermal stability and low volatility. Although longer PEC chains generally confer enhanced thermal stability and reduced volatility, they also tend to exhibit decreased solubility [40]. Consequently, for this simulation, PEC with medium-long chain lengths was selected. Medium-long chain PEC are typically characterized by ester side chains containing between six and ten carbon atoms. Within this range, PECs exhibit an optimal balance of properties, offering effective lubrication performance while maintaining appropriate viscosity and thermal stability [41]. As shown in Fig. 2, we select PEC9 with the number of ester side chain C atoms of 9. We use CVFF force fields [42] for refrigerant R1336mzz(Z) and PEC9. The interaction between the atoms is described by the Lennard-Jones 12-6 potential as Formula (1).

E=4ε((εr)12−(σr)6)(1)

where, ε and σ represent the characteristic energy and van der Waals radius, respectively. The detailed parameters of the LJ potential of each atom are shown in Table 1. The arithmetic combining rule is adopted for the interaction parameters between atoms. In addition, the cutoff for the electrostatic part of the potential and the van der Waals force is 10.0 Å. The long-distance electrostatic force and van der Waals interactions are treated by the particle-particle-particle-grid (PPPM) [43] method with a RMS accuracy of 10−4.

Figure 1: R1336mzz(Z) molecular

Figure 2: PEC9 molecular

In the first part of the work, the pure refrigerant R1336mzz(Z) was simulated. As shown in Fig. 3a, we build a simulation box containing 2250 R1336mzz(Z) molecules that have periodic boundary changes in the x, y and z directions. Under the initial conditions, all 2250 molecules are arranged repeatedly in a simple cubic structure with an initial velocity of zero. At the beginning of the simulation, the software directly changes the molecular velocity until a specified equilibrium temperature is reached. First, the pure R1336mzz(Z) system is preliminarily equilibrated under the NVT ensemble (constant N, constant V, constant T), and the molecular dynamics process of 500ps is performed at 270, 300, 330, and 360K, respectively. Subsequently, a 3ns system equilibrium process was carried out under the NPT ensemble (constant N, constant P, constant T). The NVT ensemble holds the volume of the simulation box constant until equilibrium is reached. The NPT ensemble continuously adjusts the volume of the simulation box to achieve the specified equilibrium pressure. The temperature is 270, 300, 330, 360 K, and its pressure is set to the saturation pressure at these four simulated temperatures. After the system reaches equilibrium, the relevant parameters are statistically measured within 500 ps, such as radial distribution function (RDF), mean square displacement (MSD), diffusion coefficient. In this simulation system, the timestep is set to 1 fs.

Figure 3: System model

In the second part of the work, in order to reveal the mechanism of interaction between R1336 and lubricating oil, we mixed 8 and 18 PEC9 in 2250 R1336mzz(Z) moleculars systems respectively, to form two mixing systems with different lubricating oil content. This mixing system is shown in Fig. 3b. The same temperature and pressure as the pure R1336mzz(Z) system were simulated. Each system first relaxes 500 ps under the NVT ensemble, then reaches equilibrium within 3 ns simulation time under the NPT ensemble, the timestep of the first 1 ns simulation is 0.1 fs, and the timestep of the second 2 ns simulation is 1 fs, and finally performs data acquisition through the dynamic process of 500 ps.

3  Result Discussion

First, In this section, the simulated density is compared with the date of NIST REFPROP to verify the correctness of the model. Second, the molecular interaction mechanism of R1336mzz(Z) and PEC9 is discussed from four aspects: radial distribution function, mean azimuth-shift, diffusion coefficient, binding energy and solubility parameter. The effect of the addition of PEC9 on the diffusion of R1336mzz(Z) and the change of compatibility between the two at different temperatures were analyzed. This provides a theoretical basis for compatibility between refrigerants and lubricants.

3.1 Model Verification

Density is a fundamental characteristic of the working medium, primarily influenced by its temperature and pressure conditions [44]. On a microscopic scale, variations in molecular interaction forces is related to changes in density closely. Consequently, the precision of density measurements serves as an indicator of the model’s validity, and accurate density values are crucial for ensuring the reliability of simulation outcomes. Therefore, prior to conducting any simulations, we need to validate the constructed simulation model. For this purpose, we select ten pressure points at a constant temperature of 300 K, with pressures ranging from 1 to 2.8 atm, over a duration of 3 ns. And the mean density of R1336mzz(Z) within the balanced system was determined. This calculated density was then compared to the density data from NIST REFPROP 10.0, as illustrated in Fig. 4, where 4 a shows the comparison between the calculated density of R1336 and the NIST REFPROP 10.0 density. The analysis of the figure indicates that the computed density aligns closely with the experimentally determined density. Fig. 4b illustrates the relative error between these two densities, revealing that the relative errors are minimal. Notably, the largest relative error occurs at 1.2 atm; however, it remains below 5%, which is considered acceptable within the established parameters [45].

Figure 4: Comparison of simulated and NIST REFPROP 10.0 data for liquid R1336mzz(Z) density

3.2 Interaction Mechanism between R1336mzz(Z) and PEC9

3.2.1 Radial Distribution Function

The Radial Distribution Function (RDF) quantifies the spatial distribution and interactions between molecules or atoms in a system. The arrangement of molecules, the range of interactions between molecules, and the density change of the system at different distances are revealed. As shown in Formula (2).

G(r)=14πr2ρ0⋅(dNdr)(2)

ρ0 is the mean density of the system, that is, the global mean density of the particles, and dN is the number of particles in a spherical shell of radius r and thickness dr number. The calculation of the RDF for specific atomic pairs provides insights into the microstructure and intermolecular interactions within the refrigerant. For both the pure R1336mzz(Z) system and the mixing systems, the RDFs between F and C atoms on the molecular backbone, C atoms on the side chain and H atoms, as well as between H and F atoms, have been calculated. The statistical outcomes of these calculations are presented in Figs. 5–7.

Figure 5: The RDF of each atom pairs at different temperatures in Pure R1336mzz(Z)

Figure 6: The RDF of each atom pairs at different temperatures in mixing system, where PEC9 chain number equal to 8

Figure 7: The RDF of each atom pairs at different temperatures in mixing system, where PEC9 chain number equal to 18

In the first diagram of Fig. 5, for the pure refrigerant system, the RDF between the skeleton C and F atoms exhibits three distinct peaks. The first peak is observed near 3 Å, while the third peak is notably higher than the first two. This phenomenon can be attributed to the molecular geometry of R1336mzz(Z), which leads to a weaker short-range direct interaction between the F atom and the skeleton C atom, thereby resulting in lower amplitudes for the first two peaks. At greater distances, however, molecular aggregation occurs in a manner that facilitates significant intermolecular arrangements, thereby increasing the likelihood of C-F pair occurrences at the third peak. As the first diagram of Figs. 6 and 7 show, in the mixing systems, the initial peak is observed near 3Å; however, only two peaks are present. The incorporation of PEC9 potentially induces greater chaos or disorder within the system, resulting in a less distinct distribution of C and F atoms at certain distances. This phenomenon contributes to the attenuation or even the disappearance of the peaks. In the second diagram of Figs. 5–7, the RDF exhibits a peak corresponding to the interaction between the side chain C and H atoms. Specifically, the peak occurs at 3 Å in the pure refrigerant system and shifts to approximately 2.7 Å in the mixing systems. This shift can be attributed to the substantial size and intricate spatial configuration of PEC9 molecules, which occupy significant volume within the mixed system. Consequently, this spatial occupation facilitates a closer arrangement of R1336mzz(Z) molecules, thereby reducing the distance between the side chain C and H atoms. The RDF between H and F provides insights into polar interactions. Due to high electronegativity of F atom, it acquires a substantial negative charge, while H atoms exhibit a partial positive charge. By examining the spatial distribution of H and F atoms through RDF analysis, one can predict the presence of weak hydrogen bonds or electrostatic attractions within or between molecular entities. In the third diagram of Fig. 5, pertaining to a pure R1336mzz(Z) system, the initial peak is observed near 2.5 Å. This suggests that interactions between H and F atoms, potentially including electrostatic interactions or van der Waals forces, are most pronounced at this distance. Consequently, H atoms are more likely to be in proximity to F atoms at this range. The significant interaction between H and F atoms within this distance indicates a distinct local polar interaction, which is analogous to the characteristics of weak hydrogen bonds. In the second diagram of Figs. 5–7, the initial peak is observed near 2.5 Å; however, only three distinct peaks are present. This suggests that PEC9 occupies a specific spatial region upon mixing, altering the molecular distances from their original peak values. The polar groups of PEC9 exhibit strong interactions with the H or F atoms of R1336mzz(Z), thereby disrupting the intermolecular arrangement of R1336mzz(Z). Consequently, this interference results in the elimination of some peaks, effectively obliterating the original four specific distances. Furthermore, the data presented in Figs. 5–7 indicate that the atomic pairs of R1336mzz(Z) exhibit comparable radial distribution function trends across varying temperatures. Specifically, for the C and F atomic pairs within the molecular framework, peak intensities increase with rising temperatures across all three systems. Conversely, for the C-H and H-F atomic pairs, temperature exerts a significant influence in the pure refrigerant system, whereas the mixed system demonstrates a reduced sensitivity to temperature variations.

Furthermore, the RDFs of atomic pairs between R1336mzz(Z) and PEC9 were also calculated. The RDF of side chain C atom and the O atom in PEC9 can reflect the relative position and contact probability between R1336mzz(Z) and PEC9 molecules, revealing the local arrangement of the two in the liquid phase. The F atom exhibits high electronegativity and interacts with the O atom within the PEC9 molecule. The RDF of the F-O atom pair can elucidate short-range electrostatic interactions or potential hydrogen bonding effects between these two entities. Additionally, the RDF of the H-O atom pairs can indicate the presence of weak hydrogen bond interactions, which is instrumental in comprehending the interactions between R1336mzz(Z) and PEC9. As illustrated in Figs. 8 and 9, various mixing systems exhibit a similar trend, with the peak value approximately at 4 Å. Variations in temperature can result in differences in the peak value. Conversely, the quantity of PEC9 chains exerts minimal influence on this parameter. Notably, the RDF does not display a significant peak value for the F atoms of R1336mzz(Z) and the O atoms of PEC9. Despite the high electronegativity of the F atom, its chemical reactivity and interaction potential are contingent upon the molecular environment. The molecular configuration and side chain geometry of R1336mzz(Z) likely displace F atoms away from the O atoms in PEC9, leading to an absence of significant pairing distribution in the RDF. The RDF between the H atoms of R1336mzz(Z) and the O atoms in PEC9 peaks at a distance of 4.5 Å, which is typically beyond the range associated with hydrogen bonding or strong dipole-dipole interactions. This observation implies that the interaction between H and O is predominantly characterized by van der Waals forces or weak polar interactions, rather than a direct strong attraction.

Figure 8: The RDF of atomic pairs between R1336mzz(Z) and pentaerythrito, in mixing system, where PEC9 chain number equal to 8

Figure 9: The RDF of atomic pairs between R1336mzz(Z) and pentaerythrito in mixing system, where PEC9 chain number equal to 18

3.2.2 Mean Squared Displacement and Diffusion Coefficient

Mean Squared Displacement (MSD) and the Diffusion Coefficient (D) are two critical physical parameters characterizing molecular motion. MSD quantifies the average squared displacement of a molecule from its initial position over a specified time interval, and the diffusion coefficient measures the rate at which a molecule diffuses within a system, serving as a metric for the extent of molecular dispersion over time. MSD is typically defined as Formula (3).

MSD(t)=(r(t)−r(0))2¯(3)

The relation between D and MSD is as Formula (4).

D=16limt→∞dMSD(t)dt(4)

where, r(t) represents the position of the molecule at a given time t, and r0 is the position at the initial time. Fig. 10 illustrates the MSD curves for pure R1336mzz(Z) and mixing systems across different temperatures. The figure clearly demonstrates that the MSD of these systems is influenced by temperature. At 360 K, the MSD exhibits a rapid increase over time, indicating enhanced molecular mobility. Conversely, at 270 K, the MSD progresses more slowly, suggesting restricted molecular movement. The D presented in Fig. 11 further illustrates this phenomenon. Specifically, the diffusion coefficient at 360 K surpasses those observed under the other three conditions, suggesting a greater degree of molecular dispersion within the liquid phase. This indicates that the motion of R1336mzz(Z) intensifies with rising temperature. Furthermore, Fig. 11 reveals that within the pure R1336mzz(Z) system, the diffusion coefficients at various temperatures exceed those of the mixed system. Additionally, the two systems with differing PEC9 numbers exhibit similar results.

Figure 10: The MSD of the different system at different temperatures

Figure 11: The D of the different system at different temperatures

3.2.3 Binding Energy

Binding energy represents a fundamental physical parameter that quantifies the strength of the interaction between two molecules. Typically articulated as the energy released during the molecular binding process, a higher value of binding energy indicates a more stable interaction between the molecules. The binding energy between the refrigerant R1336mzz(Z) and PEC9 can be determined using the Formula (5).

Ebinding=ER1336mzz(Z)+EPEC9−ER1336mzz(Z),PEC9(5)

where, Ebinding is the binding energy between PEC9 and R1336mzz(Z) and its value is equal to the opposite number of the interaction energy. ER1336mzz(Z),PEC9 is the total energy of the R1336mzz(Z)/PEC9 mixing system, ER1336mzz(Z) and EPEC9 are the energies of the R1336mzz(Z) and PEC9, respectively.

As illustrated in Fig. 12, the binding energy between R1336mzz(Z) and PEC9 within the mixing system was calculated, primarily to quantify the non-bonding energy interactions. The findings indicate that the binding energy diminishes as temperature rises, attributable to the intensification of molecular thermal motion at elevated temperatures. The non-bonding interaction between molecules (such as van der Waals force and electrostatic force), are typically contingent upon the relative distance and orientation of the interacting molecules. As the temperature rises, the extent of molecular thermal motion expands, complicating the maintenance of a stable relative arrangement among molecules and thereby weakening the intermolecular forces. Furthermore, an increase in the quantity of PEC9 polymer leads to a significant rise in the system’s binding energy, suggesting that the number of polymers positively influences the enhancement of intermolecular interactions. This finding elucidates the correlation between the interaction of R1336mzz(Z) and PEC9 with temperature and the concentration of PEC9 polymer within the system. Furthermore, it offers a theoretical foundation for enhancing the compatibility between refrigerant and lubricant.

Figure 12: Binding energy between R1336mzz(Z) and PEC9 in mixed systems at different temperatures

3.2.4 Solubility Parameter

The solubility parameter (δ) is a critical physical quantity employed to quantitatively assess intermolecular compatibility, primarily derived from the energy distribution of intermolecular interactions within a system. Fundamentally, the solubility parameter is defined based on the cohesive energy density of the system. It can be calculated using the Formula (6).

δ=EcohV(6)

where, Ecoh represents the cohesive energy density of the system, defined as the cohesive energy per unit volume, while V denotes the total volume of the system. We have calculated the solubility parameters of the refrigerant and pentaerythritol ester PEC9 within the mixing system, analyzed the temperature dependence and the disparity in solubility parameters between the two components, and assessed their compatibility. Table 2 and Fig. 13 present the solubility parameters of R1336mzz(Z) and PEC9 within the mixing system, as well as the difference between the solubility parameters of the two components, respectively. The data in the table illustrate the temperature dependence of the solubility parameters. Specifically, an increase in temperature corresponds to a decrease in the solubility parameters of both components in the mixing system. Moreover, the difference in solubility parameters can indicate the compatibility between two components. A larger disparity in these parameters typically corresponds to a more pronounced difference in the interaction forces between the components, potentially resulting in altered compatibility within the mixture. In practical applications involving refrigerants and lubricants, such differences can significantly impact the thermodynamic properties and stability of the mixture. The solubility parameter difference between two components is inversely related to the similarity of their intermolecular forces; a smaller difference indicates more similar intermolecular forces and thus a stronger propensity for dissolution and mixing. Conversely, a larger solubility parameter difference signifies more pronounced differences in intermolecular forces, resulting in reduced compatibility. Analysis of the solubility difference diagram reveals that the compatibility of the two components diminishes as temperature increases. Additionally, an increase in the quantity of PEC9 correlates with a greater solubility parameter difference, worse the compatibility.

Figure 13: Difference of solubility parameters between R1336mzz(Z) and PEC9

4  Conclusion

In this study, the interaction mechanism between the fourth-generation refrigerant R1336mzz(Z) and PEC, the principal component of lubricating oil, was investigated using molecular dynamics simulations. Specifically, PEC with a chain length of nine (PEC9) was selected for analysis. The parameters, including the RDF of relevant atomic pairs, MSD, diffusion coefficient, binding energy, and solubility parameter were calculated. A comparative analysis of the diffusion coefficients of R1336mzz(Z) in both pure and mixing systems was conducted. The findings indicate that the incorporation of PEC9 impedes the diffusion of the refrigerant. Furthermore, an analysis was conducted on the variations in binding energy and solubility parameters between the refrigerant and PEC9 across different mixing systems. The findings indicate that the binding energy of the mixed system containing a higher concentration of PEC9 polymer chains is greater than that of systems with fewer polymer chains. Additionally, the binding energy decreases as the temperature rises. The analysis of solubility parameter differences reveals that this difference increases with rising temperature, and the results also indicate that compatibility deteriorates as the PEC9 content increases.

Based on the mechanism of microscopic molecular interaction, the macroscopic thermodynamic parameters such as molar heat capacity at constant volume, molar heat capacity at constant pressure, compression factor and thermodynamic energy can be calculated by molecular dynamics method in future studies. Compared with direct experiments, this is a faster and quite accurate method to obtain hydrodynamic parameters for new mixed working fluids that are not included in physical property calculation software or not developed by traditional thermodynamic property calculation models. Additionally, it is helpful to optimize the efficiency and power of ORC to analyze the thermodynamic properties of the phase transition process of the mixed working medium. The thermodynamic properties of non-azeotropic mixtures used in ORC can also be predicted by molecular dynamics method. The choice of mixing medium in ORC is an important factor affecting its performance. It is a promising innovation direction to combine molecular dynamics prediction of working medium physical properties with ORC thermal cycle simulation. This enables the development of a mature mixture selection scheme for optimizing ORC.

Acknowledgement: None.

Funding Statement: This work was supported by Hainan Provincial Natural Science Foundation of China (No. 422CXTD509).

Author Contributions: The authors confirm contribution to the paper as follows: Study conception and design: Haoyuan Jing, Zhongye Wu; Data collection: Haoyuan Jing, Xiaoyang Jiang; Analysis and interpretation of results: Haoyuan Jing, Xiaoyang Jiang; Valuable guidance on the simulation aspects: Qingfen Ma; Draft manuscript preparation: Haoyuan Jing, Zhongye Wu, Xiaoyang Jiang. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

References

1. Park BS, Usman M, Muhammad I, Pesyridis A. Review of organic rankine cycle experimental data trends. Energy Conver Manag. 2018;173:679–91. doi:10.1016/j.enconman.2018.07.097. [Google Scholar] [CrossRef]

2. Quoilin S, Van Den Broek M, Declaye S, Dewallef P, Lemort V. Techno-economic survey of organic rankine cycle (ORC) systems. Renew Sustain Ener Rev. 2013;22(3):168–86. doi:10.1016/j.rser.2013.01.028. [Google Scholar] [CrossRef]

3. Lecompte S, Huisseune H, van den Broek M, Vanslambrouck B, De Paepe M. Review of organic rankine cycle (ORC) architectures for waste heat recovery. Renew Sustain Ener Rev. 2015;47:448–61. doi:10.1016/j.rser.2015.03.089. [Google Scholar] [CrossRef]

4. Bao J, Zhao L. A review of working fluid and expander selections for organic rankine cycle. Renew Sustain Ener Rev. 2013;24(8):325–42. doi:10.1016/j.rser.2013.03.040. [Google Scholar] [CrossRef]

5. Ziviani D, Beyene A, Venturini M. Advances and challenges in orc systems modeling for low grade thermal energy recovery. Appl Ener. 2014;121(15):79–95. doi:10.1016/j.apenergy.2014.01.074. [Google Scholar] [CrossRef]

6. Dixon RK. Global environment facility investments in the phase-out of ozone-depleting substances. Mitigat Adaptat Strat Global Change. 2011;16(5):567–84. doi:10.1007/s11027-011-9281-2. [Google Scholar] [CrossRef]

7. Sicard AJ, Baker RT. Fluorocarbon refrigerants and their syntheses: past to present. Chem Rev. 2020;120(17):9164–303. doi:10.1021/acs.chemrev.9b00719. [Google Scholar] [CrossRef]

8. Albà CG, Alkhatib III, Llovell Fèlix, Vega LF. Assessment of low global warming potential refrigerants for drop-in replacement by connecting their molecular features to their performance. ACS Sustain Chem Eng. 2021;9(50):17034–48. doi:10.1021/acssuschemeng.1c05985. [Google Scholar] [PubMed] [CrossRef]

9. Yang J, Ye Z, Yu B, Ouyang H, Chen J. Simultaneous experimental comparison of low-GWP refrigerants as drop-in replacements to R245fa for organic rankine cycle application: R1234ze(zR1233zd(eand R1336mzz(e). Energy. 2019;173(1):721–31. doi:10.1016/j.energy.2019.02.054. [Google Scholar] [CrossRef]

10. E.Elahi A, Mahmud T, Alam M, Hossain J, Biswas BN. Exergy analysis of organic rankine cycle for waste heat recovery using low Gwp refrigerants. Int J Thermofluids. 2022;16(15):100243. doi:10.1016/j.ijft.2022.100243. [Google Scholar] [CrossRef]

11. Wagner W, Pruß A. The LAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J Phys Chem Ref Data. 2002;31(2):387–535. doi:10.1063/1.1461829. [Google Scholar] [CrossRef]

12. Wagner W, Kruse A. IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. Berlin/Heidelberg: Springer; 1998. p. 7–37. doi:10.1007/978-3-662-03529-0_3. [Google Scholar] [CrossRef]

13. Kunz O, Group EGR editors. The GERG-2004 wide-range equation of state for natural gases and other mixtures. Als ms. gedr ed. No. 557 in Fortschritt-Berichte VDI Reihe 6, Energietechnik. VDI-Verl; 2007. [cited 2025 Feb 20]. Available from: https://www.osti.gov/etdeweb/biblio/20924249. [Google Scholar]

14. Kunz O, Wagner W. The GERG-2008 wide-range equation of state for natural gases and other mixtures: an expansion of GERG-2004. J Chem Eng Data. 2012;57(11):3032–91. doi:10.1021/je300655b. [Google Scholar] [CrossRef]

15. Jacobsen RT, Stewart RB. Thermodynamic properties of nitrogen including liquid and vapor phases from 63 K to 2000 K with pressures to 10,000 bar. J Phys Chem Ref Data. 1973;2(4):757–922. doi:10.1063/1.3253132. [Google Scholar] [CrossRef]

16. Outcalt SL, McLinden MO. Equations of state for the thermodynamic properties of R32 (difluoromethane) and R125 (pentafluoroethane). Int J Thermophys. 1995;16(1):79–89. doi:10.1007/BF01438959. [Google Scholar] [CrossRef]

17. Lemmon EW, Bell IH, Huber ML, McLinden MO. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0, National Institute of Standards and Technology; 2018. [cited 2025 Feb 20]. doi:10.18434/T4/1502528. [Google Scholar] [CrossRef]

18. Huber ML, Lemmon EW, Bell IH, McLinden MO. The nist refprop database for highly accurate properties of industrially important fluids. Indust Eng Chem Res. 2022;61(42):15449–72. doi:10.1021/acs.iecr.2c01427. [Google Scholar] [PubMed] [CrossRef]

19. Bell IH, Wronski J, Quoilin S, Lemort V. Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolprop. Indus Eng Chem Resh. 2014;53(6):2498–508. doi:10.1021/ie4033999. [Google Scholar] [PubMed] [CrossRef]

20. Thorade M, Saadat A. Helmholtzmedia-a fluid properties library. In: Linköping Electronic Conference Proceedings; 2012; Linköping, Sweden: Linköping University Electronic Press. p. 63–70. doi:10.3384/ecp1207663. [Google Scholar] [CrossRef]

21. McLinden MO, Akasaka R. Thermodynamic properties of cls-1,1,1,4,4,4-Hexafluorobutene [r-1336mzz(z)]: vapor pressure, (p, ρ, T) behavior, and speed of sound measurements and equation of state. J Chem Engg Data. 2020;65(9):4201–14. doi:10.1021/acs.jced.9b01198. [Google Scholar] [PubMed] [CrossRef]

22. Bureau international des poids et mesures. Proceedings of the 26th meeting of the CGPM; 2018 [cited 2025 Feb 20]. Available from: https://www.bipm.org/committees/cg/cgpm/26-2018. [Google Scholar]

23. Khan M, Wen J, Shakoori MA, Liu Y. Condensation and thermophysical properties of R1336mzz(z) through molecular dynamics simulations. Int J Refrig. 2023;146:290–9. doi:10.1016/j.ijrefrig.2022.11.008. [Google Scholar] [CrossRef]

24. Dawo F, Fleischmann J, Kaufmann F, Schifflechner C, Eyerer S, Wieland C, et al. R1224yd(zR1233zd(e) and R1336mzz(z) as replacements for R245fa: experimental performance, interaction with lubricants and environmental impact. Appl Ener. 2021;288:116661. doi:10.1016/j.apenergy.2021.116661. [Google Scholar] [CrossRef]

25. Molés F, Navarro-Esbrí J, Peris B, Mota-Babiloni A, Barragán-Cervera Á, Kontomaris KK. Low GWP alternatives to HFC-245fa in organic rankine cycles for low temperature heat recovery: HCFO-1233zd-E and HFO-1336mzz-Z. Appl Therm Eng. 2014;71(1):204–12. doi:10.1016/j.applthermaleng.2014.06.055. [Google Scholar] [CrossRef]

26. Li S, Xu L, Liu H, Yang Z, Duan Y. Vapor pressure measurements and correlation for cls-1,1,1,4,4,4-hexafluoro-2-Butene (HFO-1336mzz(z)). J Chem Eng Data. 2020;65(9):4223–9. doi:10.1021/acs.jced.0c00024. [Google Scholar] [CrossRef]

27. Tanaka K, Akasaka R, Sakaue E, Ishikawa J, Kontomaris KK. Thermodynamic properties of cls-1,1,1,4,4,4-hexafluoro-2-butene (HFO-1336mzz(z)measurements of the Pρ t property and determinations of vapor pressures, saturated liquid and vapor densities, and critical parameters. J Chem Eng Data. 2016;61(7):2467–73. doi:10.1021/acs.jced.6b00169. [Google Scholar] [CrossRef]

28. Tanaka K, Akasaka R, Sakaue E, Ishikawa J, Kontomaris KK. Measurements of the critical parameters for cls-1,1,1,4,4,4-Hexafluoro-2-Butene. J Chem Eng Data. 2017;62(3):1135–8. doi:10.1021/acs.jced.6b00990. [Google Scholar] [CrossRef]

29. Sun Y, Li X, Meng X, Wu J. Measurement and correlation of the liquid density and viscosity of Hfo-1336mzz(z) (cis-1,1,1,4,4,4-Hexafluoro-2-Butene) at high pressure. J Chem Eng Data. 2019;64(2):395–403. doi:10.1021/acs.jced.8b00713. [Google Scholar] [CrossRef]

30. Perkins RA, Huber ML. Measurement and correlation of the thermal conductivity of cls-1,1,1,4,4,4-Hexafluoro-2-Butene. Int J Thermophys. 2020;41(7):103. doi:10.1007/s10765-020-02681-0. [Google Scholar] [PubMed] [CrossRef]

31. Hu J, Liu C, Liu L, Li Q. Thermal energy storage of R1234yf, R1234ze, R134a and R32/mof-74 nanofluids: a molecular simulation study. Materials. 2018;11(7):1164. doi:10.3390/ma11071164. [Google Scholar] [PubMed] [CrossRef]

32. Khan M, Wen J, Shakoori MA, Tao W. Homogeneous condensation and thermophysical properties of R450a, R513a and R515a using molecular dynamics simulations. J Mol Liq. 2022;353(7):118795. doi:10.1016/j.molliq.2022.118795. [Google Scholar] [CrossRef]

33. Xu Z, Wang Z, Xia X, Li X, Chen Y, Yi Q. Molecular dynamics simulation of homogeneous condensation and thermophysical properties of R245fa/R141b. Appl Therm Eng. 2024;236(19):121627. doi:10.1016/j.applthermaleng.2023.121627. [Google Scholar] [CrossRef]

34. Conde MR. Estimation of thermophysical properties of lubricating oils and their solutions with refrigerants: an appraisal of existing methods. Appl Therm En. 1996;16(1):51–61. doi:10.1016/1359-4311(95)00011-2. [Google Scholar] [CrossRef]

35. Zhang Y, Yang Z, Feng B, Shi Y, Lv Z, Chen Y. Study on the inhibition of lubricating oil in conjunction with a flame retardant on the flammability of propane. Combust Flame. 2023;255:112867. doi:10.1016/j.combustflame.2023.112867. [Google Scholar] [CrossRef]

36. Marsh KN, Kandil ME. Review of thermodynamic properties of refrigerants + lubricant oils. Fluid Phase Equilib. 2002;199(1–2):319–34. doi:10.1016/S0378-3812(02)00025-0. [Google Scholar] [CrossRef]

37. Emel'ianov VV, Krasnykh EL, Portnova SV, Levanova SV. Synthetic oils based on pentaerythritol esters. vapor pressure and enthalpy of vaporization. Fuel. 2022;312(1):122908. doi:10.1016/j.fuel.2021.122908. [Google Scholar] [CrossRef]

38. Thompson AP, Aktulga HM, Berger R, Bolintineanu DS, Brown WM, Crozier PS, et al. LAMMPS—a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comput Phys Commun. 2022;271(4):108171. doi:10.1016/j.cpc.2021.108171. [Google Scholar] [CrossRef]

39. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys. 1995;117(1):1–19. doi:10.1006/jcph.1995.1039. [Google Scholar] [CrossRef]

40. Sarpal AS, Sastry MIS, Kumar R, Bhadhavath S, Rai K, Bansal V, et al. Molecular dynamics of synthetic-based lubricant system by spectroscopic techniques—part 1. Tribol Trans. 2013;56(3):442–52. doi:10.1080/10402004.2011.651772. [Google Scholar] [CrossRef]

41. Han ZH, Yu YD. Selection of working fluids for low-temperature power generation organic rankine cycles system. Adv Mater Res. 2012;557–559:1509–13. doi:10.4028/www.scientific.net/AMR.557-559.1509. [Google Scholar] [CrossRef]

42. Dauber-Osguthorpe P, Roberts V, Osguthorpe D, Wolff J, Genest M, Hagler A. Structure and energetics of ligand binding to proteins: escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system. Proteins. 1988;4(1):31–47. doi:10.1002/prot.340040106. [Google Scholar] [PubMed] [CrossRef]

43. Hockney RW, Goel SP, Eastwood JW. Quiet high-resolution computer models of a plasma. J Comput Phys. 1974;14(2):148–58. doi:10.1016/0021-9991(74)90010-2. [Google Scholar] [CrossRef]

44. Thol M, Lemmon EW. Equation of state for the thermodynamic properties of trans-1,3,3,3-tetrafluoropropene [r-1234ze(e)]. Int J Thermophys. 2016;37(3):28. doi:10.1007/s10765-016-2040-6. [Google Scholar] [CrossRef]

45. Dang T, Nguyen H. A study on the simulation and experiment of evaporative condensers in an R744 air conditioning system. Micromachines. 2023;14(10):1826. doi:10.3390/mi14101826. [Google Scholar] [PubMed] [CrossRef]