Molecular dynamics-based method for predicting thermophysical properties of polar surface liquid membranes

Abstract

The application discloses a polar surface liquid film thermal property prediction method based on molecular dynamics, relates to the technical field of microcosmic thermal property simulation, and comprises the following steps: constructing an initial simulation system containing a polar solid surface and polar liquid molecules, obtaining a stable system through temperature and pressure control and energy minimization, and running to thermodynamic equilibrium; extracting molecular trajectory data from the equilibrium system, calculating density distribution and molecular orientation distribution of a liquid film region, and obtaining microstructure characterization data; based on the same trajectory data, applying a Green-Kubo formula to calculate a heat flow autocorrelation function of the liquid film region, and obtaining a thermal conductivity prediction value of the region through time integral processing; and combining the microstructure data and the thermal conductivity prediction value, and analyzing an influence mechanism of a polar surface effect on liquid film thermal properties.

Description

Technical Field

[0001] This invention belongs to the field of microscopic thermal property simulation technology, specifically a method for predicting the thermal properties of polar surface liquid films based on molecular dynamics. Background Technology

[0002] Accurate prediction of the thermal properties, especially thermal conductivity, of extremely thin liquid films near solid surfaces is a key challenge in micro- and nano-scale heat transfer research. Macroscopic experimental methods are limited by spatial resolution, making it difficult to directly measure the thermal properties of nanoscale liquid films. While conventional molecular dynamics simulations provide atomic-scale research tools, they typically employ simplified force field models when dealing with such problems, neglecting atomic details and real electrostatic interactions. This leads to severe distortions in the description of key microstructures such as the arrangement of liquid molecules and hydrogen bond networks at the interface, failing to accurately reflect the complex dipole-dipole and charge-dipole interactions between polar molecules and polar surfaces.

[0003] Existing molecular simulation techniques also face difficulties in calculating the thermal conductivity of such non-homogeneous systems. Common methods either calculate the bulk thermal conductivity of the entire simulation chamber, which is diluted by the substrate and a large amount of bulk liquid, failing to represent the true properties of the interfacial liquid film; or employ non-equilibrium molecular dynamics methods to apply a temperature gradient in the system, but constructing a stable linear temperature field in the extremely thin liquid film region is very difficult, and strong temperature gradients may induce non-physical effects. How to directly and reliably extract the heat transport coefficient corresponding only to the specific spatial domain of the liquid film region from the equilibrium simulation is a bottleneck that existing conventional techniques have failed to effectively solve. Therefore, there is a need to develop a molecular simulation method that can accurately describe polar interfacial interactions and spatially resolve local thermophysical properties. Summary of the Invention

[0004] This invention aims to solve at least one of the technical problems existing in the prior art; Therefore, this invention proposes a method for predicting the thermal properties of polar surface liquid films based on molecular dynamics, including: Construct an initial molecular dynamics simulation system that includes polar solid surfaces and polar liquid molecules; Temperature and pressure controls were applied to the initial molecular dynamics simulation system, and energy minimization was performed to generate a stable pre-equilibrium system. The pre-equilibrium system is run until it reaches thermodynamic equilibrium, forming an equilibrium molecular dynamics simulation system for sampling; Extract molecular coordinates and velocity trajectory data for a specified time period from the equilibrium molecular dynamics simulation system; Based on the molecular coordinates and velocity trajectory data, the density distribution and molecular orientation distribution of the liquid film region are calculated to generate liquid film microstructure characterization data. Based on the molecular coordinates and velocity trajectory data, the heat flux autocorrelation function of the liquid film region is calculated using the Green-Kubo formula. The heat flux autocorrelation function is integrated over time to obtain the predicted thermal conductivity of the liquid film region; By combining the microstructure characterization data of the liquid film with the predicted thermal conductivity, the influence mechanism of polar surface effects on thermal properties is analyzed.

[0005] Furthermore, the construction of the initial molecular dynamics simulation system comprising a polar solid surface and polar liquid molecules includes: Based on the lattice parameters and surface polarity characteristics of the target polar solid surface, a solid surface cell model under periodic boundary conditions is established. Based on the force field parameter file of the target polar liquid molecule, determine the atomic charge distribution, bonding and non-bonding interaction parameters of the polar liquid molecule; Above the solid surface cell model, polar liquid molecules are arranged with a specified thickness and initial density to form a liquid region; A vacuum layer is provided above the liquid region to eliminate spurious interactions of periodic boundary conditions in the direction perpendicular to the surface; The solid surface cell model, the liquid region, and the vacuum layer are combined to form an initial molecular dynamics simulation system with three-dimensional periodicity.

[0006] Furthermore, the step of applying temperature and pressure control and performing energy minimization on the initial molecular dynamics simulation system to generate a stable pre-equilibrium system includes: A position constraint algorithm is applied to all atoms in the initial molecular dynamics simulation system to prevent system collapse caused by an unreasonable initial configuration; The conjugate gradient algorithm is invoked to search the potential energy surface of the initial molecular dynamics simulation system, find the local potential minimum point, and complete the energy minimization process; The algorithm releases the positional constraints on atoms in the liquid region while keeping atoms in the solid region fixed or applying flexible constraints. A constant temperature control algorithm is used to adjust the system temperature to the target temperature value, and a pressure control algorithm is used to adjust the pressure to the target pressure value; Under these conditions, a short-term molecular dynamics pre-run is performed to allow the system to initially relax and generate the stable pre-equilibrium system.

