Molecular dynamics simulation of thermal properties of modified graphene / n-octadecane composite phase change material

Phase change materials have been widely used in building energy-saving and off-peak energy storage systems, but most phase change materials limit their application because of their low thermal conductivity. Improving their heat transfer performance and revealing the heat transfer mechanism from a microscopic perspective are the keys to practical applications. The PCM system composed of n-octadecane and modified graphene was established by molecular dynamics simulation and the influence of macro thermal properties was analyzed from a microscopic perspective including end-to-end distance and torsional angle distribution, radial distribution function, and self-diffusion coefficient. The results show that the thermal conductivity of PCMS increases with increased modified graphene content. When the mass fraction is 26.88%, the thermal conductivity increases by 44.5%, which is consistent with the trend of experimental values. The addition of modified graphene resulted in a more concentrated distribution of n-octadecane molecules, indicating that the arrangement of n-octadecane molecules was affected, increasing the thermal conductivity of n-octadecane. The deviation of phase transition temperature from the experimental value is less than 1%.


Introduction
Solar energy, as a non-polluting and endless supply of energy, is considered to be the best choice to solve energy problems at present [1].However, due to the intermittent and unstable energy supply of solar energy, thermal energy storage is considered one of the key technologies for the development of solar energy [2,3].In addition, thermal energy is a low-grade energy source, which is often regarded as 2 waste heat and difficult to use effectively [4].Energy storage technology is a way to reduce the unevenness between energy supply and demand through energy storage, thereby improving efficiency and reducing the cost of renewable energy technologies [5].Among various thermal energy storage methods, latent heat storage based on phase change materials (PCMs) is the most effective technology [6,7].PCM has many excellent characteristics, such as extremely high heat storage density, strong chemical stability, and low cost [8,9].It can be used in industrial waste heat recovery, solar energy systems, and construction [10][11][12][13].PCM has almost no temperature change during phase change, and can be used for temperature control [14][15][16].
The low thermal conductivity is the most common disadvantage of most PCMs, which limits the thermal response rate of thermal energy storage systems [17].Carbon-based materials, including graphene [18,19], expanded graphite [20], and multi-walled carbon nanotubes [21,22], are the most promising high thermal conductivity fillers [23].Due to their good thermal conductivity, they can be used to enhance the heat transfer of phase change heat storage materials.
Graphene has extremely high intrinsic thermal conductivity, large specific surface area, and excellent chemical stability.Based on previous research, it was confirmed that adding nanographene sheets to n-octadecane can effectively enhance its thermal conductivity.When the mass fraction of nanographene sheets is 2wt.%, the thermal conductivity of phase change materials increases by 89.4% [24].
Although composite PCMs have good application prospects, cyclic stability is still a problem [25].In composite PCM systems, graphene materials will agglomerate and settle due to the van der Waals force, resulting in a decrease in thermal conductivity.Although surfactants can be added to inhibit agglomeration [26], they may affect the physical properties of composite PCMs.Some studies have strengthened the dispersibility in EG by grafting functional groups on the edges of graphene [27].Our previous experiments have shown that adding modified graphene (MG) to n-octadecane can greatly improve the energy storage efficiency of PCMs [28].
As an effective method for the microscopic analysis of materials, molecular dynamics methods have made great progress in the study of the thermophysical properties of materials.With the development of molecular dynamics methods, some meaningful explorations have appeared in the microscopic research of alkane PCMs.Sosso et al. used molecular dynamics software to simulate the rapid crystallization process of composite PCMs [29].Yin et al. used molecular dynamics simulation to calculate the heat capacity and self-diffusion coefficient of the n-hexadecane system, and the freezing point of n-hexadecane was obtained accordingly [30].Babaei et al. used molecular dynamics simulations to study the phase transition process of PCM systems with carbon nanotubes (CNTs) or graphene, and calculated the thermal conductivity of the system using balanced molecular dynamics simulations.The results shows that CNTs make the thermal conductivity of PCMs increased by 66% (solid) and 48% (liquid), graphene increased the thermal conductivity of PCM by 87% and 52% [31].Rao et al. calculated the phase change behavior of alkane phase change materials, and the simulation results were close to the experiments, which indicating that the addition of water molecules into the alkane molecules would reduce the system fluidity and increase the heat storage capacity [32].Zabihi et al. found that the functional groups of hydrogenation, methyl and phenyl on the surface of graphene were helpful to improve the thermal conductivity.As the coverage of H, CH3 and C6H5 functional groups on the surface of graphene increased, the thermal resistance of the graphene-paraffin interface decreased, and the thermal conductivity of the nanocomposite increased [33].Wang et al found that the hydrogenation of graphene can greatly reduce the interfacial thermal resistance.With the increase of hydrogen atom coverage, the relative interfacial thermal resistance decreases [34].Luo et al. used non-equilibrium molecular dynamics methods to study the factors affecting the thermal resistance of the interface between graphene and polymer matrix.They proposed that increasing the polymer density, increasing the size of graphene, and forming covalent bonds are conducive to interfacial heat transfer [35].
In this paper, we will study the thermophysical properties of the composite system of modified graphene (MG) and n-octadecane through molecular dynamics simulation and compare the simulation results with the experimental test results.

Structure models and simulation methods
Due to the limitation of simulation conditions, the MG model used in the simulation is idealized and simplified.A supercell structure was established to obtain a single-layer graphene sheet containing 160 carbon atoms.Four functional groups are connected at the edges of the graphene sheet near the four corners.A model of the MG is shown in Fig 1 .The gray atoms in the picture represent carbon atoms, the white atoms represent hydrogen atoms, and the red atoms are oxygen atoms.The Amorphous cell module of Materials Studio software was used to construct the initial model of the MG/n-octadecane composite PCM.In order to reflect the difference in the amount of MG, a total of 6 models have been established, in which the number of MG molecules is 0, 1, 2, 3, 4, 5, respectively.For the convenience of description, these sets of models were labeled as MG0, MG1, MG2, MG3, MG4 and MG5, respectively.At the same time, in order to better simulate the distribution of MG in noctadecane, the MG was fixed parallel to each other in the center of the simulation box.After the initial model is built, there are some irrationals in the system, so it is necessary to perform a relaxation process and optimize the structure to obtain a system as close to the actual situation as possible.The relaxation process includes energy minimization, annealing, and kinetic equilibrium processing.After the relaxation process is completed, from the energy and temperature curves of the entire process, we can see that the system is basically stable, as shown in Fig 3.After the system reached equilibrium, the COMPASS force field was selected to start dynamic calculations at different temperatures between 280K and 330K.Within the set temperature interval, calculations are performed every 5K intervals.The 1st World Energy Materials Conference Journal of Physics: Conference Series 2749 (2024) 012008

3.1Results of thermal conductivities of composites
The thermal conductivities of six groups of MG composite PCM systems were calculated using the RNEMD method [36], the results are shown in Fig 4 .As the temperature increases, the thermal conductivity of each system generally decreases.It can be seen that when the temperature exceeds 300K, the thermal conductivity decreases significantly.Taking MG0 as an example, the phase transition temperature of n-octadecane is 28℃ (301K), and the region where the thermal conductivity significantly decreases can be considered as the phase transition occurrence region (295K-305K).The calculated results show that the average thermal conductivity of solid and liquid n-octadecane are 0.222W/(m•K) and 0.198W/(m•K), respectively.It's worth noting that the thermal conductivity increases slightly between 305 and 310K, which is the end of the phase transition.The variation of the thermal conductivity of the other groups of MG PCM systems at different temperatures is basically consistent with MG0.The decrease of the thermal conductivity and the slight increase of the thermal conductivity after the phase transition is related to the phase transition behavior of n-octadecane, which is the result of the effect of the phase transition process on the molecular order and molecular motion in the system.The mass fraction of MG in six composite PCM systems was calculated as 0%, 5.40%, 10.80%, 16.69%, 21.54% and 26.88%, respectively.The thermal conductivity calculated at the temperature of 293K is compared with the thermal conductivity measured experimentally, and the results are shown in Fig 5 .The experimental and simulation results show that the thermal conductivity of the composite PCM increases with increased MG content.However, the experimentally measured thermal conductivity increases much faster with the increased modified graphene mass fraction than the simulation results.This is because it is difficult to simulate the space heat conducting mesh structure formed by the modified graphene in the n-octadecane system under the existing simulation conditions.