[0007] Furthermore, the process of running the pre-equilibrium system until it reaches thermodynamic equilibrium, forming an equilibrium molecular dynamics simulation system for sampling, includes: Long-term molecular dynamics simulations were performed on the pre-equilibrium system under isothermal and isobaric or isothermal and isochoric conditions. Real-time monitoring of the potential energy, temperature, pressure, and time-history changes of the system dimensions of the pre-equilibrium system; When the average values ​​of the time-history changes in potential energy, temperature, pressure, and system size remain stable and the fluctuation amplitude is less than a preset threshold, the system is determined to have reached thermodynamic equilibrium. The continuous simulated trajectory after reaching thermodynamic equilibrium is extracted, and the corresponding system state constitutes the equilibrium molecular dynamics simulation system used for sampling.

[0008] Further, the extraction of molecular coordinates and velocity trajectory data for a specified time period from the equilibrium molecular dynamics simulation system includes: A time length much longer than the molecular vibration period is set as the specified time period; During the molecular dynamics simulation, the three-dimensional coordinate vectors and three-dimensional velocity vectors of all atoms in the equilibrium molecular dynamics simulation system are recorded at fixed time intervals. The recorded three-dimensional coordinate vectors and three-dimensional velocity vectors are organized in chronological order to form the molecular coordinate and velocity trajectory data.

[0009] Further, the step of calculating the density distribution and molecular orientation distribution of the liquid film region based on the molecular coordinates and velocity trajectory data to generate liquid film microstructure characterization data includes: The liquid film region is divided into a series of thin layers of equal thickness along a direction perpendicular to the surface of the polar solid; The density distribution is calculated by statistically analyzing the average change in the number of atoms or centroids within each thin layer over time. For polar liquid molecules, calculate the dipole moment vector for each molecule; Based on the angle between the dipole moment vector and the normal vector of the polar solid surface, the statistical distribution of molecular orientation within each thin layer is analyzed to obtain the molecular orientation distribution; The density distribution and the molecular orientation distribution are combined to form the microstructure characterization data of the liquid film.

[0010] Furthermore, the step of calculating the autocorrelation function of the heat flux in the liquid film region using the Green-Kubo formula based on the molecular coordinates and velocity trajectory data includes: Based on the molecular coordinates and velocity trajectory data, the total potential energy of the equilibrium molecular dynamics simulation system and the force on each atom at each time step are calculated using force field parameters. Using the atomic velocities, total potential energy, and forces acting on the atoms from the molecular coordinates and velocity trajectory data, the instantaneous heat flux vector at each time step is calculated using the heat flux definition formula; Extract the heat flow component parallel to the polar solid surface from the instantaneous heat flow vector; The autocorrelation average value of the heat flow component parallel to the polar solid surface at different time intervals is calculated to obtain the heat flow autocorrelation function.

[0011] Further, the step of performing time integration on the heat flux autocorrelation function to obtain the predicted thermal conductivity value of the liquid film region includes: Integrate the heat flow autocorrelation function from time zero to a sufficiently long correlation cutoff time; Integration stops when the heat flow autocorrelation function decays to zero or fluctuates slightly around its mean after the correlation cutoff time. Substituting the integral result into the thermal conductivity calculation formula, and combining it with the system's volume and temperature parameters, the predicted thermal conductivity of the liquid film region in the direction parallel to the surface is calculated.

[0012] Furthermore, the analysis of the influence mechanism of polar surface effects on thermal properties by combining the microstructure characterization data of the liquid film with the predicted thermal conductivity includes: The density distribution curve and molecular orientation distribution curve in the microstructure characterization data of the liquid film are correlated with the predicted thermal conductivity value. Identify the location and intensity of the oscillation peaks that appear near the surface of the polar solid in the density distribution curve; Identify the preferred orientation angles of molecules in specific regions as shown in the molecular orientation distribution curve; By analyzing the correspondence between the oscillation peak position, the oscillation peak intensity, the preferred orientation angle, and the predicted thermal conductivity, the correlation between the interface layer structure and the heat transfer efficiency can be inferred.

[0013] Furthermore, after analyzing the correspondence between the oscillation peak position, the oscillation peak intensity, the preferred orientation angle, and the predicted thermal conductivity value, the method further includes: By changing the simulation condition parameters, a series of molecular dynamics simulation systems with different condition parameters can be constructed; For each molecular dynamics simulation system with different condition parameters, the entire process from energy minimization to thermal property analysis is repeated to obtain a series of corresponding predicted thermal conductivity values ​​and liquid film microstructure characterization data. Establish a quantitative relationship graph between the simulation condition parameters, the liquid film microstructure characterization data, and the predicted thermal conductivity value.

[0014] Compared with the prior art, the beneficial effects of the present invention are: By constructing a polar solid surface model with atomic-level details and parameterizing polar liquid molecules using force fields that accurately describe electrostatic interactions, this method fully preserves key physical properties at the interface, such as dipole moments and partial charges. This construction approach allows the simulated system to spontaneously form microstructures consistent with reality, including interfacial electric double layers, specific adsorption, and ordered molecular layers. The resulting density and molecular orientation distribution data directly reveal how polar surfaces influence the molecular arrangement and hydrogen bond dynamics of near-wall liquids, providing reliable atomic-level data for understanding how interfacial forces modulate microstructures—an effect unattainable with simplified or nonpolar models.

[0015] By employing the Green-Kubo formula and decomposing the heat flux vector into spatial regions to specifically calculate the autocorrelation function of heat flux in the liquid film region, this method directly extracts the spatially localized thermal conductivity from the equilibrium molecular dynamics trajectory. This approach avoids the challenge of constructing macroscopic temperature gradients within nanoscale gaps and eliminates the averaging interference of the bulk liquid phase on the calculation results. The resulting predicted thermal conductivity directly corresponds to the liquid film region itself affected by the surface, more realistically reflecting the enhancement or inhibition of heat transport by confinement effects and interfacial interactions. Correlation analysis of this predicted value with the aforementioned microstructure characterization data allows for the quantitative establishment of structure-property relationships between structural order parameters such as molecular arrangement and orientational order and changes in thermal conductivity, thus clearly elucidating the microscopic physical mechanisms by which polar surfaces influence thermal properties. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the steps of the method for predicting the thermal properties of polar surface liquid films based on molecular dynamics as described in this invention. Figure 2 Flowchart for the initial simulation system construction; Figure 3 Flowchart for the pre-equilibrium treatment of the system; Figure 4 This is a monitoring curve of the energy minimization process in a molecular dynamics simulation system. Figure 5 This is a graph showing the correlation between liquid film thickness and predicted thermal conductivity. Detailed Implementation

[0017] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0018] See Figure 1An initial molecular dynamics simulation system was constructed, simultaneously incorporating a polar solid surface and polar liquid molecules. Temperature and pressure were controlled, and energy minimization was applied to this initial system to generate a structurally stable pre-equilibrium system. This pre-equilibrium system was run until it reached thermodynamic equilibrium, forming a sampleable equilibrium molecular dynamics simulation system. Molecular coordinate and velocity trajectory data over a specified time period were extracted from this equilibrium system. Based on these trajectory data, the density distribution and orientation distribution of polar molecules along the direction perpendicular to the surface in the liquid film region were calculated, generating characterization data describing the liquid film's microstructure. Simultaneously, based on the same trajectory data, the autocorrelation function of the heat flux vector in the liquid film region was calculated using the Green-Kubo formula. Time integration was performed on this autocorrelation function to obtain the predicted value of the liquid film's thermal conductivity. Finally, the obtained liquid film microstructure characterization data and predicted thermal conductivity values ​​were comprehensively analyzed to explore the influence mechanism of the polar solid surface on the thermal properties of the liquid film.

[0019] In one embodiment of the present invention, see [reference] Figure 2 The specific method for constructing an initial molecular dynamics simulation system containing a polar solid surface and polar liquid molecules is as follows: A solid surface unit cell model under periodic boundary conditions is established based on the lattice parameters and surface polarity characteristics of the target polar solid surface. The target polar solid surface is an α-quartz (001) crystal plane, and its lattice parameters are obtained from a crystallography database. The lattice constant... and Both are 0.491 nanometers. The surface polarity is 0.540 nm, and the surface polarity is characterized by the presence of suspended hydroxyl groups at the surface terminals. The surface cell model is constructed by periodically replicating and expanding along the lattice vector direction in three-dimensional space to form a solid substrate model with a clear surface orientation and periodic boundary conditions. The solid substrate model is periodic in the xy plane parallel to the surface and has a definite surface boundary in the z direction perpendicular to the surface. Based on the force field parameter file of the target polar liquid molecule, the atomic charge distribution, bonding and non-bonding interaction parameters of the polar liquid molecule are determined. The target polar liquid molecule is a water molecule, and the selected force field is the TIP4P / 2005 force field model. The force field parameter file clearly specifies the atomic type and mass of oxygen and hydrogen atoms, as well as the charge of -1.1128e for each oxygen atom and +0.5564e for each hydrogen atom. The bonding interaction parameters include the equilibrium bond length and bond strength of the OH bond and the equilibrium bond angle and angular strength of the HOH bond. The non-bonding interaction parameters include the σ and ε parameters of the Lennard-Jones potential of the oxygen atom. Hydrogen atoms only participate in electrostatic interactions.

[0020] In some embodiments, polar liquid molecules are arranged above a solid surface cell model at a specified thickness and initial density to form a liquid region. The specified thickness of the liquid region is 5 nanometers, and the initial density is set to 33.3 molecules per cubic nanometer, referencing the bulk density of the TIP4P / 2005 water model at 300 Kelvin. The arrangement process uses the filling tool of a molecular simulation software package to initialize the positions and orientations of water molecules within a specified spatial region above the solid surface, either randomly or according to a simple cubic lattice, generating a liquid region configuration containing thousands of water molecules. A vacuum layer is set above the liquid region to eliminate spurious interactions of periodic boundary conditions in the direction perpendicular to the surface. The thickness of the vacuum layer is set to be no less than 2 nanometers, such that the minimum distance between the liquid films of the upper and lower periodic mirror images is greater than twice the cutoff radius of the intermolecular forces.

[0021] Optionally, the solid surface unit cell model, liquid region, and vacuum layer are combined to form an initial molecular dynamics simulation system with three-dimensional periodicity. This combination is performed in the simulation software by defining the dimensions of a three-dimensional periodic simulation box. The dimensions of the simulation box in the x and y directions are strictly equal to the periodic length of the solid surface unit cell model, and the dimension in the z direction is equal to the sum of the thickness of the solid surface unit cell model, the thickness of the liquid region, and the thickness of the vacuum layer. The constructed solid atomic coordinates and liquid molecular coordinates are placed within this simulation box, and the system is declared to follow three-dimensional periodic boundary conditions. This generates an initial simulation system coordinate file that can be read and calculated by the molecular dynamics program. It can be understood that the lattice vector of the solid surface unit cell model is defined by the lattice constant and the Miller index, expressed as:

[0022] Where: symbol The lattice vector representing the unit cell model of a solid surface, with the symbol... , , It is the expansion factor of the supercell selected when constructing the surface model in the three lattice directions, with the symbol... , , It is the lattice constant of the target polar solid material, with the symbol... , , is the unit cell basis vector of the crystal, and this expression is used to calculate the actual size of the solid substrate in the periodic direction that is ultimately used for molecular dynamics simulations.

[0023] In one embodiment of the present invention, see [reference] Figure 3The specific method for generating a stable pre-equilibrium system by applying temperature and pressure control and minimizing energy in the initial molecular dynamics simulation system is as follows: A position constraint algorithm is applied to all atoms in the initial molecular dynamics simulation system to prevent system collapse caused by unreasonable initial configurations. The position constraint algorithm is achieved by applying a rigidly fixed or highly elastic harmonic potential to the atomic positions, freezing the positions of all atoms on their initial coordinates in the initial stage of the simulation. The conjugate gradient algorithm is then used to search the potential energy surface of the initial molecular dynamics simulation system to find the local potential minimum point and complete the energy minimization process. The conjugate gradient algorithm iteratively calculates the direction of the force on the atoms and moves the atoms along the conjugate direction to reduce the total potential energy of the system. The iterative process continues until the change in the total potential energy of the system is continuously less than the set tolerance of 1 kilojoule per mole multiple times.

[0024] In some embodiments, the algorithm for releasing the positional constraints on atoms in the liquid region simultaneously keeps atoms in the solid region fixed or applies flexible constraints. Water molecules in the liquid region are allowed to move freely, while α-quartz substrate atoms in the solid region remain completely fixed in their lattice positions to simulate a rigid substrate. Alternatively, flexible constraints are applied to the atoms in the solid region by anchoring them to a harmonic oscillator potential near their equilibrium positions. The force constant of the harmonic oscillator potential is set to 1000 kJ / mol per square nanometer to allow for minute thermal vibrations in the solid region atoms. A isothermal control algorithm is used to adjust the system temperature to the target temperature value, and a pressure control algorithm is used to adjust the pressure to the target pressure value. The isothermal control algorithm uses a Nose-Hoover thermostat with a target temperature of 300 Kelvin and a relaxation time of 0.1 picoseconds. The pressure control algorithm uses a Parrinello-Rahman pressure bath with a target pressure of 0.1 MPa and a relaxation time of 1 picosecond. Isotropic pressure coupling is applied in all directions of the simulation box.

[0025] Optionally, under these conditions, a short molecular dynamics pre-run is performed to allow the system to initially relax and generate a stable pre-equilibrium system. The pre-run duration is set to 50 picoseconds, with an integration step size of 1 femtosecond. It is executed under the joint control of the Nose-Hoover thermal bath and the Parrinello-Rahman pressure bath. During the pre-run, the total potential energy, temperature, and pressure parameters of the system gradually stabilize from initial violent fluctuations. The system dimensions are moderately adjusted under pressure coupling to adapt to the set target pressure. Liquid molecules rearrange themselves near the solid surface, thus obtaining a pre-equilibrium system configuration and corresponding state file with preliminary optimization in both energy and structure. It can be understood that the Nose-Hoover thermal bath adjusts the system's kinetic energy by introducing a virtual thermal bath variable ξ. Its equation of motion extends the Hamiltonian of the system. The evolution of the thermal bath variable ξ is driven by the difference between the target temperature T and the instantaneous system temperature. The evolution equation of the thermal bath variable ξ is:

[0026] Where: symbol The derivative of the heat bath variable ξ with respect to time is represented by the sign ξ. It is an effective mass parameter related to the thermal bath coupling strength, with the symbol... Represents the sum of the kinetic energies of all N atoms in the system, with the symbol... It is the number of degrees of freedom of the system, with the symbol... It is the Boltzmann constant, symbol That is the target temperature value.

[0027] The specific method for running the pre-equilibrium system until it reaches thermodynamic equilibrium to form an equilibrium molecular dynamics simulation system for sampling is as follows: A long-term molecular dynamics simulation is performed on the pre-equilibrium system under isothermal and isobaric or isothermal and isochoric ensembles. Based on the pre-equilibrium system, a long-term molecular dynamics simulation of 5 nanoseconds is continued under the NPT ensemble to fully eliminate the influence of the initial configuration and traverse the phase space. The temporal changes of potential energy, temperature, pressure, and system dimensions of the pre-equilibrium system are monitored in real time. During the simulation, the total potential energy, instantaneous temperature, instantaneous pressure, and the length values ​​of the simulation box in the x, y, and z directions are recorded every 100 steps. The system is considered to have reached thermodynamic equilibrium when the average values ​​of the temporal changes of potential energy, temperature, pressure, and system dimensions remain stable and the fluctuation amplitude is less than a preset threshold. The preset threshold is set as follows: the standard deviation of the average value fluctuation of the physical quantity in its trajectory in the last nanosecond is less than 2% of the average value, and the fluctuation amplitude of the simulation box size in the last nanosecond is less than 0.5%. The system state corresponding to the continuous simulation trajectory after reaching thermodynamic equilibrium is used to form the equilibrium molecular dynamics simulation system for sampling. The simulation trajectory of the last 2 nanoseconds is used as the equilibrium molecular dynamics simulation system trajectory for subsequent analysis and sampling. During this time period, all physical quantities of the system remain stable in a statistical sense, and the system is in a state of full equilibrium.