3.2Discussions about results of thermal conductivities of composite
During the phase transition, when the solid material absorbs heat and melts into a liquid state, the entropy value of the system increases, and the degree of order decreases.Through molecular dynamics simulation, we can observe that with the change of temperature, the atoms in the system move continuously, and their arrangement also changes, which leads to changes in the molecular arrangement and configuration in the model.We can research the changes of order degree in the composite PCMs system by end-to-end distance and torsion angle of the n-octadecane molecular chain at different temperatures, as shown in Fig 6.  distribution of molecules with end-to-end distance less than 15 Å are more random, with no obvious regularity, means that the degree of disorder of the n-octadecane molecule is getting higher and more random.
Fig 7 is a distribution diagram of torsion angles of n-octadecane molecules at different temperatures.The peak value in the coplanar ±180° region gradually decreases as the temperature increases, and the peak value in the ±160° and ±75° regions gradually appears, indicating that the degree of twisting and bending of n-octadecane molecules increases under the action of heating up, and the consistency and order degree of molecules decrease.When the temperature exceeds the phase transition point, it can be considered that the ordered arrangement of the solid n-octadecane system and the consistency of the molecules are destroyed.According to the principle of solid heat conduction, there are two types of heat transfer inside solid materials, which are lattice heat conduction and electron heat conduction.For n-octadecane, it mainly depends on lattice conduction.The normal mode energy quantum of the lattice vibration is called a phonon.The heat conduction in the crystal depends on the motion of the phonon.The microscopic expression of its thermal conductivity is shown in formula (1) In the formula, cv represents the heat capacity per unit volume, v ̅ represents the average group velocity of the phonon, and l represents the mean free path of the phonon.
When the orderness of the n-octadecane PCM decreases due to melting, the mean free path of the phonon also decreases.According to formula (1), the thermal conductivity decreases.
In addition, during the phase transition, molecules in the solid and liquid states coexist in the noctadecane system.At the interface between the solid and liquid, the molecule's movement is hindered.After the phase transition is completed, the blocked state disappears, thermal conductivity will rise  , the end-to-end distance of the pure octadecane system is mostly distributed between 15Å and 20Å.As the amount of MG in the system increases, the peaks around 17.5Å become higher, and the peaks in other ranges gradually become shorter.In Fig 9 (b), as the amount of MG increases, the peak around ± 180 ° rises, indicating that the degree of bending of the n-octadecane molecules in these systems are getting lower, until MG5, which exceeds 80% of n-octadecane.The torsion angle of the alkane molecule is not less than ± 170°, showing a high consistency.According to the previous analysis, this is the expression of increasing order of n-octadecane molecules in these systems.In other words, the addition of MG can improve the order degree of molecules in the PCM system.

Radial distribution function
The physical meaning of the Radial Distribution Function (RDF) is the ratio of the probability of another particle appearing at a distance "r" relative to the random distribution of a reference particle.The formula is as follows [37] In the formula, N represents the number of particles, and ρ represents the density of the entire system.
We can use the RDF to reflect the aggregation and binding of molecules in the n-octadecane PCM system.The results are shown in Fig10.There are peaks at 1.125 Å and 1.145 Å, and smaller peaks near 2Å to 3Å.The higher the peak value, the higher the proportion of n-octadecane molecules maintaining the corresponding distance.The positions of the two highest peaks did not change substantially with temperature, indicating that temperature did not significantly affect the aggregation and binding mode of n-octadecane molecules.As the temperature increases, we can observe that all peaks become shorter, which means that the number of molecules that maintain the corresponding distance gradually decreases.This may be due to the higher kinetic energy of n-octadecane at higher temperatures, the greater the relative distance between molecules, and the increased mobility of the molecules in the system.At the same time, the decrease of the number of molecules keeping the same relative distance can also be interpreted as the decrease of the order degree of the n-octadecane system.molecules, but the change is not obvious.This is because only n-octadecane occurs in the phase change process of the composite PCMs, and the MG does not significantly affect the process.Since the MG was fixed in the center of the simulation box, the thermal movement of n-octadecane molecules was limited, and their mobility was reduced, resulting in a slightly higher radial distribution function.