[0028] In one embodiment of the present invention, the specific method for extracting molecular coordinate and velocity trajectory data for a specified time period from an equilibrium molecular dynamics simulation system is as follows: a time length much longer than the molecular vibration period is set as the specified time period. The molecular vibration period is referenced to the stretching vibration period of the OH bond in water molecules, which is approximately 10 femtoseconds. Therefore, the length of the specified time period is set to 10 nanoseconds, which is also much longer than the diffusion characteristic time of water molecules. During the molecular dynamics simulation, the three-dimensional coordinate vectors and three-dimensional velocity vectors of all atoms in the equilibrium molecular dynamics simulation system are recorded at fixed time intervals. The time integration step used in the simulation is 1 femtosecond, and the trajectory recording time interval is set to record once every 10 integration steps. That is, every 10 femtoseconds, the coordinate and velocity information of all atoms in the system, including α-quartz substrate atoms and oxygen and hydrogen atoms of water molecules, are written into the trajectory file. The recorded three-dimensional coordinate vectors and three-dimensional velocity vectors are organized in chronological order to form molecular coordinate and velocity trajectory data. The trajectory data is stored in binary or plaintext format frame by frame. Each frame contains the identifier, coordinate, and velocity information of all atoms at one time step. Continuous multiple frames of data constitute a complete molecular coordinate and velocity trajectory dataset for subsequent analysis.

[0029] In some embodiments, the specific method for generating liquid film microstructure characterization data by calculating the density distribution and molecular orientation distribution of the liquid film region based on molecular coordinates and velocity trajectory data is as follows: The liquid film region is divided into a series of thin layers of equal thickness along a direction perpendicular to the polar solid surface. Using the average position of the α-quartz solid surface as a reference, the entire liquid region is divided into thin layers with a thickness of 0.01 nanometers along the surface normal direction. The density distribution is calculated by statistically analyzing the average number of atoms or centroids in each thin layer over time. For the water molecule system, the number of oxygen atoms appearing in each water molecule within each thin layer is counted. The statistical results of all time frames in the trajectory are averaged to obtain the average number of oxygen atoms in each thin layer. The oxygen atom number density distribution curve along the normal direction is calculated by combining the cross-sectional area and thickness of the thin layer. The dipole moment vector of each molecule is calculated for the polar liquid molecule. For the TIP4P / 2005 water model, the dipole moment vector is determined by the positions of two positively charged hydrogen atoms and one negatively charged virtual charge. The dipole moment vector of each water molecule in each frame is calculated by using the internal geometric relationship of the water molecule and the atomic charge parameters.

[0030] Optionally, the molecular orientation distribution can be obtained by analyzing the statistical distribution of molecular orientation within each thin layer based on the angle between the dipole moment vector and the surface normal vector of the polar solid. For water molecules within each thin layer, the angle θ between their dipole moment vector and the surface normal vector pointing into the liquid is calculated. The angle θ ranges from 0 to 180 degrees. The angle values ​​counted across all time frames within each thin layer are statistically analyzed to obtain a histogram or probability density function of the angle θ distribution. The density distribution and molecular orientation distribution are combined as characterization data of the liquid film microstructure. The density distribution is presented as a curve with the position along the surface normal direction as the abscissa and density as the ordinate. The molecular orientation distribution is presented as a family of curves indexed along the normal position, with the angle θ as the abscissa and probability density as the ordinate. These data together constitute a complete characterization of the liquid film microstructure. It can be understood that the process of calculating the density distribution by statistically analyzing the average change in the number of atoms or centroids within each thin layer over time follows the following formula:

[0031] Where: symbol Represents the average number density within the k-th thin layer, denoted by [symbol]. The symbol represents the ensemble average of the number of atomic or molecular centroids obtained statistically within the k-th thin layer over the entire sampling time. Represents the cross-sectional area of ​​the simulated system in the direction parallel to the surface, with the symbol... This represents the thickness of each thin layer.

[0032] In one embodiment of the present invention, the specific method for calculating the autocorrelation function of heat flow in the liquid film region using the Green-Kubo formula based on molecular coordinates and velocity trajectory data is as follows: Based on the molecular coordinates and velocity trajectory data, the total potential energy of the equilibrium molecular dynamics simulation system and the force on each atom at each time step are calculated using force field parameters. For the system composed of TIP4P / 2005 water molecules and an α-quartz surface, the force field parameters include bonding and non-bonding parameters. At each trajectory frame time with recorded coordinates, all bond stretching energies, angular bending energies, dihedral torsional energies, and van der Waals and electrostatic interaction potential energies are calculated based on the relative positions between atoms, and the total potential energy of the system is obtained by summing them. Simultaneously, the resultant force vector on each atom at this moment is calculated based on the negative gradient of the potential energy with respect to the atomic coordinates. Using the atomic velocities, total potential energy, and forces on atoms from the molecular coordinates and velocity trajectory data, the instantaneous heat flow vector at each time step is calculated using the heat flow definition formula. For a system containing N atoms, the instantaneous heat flow vector... The calculation is based on the microscopic heat flow expression in molecular dynamics, which includes the summation of the contributions of the kinetic energy term, the potential energy term, and the virial term to the heat flow.