Results and discussions of Phase change characteristics 3.4.1 Self-diffusion coefficient
The phase change characteristics of PCMs in molecular dynamics simulation can be studied by the self-diffusion coefficient.The self-diffusion coefficient is a single-valued function of temperature.It is generally used to describe the diffusion behavior of atoms or molecules in molecular dynamics simulation.Its numerical value reflects the strength of the mobility of atoms in the system.If the value of the self-diffusion coefficient changes abruptly, it indicates that the system has undergone a phase change.
Fig 12 is a curve of the mean square displacement (MSD) of the n-octadecane PCM system at different temperatures during the simulated temperature rise.With the development of the simulation, the kinetic energy of molecules in the system gradually increased, and MSD showed an increasing trend.At the same time, the higher the simulated temperature, the higher the kinetic energy obtained by the molecules in the same time, so the MSD value is higher.The slope of n-octadecane molecule MSD at each temperature is calculated.After calculation and processing, the self-diffusion coefficient of n-octadecane PCM systems in the temperature range of 280K to 330K is obtained [38,39], and the result is shown in Fig 13 .The value of the self-diffusion coefficient increases with increasing temperature, indicating that the molecular mobility in the noctadecane system increases with increasing temperature [41] .The self-diffusion coefficient scatter points of two temperature ranges from 280K to 300K and from 305K to 330K were linearly fitted, and two straight lines were obtained.The temperature at the intersection was 302.6K.This point can be regarded as the phase transition temperature of the composite PCM [32].MG promotes the crystallization of noctadecane molecules at this temperature [31].Compared with the experimentally measured phase transition temperature (301.6K) of n-octadecane, the error does not exceed 1%, which proves that the simulation results are reliable.

Specific heat capacity
In addition to the latent heat of phase change, the specific heat capacity of the PCM affects its heat storage capacity.In molecular dynamics simulations, we calculate the derivative of the total potential energy of the system with respect to temperature, and obtain the constant heat capacity CV () The potential energy curve of the n-octadecane PCM system during the temperature rise process is drawn.As shown in Fig 16, the potential energy of the system increases with increasing temperature, which approximates a linear increase.When the PCM does not undergo a phase change, the heat capacity is a single value function of temperature.If the phase change occurs, the value of the heat capacity will change greatly, but the temperature will not change.Therefore, if the curve of specific heat capacity with respect to temperature changes abruptly, it indicates that phase change occurs in the PCM.As shown in Fig 17, the specific heat capacity at 300K increased significantly, which is close to the experimental value of the phase transition temperature of n-octadecane (301.6 K), which proves that the simulation results are reliable.Although we cannot obtain the latent heat value of the n-octadecane system through simulation, we can see that the specific heat capacity value near the phase transition temperature jumps significantly due to the latent heat, which is obviously conducive to thermal energy storage.The experimental measured value is compared with the simulated specific heat capacity value.When the temperature is 293K, we measured the experimental value as 1.658J/(g•K) and the simulated value as 2.083J/(g•K).The difference between the two values may come from the COMPASS force field, which is the most suitable force field for describing condensed-phase materials, but it is also difficult to accurately describe the accurate thermodynamic behavior of complex organic compound molecules [40].Nevertheless, the simulation results can still be used as a reference for studying the thermophysical properties of PCMs.

Conclusion
In this paper, the molecular dynamics simulation was used to calculate the n-octadecane and MG composite PCM system, and the results were compared with the experimental results.The following conclusions are obtained: (1) Comparing the simulation results of the composite phase change system with the n-octadecane system, it is found that the trend of thermal conductivity with temperature changes is the same as the experimental results, and the value of thermal conductivity increases with the increase in the content of MG.
(2) Analysis of the parameters such as the end-to-end distance, torsion angle, RDF, and selfdiffusion coefficient of each system shows that the reason for the increase in thermal conductivity of the composite PCMs is that the interaction between MG and n-octadecane caused an increase in the order of n-octadecane, leading to an increase in the thermal conductivity of n-octadecane itself.This is the main factor affecting the thermal conductivity of composite PCMs in molecular dynamics simulation.
(3) Simulation results show that the phase transition temperature of PCM decreases with the addition of MG.The difference between the calculated value of the phase transition temperature and the experimental value is less than 1%.
Fig 2 is the initial model of MG3.The yellow part is n-octadecane molecule and the blue part is MG.

Fig 3 .
Fig 3.Energy variation of relaxation process (a)Energy variation of annealing process; (b)Energy variation of kinetic equilibrium process.

Fig 4 .
Fig 4.Thermal conductivity of various systems at different temperatures.

Fig 5 .
Fig 5.Comparison of experimental and simulation results of thermal conductivity.

Fig 6 .
Fig 6.Schematic diagram of end to end distance and torsion angle (a) end to end distance; (b) torsion angle.

Fig 7
Fig 7 shows the end-to-end distance distribution of the n-octadecane PCM at different temperatures.As shown in Fig 7, at 285K, the terminal distance distribution has a highest peak between 17 Å and 20 Å, indicating that at low temperatures, the arrangement of n-octadecane is relatively neat, with a high order degree.As the temperature increases, we can clearly observe that the peaks in the range of 17Å to 20 Å are getting shorter and shorter, and the highest peak has a tendency to shift to the left.The

Fig 7 .
Fig 7.Distribution curve of end to end distance of n-octadecane molecules at different temperatures during heating.
explains the increase in thermal conductivity from 305K to 310K inFig 4.

Fig 8 .
Fig 8.Distribution curve of torsion angle of n-octadecane molecules at different temperatures during heating.

Fig 9 .
Fig 9.Distribution curve of end to end distance and torsion of n-octadecane molecules in various systems. :

Fig 11
Fig 11 shows the radial distribution function of n-octadecane molecules in MG1, MG2, MG3, MG4 and MG5.As shown in Fig 11, the value of the RDF increases with the increase of the number of MGmolecules, but the change is not obvious.This is because only n-octadecane occurs in the phase change process of the composite PCMs, and the MG does not significantly affect the process.Since the MG was fixed in the center of the simulation box, the thermal movement of n-octadecane molecules was limited, and their mobility was reduced, resulting in a slightly higher radial distribution function.

Fig 11 .
Fig 11.RDF of n-octadecane molecules in each system at 325K.

Fig 12 .
Fig 12.Mean square displacement curve of n-octadecane system at different temperatures during heating.

Fig 14 .
Fig 14.MSD of n-octadecane molecules in each system at 325K.

Fig 15 Fig 15 .
Fig 15 is a scatter diagram of the self-diffusion coefficient of each group of MG/PCM systems with the change of temperature.The phase transition temperature can be obtained by fitting the scatter points to a straight line.As shown in Fig 15, all the systems follow the rule that the self-diffusion coefficient increases with the increase of temperature.The phase transition temperature decreased slightly with the increase of the content of the MG (Fig 15(b) , and the overall trend slightly increased compared to Fig 15(a)), and the magnitude of the change was small, which was consistent with the experimental conclusion.This may mean that some interactions between the MG and n-octadecane caused such changes.

Fig 16 .
Fig 16.Potential energy curve of n-octadecane system during heating.

Fig 17 .
Fig 17.Specific heat capacity variation of n-octadecane system with temperature.

Fig 18 .
Fig 18.Specific heat capacity variation of each system with temperature.