[0033] In some embodiments, the heat flow component parallel to the polar solid surface is extracted from the instantaneous heat flow vector. The solid surface plane is defined as the xy plane, the surface normal direction is defined as the z direction, and the instantaneous heat flow vector... Two components in the xy plane and This refers to the heat flux component parallel to the surface. The autocorrelation function of the heat flux is obtained by calculating the autocorrelation average of the heat flux component parallel to the polar solid surface at different time intervals. For an equilibrium system, the x-component of the parallel heat flux is selected for calculation. The autocorrelation function is defined as the ensemble average of the product of the heat flux component at time t and the value at time t+τ. It is specifically calculated by averaging over all possible starting times t. Example data for calculating the autocorrelation function of the parallel heat flux components can be found in Table 1.

[0034] Table 1: Variation of heat flux autocorrelation function with time delay Time delay Heat flux autocorrelation function 0.00 105.32 0.05 12.47 0.10 -3.85 0.15 -1.22 0.20 0.18 0.25 0.05 Optionally, the specific method for obtaining the predicted thermal conductivity of the liquid film region by performing time integration on the heat flux autocorrelation function is as follows: Integrate the heat flux autocorrelation function from time zero to a sufficiently long correlation cutoff time. The integration operation is performed by summing the heat flux autocorrelation function values ​​at discrete time points using numerical integration methods such as the trapezoidal rule. Integration stops when the heat flux autocorrelation function decays to zero after the correlation cutoff time or fluctuates slightly around its mean. The correlation cutoff time is usually selected as the time point when the heat flux autocorrelation function value decays to within 5% of its initial absolute value or when the function value fluctuates randomly around zero. Substitute the integration result into the thermal conductivity calculation formula and combine it with the system's volume and temperature parameters to calculate the predicted thermal conductivity of the liquid film region in the direction parallel to the surface. The system volume parameter is the product of the average dimensions of the simulation box in the x, y, and z directions during the equilibrium simulation stage, and the temperature parameter is the average temperature of the equilibrium system.

[0035] It can be understood that the expression for the heat flux autocorrelation function is:

[0036] Where: symbol Represents the autocorrelation function of heat flux, which is a time delay. The function, symbol Represents the start time ensemble mean, sign Represents time The x-component of the heat flow vector parallel to the solid surface, sign... Represents time The x-component of the heat flux vector parallel to the solid surface. In some embodiments, when calculating the autocorrelation function of the heat flux, the autocorrelation functions of the x and y components of the parallel heat flux can be calculated separately and then averaged.

[0037] See Figure 4 This is a monitoring curve of the energy minimization process in a molecular dynamics simulation system. The exponential decay of potential energy is a typical characteristic of energy minimization, indicating that algorithms such as conjugate gradients have effectively found the low-energy conformation of the system, avoiding system collapse caused by an unreasonable initial configuration. The temperature stabilizes after initial fluctuations, proving that the temperature control algorithm successfully maintained the system at the target temperature, laying the foundation for subsequent equilibrium molecular dynamics simulations. The synchronous stability of potential energy and temperature indicates that the initial molecular dynamics simulation system has completed pre-equilibrium processing and can proceed to the next step of thermodynamic equilibrium simulation. By analyzing the convergence rate of potential energy and temperature, the optimal number of relaxation steps for energy minimization can be determined, avoiding the waste of computational resources caused by excessive relaxation or the instability of the system caused by insufficient relaxation.

[0038] In one embodiment of the present invention, the specific method for analyzing the influence mechanism of polar surface interaction on thermal properties by combining liquid film microstructure characterization data and thermal conductivity prediction values ​​is as follows: The density distribution curve and molecular orientation distribution curve in the liquid film microstructure characterization data are correlated with the predicted thermal conductivity values. The predicted thermal conductivity values ​​are derived from calculations of the thermal transport properties of the equilibrium molecular dynamics simulation system parallel to the surface direction. The density distribution curve provides information on the density variation of the atomic or molecular centroids along the surface normal direction, while the molecular orientation distribution curve provides information on the probability distribution of the angle between the polar molecular dipole moments and the surface normal. The position and intensity of the oscillation peaks appearing near the polar solid surface in the density distribution curve are identified. In the interface system of α-quartz surface and water, the density distribution curve exhibits obvious oscillation peaks at distances of 0.25 nm and 0.45 nm from the surface, with the peak intensity of the first oscillation peak reaching 1.8 times the bulk density. The preferential orientation angle of molecules in a specific region is shown in the molecular orientation distribution curve. In a thin layer within 0.3 nm of the α-quartz surface, the probability distribution of the angle θ between the dipole moment vector of water molecules and the surface normal shows a significant peak in the range of 70 to 110 degrees, indicating that the molecules tend to align parallel to the surface or at a large angle to the surface normal.

[0039] In some embodiments, the correlation between the position of the oscillation peak, the intensity of the oscillation peak, the preferred orientation angle and the predicted thermal conductivity is analyzed to infer the relationship between the interface layer structure and the heat transfer efficiency. For α-quartz surfaces with different surface hydroxyl densities, higher surface polarity corresponds to the first peak of the density oscillation peak being closer to the surface and having a higher peak intensity. At the same time, it corresponds to a more concentrated distribution of the preferred orientation angle of water molecules in the interface layer. These changes in microstructural characteristics are correlated with the increasing or decreasing trend of the predicted thermal conductivity values ​​parallel to the surface direction.

[0040] A series of molecular dynamics simulation systems with varying simulation parameters were constructed. These parameters included adjusting the partial charge of the hydroxyl hydrogen atoms on the α-quartz solid surface from +0.5564e to +0.4564e to simulate different surface charge densities, and adjusting the liquid film thickness from 5 nm to 3 nm and 7 nm. For each molecular dynamics simulation system with different parameters, the entire process from energy minimization to thermophysical property analysis was repeatedly executed to obtain a series of corresponding predicted thermal conductivity values ​​and liquid film microstructure characterization data. For each modified system, system construction, energy minimization, equilibrium simulation, trajectory sampling, microstructure calculation, and thermal conductivity calculation were performed independently, ultimately yielding multiple datasets containing surface charge density, liquid film thickness, peak interfacial layer density, average preferred orientation angle, and predicted thermal conductivity values. A quantitative relationship graph was established between simulation condition parameters, liquid film microstructure characterization data and predicted thermal conductivity. A curve was plotted with surface charge density or liquid film thickness as independent variables and predicted thermal conductivity as dependent variables. At the same time, the peak intensity of the first degradation oscillation peak or the average value of the molecular preferential orientation angle under the corresponding conditions were marked on the graph as auxiliary coordinates to intuitively show the interrelationship trend between multiple variables.

[0041] It is understandable that, in order to quantify the relationship between structural parameters and thermophysical properties, the intensity of the first peak of the density oscillation peak can be calculated. Compared with predicted thermal conductivity The linear correlation coefficient between them is expressed as follows:

[0042] Where: symbol Representative thermal conductivity prediction value With the intensity of the first density peak The Pearson correlation coefficient between them, sign Represents the total number of molecular dynamics simulation systems constructed under different conditions, with the symbol... Representing the Predicted thermal conductivity values ​​for the system, symbol Representing all The arithmetic mean of the predicted thermal conductivity values ​​of the system, sign Representing the The peak intensity of the first oscillation peak in the density distribution curve of the system, with sign... Representing all The arithmetic mean of the peak intensity of the first density peak of each system.

[0043] See Figure 5This is a correlation curve between liquid film thickness and predicted thermal conductivity, a core result of thermal property analysis in molecular dynamics simulation studies. When the liquid film is thin (<4 nm), due to the strong adsorption effect of the polar surface, the liquid film molecules are highly ordered, with significant density oscillations, resulting in enhanced phonon scattering and lower thermal conductivity. As the thickness increases, the liquid film region far from the surface gradually approaches the disordered structure of the bulk liquid, phonon scattering weakens, and thermal conductivity increases accordingly. This curve can provide a quantitative reference for micro / nano-scale thermal management: to improve the thermal conductivity of the liquid film, it can be achieved by optimizing the liquid film thickness; to control heat transfer efficiency, an appropriate liquid film thickness range can be selected based on this curve. This figure directly verifies the correlation mechanism of "liquid film microstructure-thickness-thermal conductivity," providing intuitive data support for revealing the influence of polar surfaces on the thermal properties of liquid films, and can serve as a visualization of core results in academic papers.

[0044] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for predicting the thermal properties of polar surface liquid films based on molecular dynamics, characterized in that, include: Construct an initial molecular dynamics simulation system that includes polar solid surfaces and polar liquid molecules; Temperature and pressure controls were applied to the initial molecular dynamics simulation system, and energy minimization was performed to generate a stable pre-equilibrium system. The pre-equilibrium system is run until it reaches thermodynamic equilibrium, forming an equilibrium molecular dynamics simulation system for sampling; Extract molecular coordinates and velocity trajectory data for a specified time period from the equilibrium molecular dynamics simulation system; Based on the molecular coordinates and velocity trajectory data, the density distribution and molecular orientation distribution of the liquid film region are calculated to generate liquid film microstructure characterization data. Based on the molecular coordinates and velocity trajectory data, the heat flux autocorrelation function of the liquid film region is calculated using the Green-Kubo formula. The heat flux autocorrelation function is integrated over time to obtain the predicted thermal conductivity of the liquid film region; By combining the microstructure characterization data of the liquid film with the predicted thermal conductivity, the influence mechanism of polar surface effects on thermal properties is analyzed.

2. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 1, characterized in that, Construct an initial molecular dynamics simulation system comprising a polar solid surface and polar liquid molecules, including: Based on the lattice parameters and surface polarity characteristics of the target polar solid surface, a solid surface cell model under periodic boundary conditions is established. Based on the force field parameter file of the target polar liquid molecule, determine the atomic charge distribution, bonding and non-bonding interaction parameters of the polar liquid molecule; Above the solid surface cell model, polar liquid molecules are arranged with a specified thickness and initial density to form a liquid region; A vacuum layer is provided above the liquid region to eliminate spurious interactions of periodic boundary conditions in the direction perpendicular to the surface; The solid surface cell model, the liquid region, and the vacuum layer are combined to form an initial molecular dynamics simulation system with three-dimensional periodicity.

3. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 1, characterized in that, Temperature and pressure controls are applied to the initial molecular dynamics simulation system, and energy minimization is performed to generate a stable pre-equilibrium system, including: A position constraint algorithm is applied to all atoms in the initial molecular dynamics simulation system to prevent system collapse caused by an unreasonable initial configuration; The conjugate gradient algorithm is invoked to search the potential energy surface of the initial molecular dynamics simulation system, find the local potential minimum point, and complete the energy minimization process; The algorithm releases the positional constraints on atoms in the liquid region while keeping atoms in the solid region fixed or applying flexible constraints. A constant temperature control algorithm is used to adjust the system temperature to the target temperature value, and a pressure control algorithm is used to adjust the pressure to the target pressure value; Under these conditions, a short-term molecular dynamics pre-run is performed to allow the system to initially relax and generate the stable pre-equilibrium system.

4. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 1, characterized in that, The pre-equilibrium system is run until thermodynamic equilibrium is reached, forming an equilibrium molecular dynamics simulation system for sampling, including: Long-term molecular dynamics simulations were performed on the pre-equilibrium system under isothermal and isobaric or isothermal and isochoric conditions. Real-time monitoring of the potential energy, temperature, pressure, and time-history changes of the system dimensions of the pre-equilibrium system; When the average values ​​of the time-history changes in potential energy, temperature, pressure, and system size remain stable and the fluctuation amplitude is less than a preset threshold, the system is determined to have reached thermodynamic equilibrium. The continuous simulated trajectory after reaching thermodynamic equilibrium is extracted, and the corresponding system state constitutes the equilibrium molecular dynamics simulation system used for sampling.

5. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 1, characterized in that, Extracting molecular coordinates and velocity trajectory data for a specified time period from the equilibrium molecular dynamics simulation system includes: A time length much longer than the molecular vibration period is set as the specified time period; During the molecular dynamics simulation, the three-dimensional coordinate vectors and three-dimensional velocity vectors of all atoms in the equilibrium molecular dynamics simulation system are recorded at fixed time intervals. The recorded three-dimensional coordinate vectors and three-dimensional velocity vectors are organized in chronological order to form the molecular coordinate and velocity trajectory data.

6. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 5, characterized in that, Based on the molecular coordinates and velocity trajectory data, the density distribution and molecular orientation distribution of the liquid film region are calculated to generate liquid film microstructure characterization data, including: The liquid film region is divided into a series of thin layers of equal thickness along a direction perpendicular to the surface of the polar solid; The density distribution is calculated by statistically analyzing the average change in the number of atoms or centroids within each thin layer over time. For polar liquid molecules, calculate the dipole moment vector for each molecule; Based on the angle between the dipole moment vector and the normal vector of the polar solid surface, the statistical distribution of molecular orientation within each thin layer is analyzed to obtain the molecular orientation distribution; The density distribution and the molecular orientation distribution are combined to form the microstructure characterization data of the liquid film.

7. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 5, characterized in that, Based on the molecular coordinates and velocity trajectory data, the heat flux autocorrelation function of the liquid film region is calculated using the Green-Kubo formula, including: Based on the molecular coordinates and velocity trajectory data, the total potential energy of the equilibrium molecular dynamics simulation system and the force on each atom at each time step are calculated using force field parameters. Using the atomic velocities, total potential energy, and forces acting on the atoms from the molecular coordinates and velocity trajectory data, the instantaneous heat flux vector at each time step is calculated using the heat flux definition formula; Extract the heat flow component parallel to the polar solid surface from the instantaneous heat flow vector; The autocorrelation average value of the heat flow component parallel to the polar solid surface at different time intervals is calculated to obtain the heat flow autocorrelation function.

8. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 7, characterized in that, The heat flux autocorrelation function is integrated over time to obtain the predicted thermal conductivity of the liquid film region, including: Integrate the heat flow autocorrelation function from time zero to a sufficiently long correlation cutoff time; Integration stops when the heat flow autocorrelation function decays to zero or fluctuates slightly around its mean after the correlation cutoff time. Substituting the integral result into the thermal conductivity calculation formula, and combining it with the system's volume and temperature parameters, the predicted thermal conductivity of the liquid film region in the direction parallel to the surface is calculated.

9. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 1, characterized in that, Combining the microstructure characterization data of the liquid film with the predicted thermal conductivity, the influence mechanism of polar surface interaction on thermal properties is analyzed, including: The density distribution curve and molecular orientation distribution curve in the microstructure characterization data of the liquid film are correlated with the predicted thermal conductivity value. Identify the location and intensity of the oscillation peaks that appear near the surface of the polar solid in the density distribution curve; Identify the preferred orientation angles of molecules in specific regions as shown in the molecular orientation distribution curve; By analyzing the correspondence between the oscillation peak position, the oscillation peak intensity, the preferred orientation angle, and the predicted thermal conductivity, the correlation between the interface layer structure and the heat transfer efficiency can be inferred.

10. The method for predicting the thermal properties of polar surface liquid films based on molecular dynamics according to claim 9, characterized in that, After analyzing the correspondence between the oscillation peak position, the oscillation peak intensity, the preferred orientation angle, and the predicted thermal conductivity value, the method further includes: By changing the simulation condition parameters, a series of molecular dynamics simulation systems with different condition parameters can be constructed; For each molecular dynamics simulation system with different condition parameters, the entire process from energy minimization to thermal property analysis is repeated to obtain a series of corresponding predicted thermal conductivity values ​​and liquid film microstructure characterization data. Establish a quantitative relationship graph between the simulation condition parameters, the liquid film microstructure characterization data, and the predicted thermal conductivity value.

Images