Construction method of molecular dynamics model for different bubble contents of ice and methane hydrate

Abstract

The invention belongs to the technical field of dynamic simulation, and discloses a construction method of a molecular dynamic model for different bubble contents of ice and methane hydrate, and the method comprises the following specific steps: 1, selecting a crystal structure of the ice and methane hydrate as a model basis for initial structure construction; by constructing a molecular dynamic model of an ice and methane hydrate coexistence system, the problem of extensive quantitative regulation and control of bubble content is solved, accurate association of microscopic molecule removal and macroscopic volume fraction is realized, force field parameter suitability of a cross-phase interface is optimized, interaction of a hydrogen bond network and molecules at an ice and methane hydrate interface is accurately described, and the method has the advantages of simple operation and low cost. The simulation precision of the molecular diffusion behavior and the structural stability of the interface area is improved, meanwhile, the bubble form diversity and the temperature and pressure coupling environment are comprehensively considered, a multi-dimensional experimental verification system is established, and the influences of different bubble contents on the thermodynamic property, the dynamic behavior and the mechanical property of the system can be systematically evaluated.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of dynamics simulation, and particularly relates to a method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents. BACKGROUND

[0002] Methane hydrate, also known as combustible ice, is a solid compound formed by wrapping a cage-like crystal structure of methane molecules with water molecules through hydrogen bonds, and its formation needs to meet the conditions of low temperature and high pressure, and is commonly found in deep-sea sediments or permafrost zones.

[0003] As an energy storage carrier and geological medium, the structural stability of ice and methane hydrate is closely related to the bubble content, and the construction of a molecular dynamics model with different bubble contents is of great significance to reveal the phase transition mechanism, mechanical properties and multiphase coupling behavior. At present, a lot of research has been carried out on the molecular dynamics model of single ice phase or methane hydrate phase, but there are still technical bottlenecks in the modeling method of different bubble contents in the coexistence system of the two. In the prior art, the quantitative control method of bubble content is relatively rough, and it is difficult to realize the accurate correlation between micro-molecular arrangement and macro-volume fraction; the description of the interaction between ice and methane hydrate across the phase interface is simplified, and the cross-system adaptability of the force field parameters is insufficient, which limits the simulation accuracy of the hydrogen bond network and molecular diffusion behavior at the interface; at the same time, the adaptability of the bubble shape and the simulation environment is not comprehensive enough in the model construction process, and the experimental verification system is also limited to single structure characterization, which is difficult to systematically evaluate the comprehensive influence of bubbles on the thermodynamic and dynamic properties of the system, so a more perfect modeling method is needed to improve the characterization ability of the model to the real complex system. SUMMARY

[0004] In view of the deficiencies of the prior art, the application provides a method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents, which has the advantages of precise multiphase modeling.

[0005] To achieve the above purpose, the application provides the following technical scheme: a method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents, comprising the following specific steps:

[0006] Step 1: initial structure construction

[0007] The crystal structure of ice and methane hydrate is selected as the basis of the model, and the most common hexagonal ice Ih structure in nature is adopted, and the cell parameters are a = 6.57 A, c = 5.43 A. Each unit cell contains 4 water molecules, and the water molecules form a hexagonal arrangement through hydrogen bonds; the structure I (sI) type of methane hydrate is selected, and the cell parameters are a = 17.30 A, c = 10.00 A. Each unit cell contains 46 water molecules and 8 methane molecules, the water molecules form two kinds of cage structures, pentagonal dodecahedron and tetrakaidecahedron, and the methane molecules fill in the cages; the initial unit cell is constructed using molecular modeling software to ensure that the atoms are correctly positioned and that the positions of the water molecules and methane molecules conform to the characteristics of the crystal structure;

[0008] Step two: bubble introduction and content control

[0009] The bubble content is expressed as a volume fraction φ, the calculation formula is φ = (V _bubble / V _total ) x 100%, where V _bubble is the bubble volume, V _total is the total volume of the unit cell, for sI type methane hydrate unit cell, V _total = a 3 ; assuming the bubble is spherical, the radius r = (3V _bubble / (4π))^(1 / 3), the bubble radius r needs to meet 2r < min(Lx, Ly, Lz), where L is the cell edge length, otherwise the Wigner-Seitz primitive cell is used to correct the periodic boundary effect, which needs to ensure that the bubble does not exceed the cell boundary; generate uniformly distributed random coordinates (x, y, z) in the unit cell through a random number generator, the range is within the cell size, check if the selected position is outside the existing bubble area to avoid overlap; remove the water molecules at this position to form a vacuum area, the gas phase methane molecules in the bubble can be selectively filled to simulate the actual environment, repeat this process until the target bubble number is reached;

[0010] Step three: force field selection and parameter setting

[0011] Ice uses the TIP4P / Ice force field, which can accurately describe the hydrogen bond interaction between water molecules, the parameters include the Lennard-Jones potential parameter of the oxygen atom ε = 78.0 kJ / mol and the charge distribution of the hydrogen atom; methane uses the OPLS-AA force field, the carbon-oxygen interaction parameters are obtained through literature data or quantum chemical calculations to ensure that the van der Waals force and electrostatic interaction between methane molecules and water molecules conform to the actual situation; for the ice-hydrate interface region, the Lorentz-Berthelot mixing rule is used to optimize the LJ parameters of different atoms, and the interface charge distribution is calibrated through quantum chemical calculations;

[0012] Step four: energy minimization

[0013] The conjugate gradient method is used to minimize the energy of the initial structure to eliminate unreasonable overlaps and high stress areas between atoms, and the convergence criterion is set to be less than In the energy minimization process, the cell size is fixed, only the atomic coordinates are optimized, the iteration number is 10,000 steps or until the total potential energy of the system no longer decreases significantly, if a local energy minimum appears, the steepest descent method can be switched to continue optimization to ensure that the system reaches the global energy minimum;

[0014] Step five: NVT equilibration simulation

[0015] After energy minimization, the structure is placed in a canonical ensemble (NVT) for equilibration simulation, the control temperature is 273K, and the time step is set to 1fs; the Andersen thermostat is used, with a collision frequency of 1ps-1 to maintain stable temperature, the simulation time is 100ps, during which the potential energy, temperature and density of the system are monitored to ensure that they tend to be stable, if the temperature fluctuation is large, the simulation time can be extended or the temperature control parameters can be adjusted; after equilibration, the atomic coordinates and velocities of the system are saved as the initial state for subsequent simulation;

[0016] Step six: NPT equilibration simulation

[0017] The structure after NVT equilibration is converted to an isothermal-isobaric ensemble (NPT), while controlling the temperature and pressure, the temperature is maintained at 273K, and the pressure is set to 1atm; the Parrinello-Rahman pressure control method is used, with a pressure coupling constant of 1ps, allowing the cell size to be freely adjusted in three directions; the simulation time is more than 3 times the diffusion characteristic time of the system, during which the pressure, volume and density of the system are monitored to ensure that it reaches a stable state; during the equilibration process, the cell volume may change due to different bubble contents, so the final cell parameters need to be recorded for production simulation;

[0018] Step seven: production simulation and data collection

[0019] Based on the structure after NPT equilibration, production simulation is carried out, the time length is 1-2ns, the time step is kept at 1fs, and the Andersen thermostat and Parrinello-Rahman pressure control method are continued to be used to ensure the stability of the system state, during the simulation process, the atomic coordinates and velocities are saved once every 100 steps for subsequent analysis; at the same time, the energy, temperature, pressure, density and other parameters of the system are monitored in real time to generate a log file, if abnormal parameters are found, the initial structure needs to be checked or the simulation parameters need to be adjusted;

[0020] Step eight: data analysis and model validation

[0021] The production simulation trajectory is analyzed, and the RDF (radial distribution function) is calculated to study the intermolecular interaction; the O-O RDF of ice should have a main peak at , corresponding to the hydrogen bond distance; the C-O RDF of methane hydrate should be at The peaks on the left and right indicate that the methane molecules are located in the cage. In addition, the MSD (mean square displacement) and diffusion coefficient are calculated to analyze the motion characteristics of water molecules and methane molecules. The simulation results are compared with experimental data to ensure the accuracy of the model.

[0022] Step nine: Construction and comparison of different bubble content models

[0023] Repeat steps two to eight to construct models with different bubble contents. The different bubble content models should be constructed from the same initial structure after energy minimization to ensure the uniqueness of the variables. During the bubble introduction process, the number of water molecules removed and the volume of methane molecules filled are adjusted to accurately control the bubble content. The same simulation process and data analysis are performed on each model to compare the structural stability of ice and methane hydrate, the molecular dynamics behavior, and the changes in mechanical properties under different bubble contents. The influence mechanism of bubble content on model performance is analyzed.

[0024] Step ten: Result optimization and model correction

[0025] Based on the data analysis results, the model is optimized and corrected. If the structure is unstable or deviates significantly from the experimental data, the force field parameters, initial structure, or simulation settings need to be checked. If the RDF peak position is shifted, the hydrogen bond parameters of the force field need to be adjusted. If the diffusion coefficient is abnormal, the equilibrium simulation time needs to be extended or the temperature control method needs to be adjusted. Through repeated iterative optimization, the model can accurately reflect the molecular dynamics characteristics of ice and methane hydrate under different bubble contents.

[0026] Step eleven: Visualization and document organization

[0027] Use molecular visualization software to visualize the simulation results, generate cell structure, bubble distribution, molecular motion trajectory images and animations, and label the Voronoi grid density in the bubble distribution image. The frame rate of the molecular motion trajectory animation is ≥1 ps / frame. Organize the simulation data, including energy change curves, RDF graphs, MSD curves, etc., to form a detailed technical report. The report should include model construction methods, simulation parameter settings, data analysis results, and conclusions to ensure the repeatability of the method and the credibility of the results.

[0028] Preferably, the cell parameters a and c of ice Ih structure in step one need to be accurate to The methane hydrate cell parameter a needs to be accurate to Ensure that the deviation of the initial structure from the real crystal structure is within 1%.

[0029] Preferably, the bubble introduction process in step two is only applicable to the ideal case of low content and uniform interfacial tension, assuming spherical bubbles. When the bubble content is high or coalescence occurs, a non-spherical bubble model needs to be used and the volume calculation method needs to be corrected. At this time, the bubble radius formula needs to introduce a shape factor to ensure that the calculation error of bubble volume ratio is not more than 0.5%.

[0030] Preferably, the convergence criterion in step four is set to be less than Only applicable to the case that the initial structure has no serious overlap, if the initial structure has an atomic distance less than the sum of the van der Waals radii (such as O-H distance ), rigid constraint optimization is needed first, and then the degrees of freedom are gradually released to prevent the energy minimization from falling into a local minimum value.

[0031] Preferably, the temperature in step five is set to 273K, and the phase state of the system needs to be determined. If the model contains liquid water regions, the temperature needs to be adjusted to 298K and the TIP3P force field is used to describe liquid water molecules. At the same time, the collision frequency needs to be dynamically adjusted according to the mass density of the system to ensure that the temperature control accuracy is within ±1K. When the system contains liquid water, the hydrogen atom charge of the TIP3P force field is adjusted to +0.417e, and the oxygen atom charge is adjusted to -0.834e.

[0032] Preferably, the pressure in step six is set to 1atm, which is only applicable to the surface or shallow layer environment. If the deep-sea high-pressure environment is simulated, the Berendsen pressure control method needs to be used and the pressure coupling constant needs to be shortened to 0.1ps, while allowing the anisotropic adjustment of the cell size in three directions.

[0033] Preferably, in step seven, the atomic coordinates and velocities are saved once every 100 steps, which needs to be adjusted according to the simulation target. When studying fast dynamics, the interval should be shortened to 10 steps; when studying long-term structural evolution, it can be extended to 500 steps.

[0034] Preferably, in step eight, the RDF (radial distribution function) calculation needs to include at least 1000 time frames to ensure statistical significance, and the cutoff radius needs to be greater than the first solvation layer. When comparing experimental data, the experimental conditions need to be clear.

[0035] Preferably, in step ten, the result optimization and model correction, when the deviation between the simulated value and the experimental value of the diffusion coefficient is more than 20%, the methane molecule charge distribution is calibrated first (±0.05e), then the ε value of the LJ potential is adjusted (±5%); Finally, by extending the simulation to a diffusion coefficient error of less than 10%, avoid the deviation of the optimization direction caused by the wrong description of the basic interaction.

[0036] Compared with the prior art, the beneficial effects of the present application are as follows:

[0037] By constructing a molecular dynamics model of the coexistence system of ice and methane hydrate, the problem of extensive quantitative control of bubble content is solved, the precise correlation between microscopic molecular removal and macroscopic volume fraction is achieved, and the adaptability of force field parameters across phase interfaces is optimized. The hydrogen bond network and molecular interactions at the interface between ice and methane hydrate are accurately described, and the simulation accuracy of molecular diffusion behavior and structural stability in the interface region is improved. At the same time, considering the diversity of bubble morphology and the temperature-pressure coupling environment, a multi-dimensional experimental verification system is established, which can systematically evaluate the influence of different bubble contents on the thermodynamic properties, kinetic behavior and mechanical characteristics of the system, providing a reliable model support for the research on issues such as bubble migration prediction in natural gas hydrate extraction and the study of deep-sea pipeline blockage mechanisms. DETAILED DESCRIPTION

[0038] Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative work shall fall within the scope of protection of the present invention.

[0039] The embodiment of the present invention provides a method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents, and the specific steps are as follows:

[0040] Step 1: Initial structure construction

[0041] The crystal structures of ice and methane hydrate are selected as the basis of the model. The most common hexagonal ice Ih structure in nature is adopted, and its unit cell parameters are Each unit cell contains 4 water molecules, which are arranged in a hexagonal pattern through hydrogen bonds. The methane hydrate selects structure I (sI) type, and its unit cell parameters are Each unit cell contains 46 water molecules and 8 methane molecules. The water molecules form two cage structures, a pentagonal dodecahedron and a tetradecahedron, and the methane molecules fill the cages. Molecular modeling software was used to construct the initial unit cell to ensure that the atomic coordinates were accurate and the positions of the water and methane molecules conformed to the crystal structure characteristics.

[0042] Step 2: Bubble introduction and content control

[0043] The bubble content is expressed as volume fraction φ, and the calculation formula is φ=(V _bubble / V _total )×100%, where V _bubble is the bubble volume, V _total is the total volume of the unit cell. For the sI type methane hydrate unit cell, V _total =a 3 ; Assume that the bubble is spherical, with radius r=(3V _bubble(4π)^(1 / 3), the bubble radius r needs to satisfy 2r < min(Lx, Ly, Lz), where L is the cell edge length, otherwise the Wigner-Seitz cell is used to correct the periodic boundary effect, and it is necessary to ensure that the bubble does not exceed the cell boundary; a random number generator is used to generate uniformly distributed random coordinates (x, y, z) in the cell, within the cell size, and check whether the selected position is outside the existing bubble area to avoid overlap; remove the water molecules at this position to form a vacuum area, and optionally fill the gas-phase methane molecules in the bubble to simulate the actual environment, and repeat this process until the target number of bubbles is reached;

[0044] Step three: force field selection and parameter setting

[0045] Ice uses the TIP4P / Ice force field, which can accurately describe the hydrogen bond interaction between water molecules, and the parameters include the Lennard-Jones potential parameter of the oxygen atom ε = 78.0 kJ / mol, and the charge distribution of the hydrogen atom; methane uses the OPLS-AA force field, and the carbon-oxygen interaction parameters are obtained through literature data or quantum chemical calculation to ensure that the van der Waals force and electrostatic interaction between methane molecules and water molecules are consistent with the actual situation; for the ice-water interface region, the Lorentz-Berthelot mixing rule is used to optimize the LJ parameters of different atoms, and the interface charge distribution is calibrated through quantum chemical calculation;

[0046] Step four: energy minimization

[0047] The conjugate gradient method is used to minimize the energy of the initial structure to eliminate unreasonable overlaps and high stress areas between atoms, and the convergence criterion is set to be less than During the energy minimization process, the cell size is fixed and only the atomic coordinates are optimized, with an iteration number of 10,000 steps or until the total potential energy of the system no longer decreases significantly, and if a local energy minimum value appears, the steepest descent method can be switched to continue optimization to ensure that the system reaches the global energy minimum value;

[0048] Step five: NVT equilibrium simulation

[0049] The structure after energy minimization is placed in the canonical ensemble (NVT) for equilibrium simulation, with a temperature of 273 K and a time step of 1 fs; the Andersen thermostat is used with a collision frequency of 1 ps-1 to maintain a stable temperature, and the simulation time is 100 ps, during which the potential energy, temperature and density of the system are monitored to ensure that they tend to be stable, and if the temperature fluctuation is large, the simulation time can be extended or the temperature control parameters can be adjusted; after the equilibrium is completed, the atomic coordinates and velocities of the system are saved as the initial state for subsequent simulation;

[0050] Step six: NPT equilibrium simulation

[0051] Convert the NVT equilibrated structure to the isothermal-isobaric ensemble (NPT) while controlling temperature and pressure, maintain temperature at 273 K and pressure at 1 atm; use Parrinello-Rahman pressure coupling method with a coupling constant of 1 ps to allow the cell size to adjust freely in three directions; simulate for more than 3 times the characteristic diffusion time of the system, monitor the pressure, volume, and density changes during the simulation to ensure the system reaches a stable state; during the equilibration process, the cell volume may change due to different bubble contents, record the final cell parameters for production simulation;

[0052] Step Seven: Production Simulation and Data Collection

[0053] Perform production simulation based on the NPT equilibrated structure, duration is 1-2 ns, time step remains 1 fs, continue using Andersen thermostat and Parrinello-Rahman pressure coupling method to ensure system stability, save atomic coordinates and velocities every 100 steps during simulation for subsequent analysis; at the same time, monitor system energy, temperature, pressure, density and other parameters in real time, generate log files, if abnormal parameters are found, need to recheck the initial structure or adjust simulation parameters;

[0054] Step Eight: Data Analysis and Model Validation

[0055] Analyze the production simulation trajectory, calculate RDF (Radial Distribution Function) to study the intermolecular interactions; the O-O RDF of ice should have a main peak at , corresponding to the hydrogen bond distance; the C-O RDF of methane hydrate should have peaks around , indicating that methane molecules are in the cage; in addition, calculate MSD (Mean Square Displacement) and diffusion coefficient to analyze the motion characteristics of water molecules and methane molecules; compare the simulation results with experimental data to ensure the accuracy of the model;

[0056] Step Nine: Construction and Comparison of Models with Different Bubble Contents

[0057] Repeat steps two to eight to construct models with different bubble contents, different bubble content models need to start from the same energy minimized initial structure to ensure the uniqueness of the variable, during the bubble introduction process, adjust the number of removed water molecules and the volume of filled methane molecules to achieve accurate control of the bubble content; perform the same simulation process and data analysis for each model, compare the structural stability, molecular dynamics behavior and mechanical properties of ice and methane hydrate under different bubble contents, analyze the influence mechanism of bubble content on model performance;

[0058] Step Ten: Result Optimization and Model Correction

[0059] According to the data analysis results, the model is optimized and corrected. If the structure is unstable or deviates greatly from the experimental data, the force field parameters, initial structure or simulation settings need to be checked. If the RDF peak position deviates, the hydrogen bond parameters of the force field need to be adjusted. If the diffusion coefficient is abnormal, the equilibrium simulation time needs to be extended or the temperature control method needs to be adjusted. Through repeated iteration optimization, it is ensured that the model can accurately reflect the molecular dynamics characteristics of ice and methane hydrate under different bubble contents.

[0060] Step eleven: visualization and document organization

[0061] Visualize the simulation results using molecular visualization software to generate images and animations of cell structure, bubble distribution, molecular motion trajectory, etc. Bubble distribution images need to be labeled with Voronoi grid density, and molecular motion trajectory animation frame rate ≥ 1 ps / frame. Organize simulation data, including energy change curve, RDF graph, MSD curve, etc. to form a detailed technical report. The report should include model construction method, simulation parameter setting, data analysis result and conclusion, to ensure the repeatability of the method and the credibility of the result.

[0062] Based on the ice Ih and sI type methane hydrate cell, the initial structure is constructed. After calculating and removing water molecules according to the target bubble volume fraction, methane molecules are inserted. TIP4P / Ice and OPLS-AA force field optimization parameters are selected. Energy minimization, NVT and NPT ensemble equilibrium simulation are carried out in turn. Monitor parameters and save data during production simulation. Analyze molecular interaction by calculating radial distribution function. Compare experimental data to verify the model. Repeat the process to build models with different bubble contents and compare the differences. According to the results, optimize and correct the visualization data and organize the report to form a systematic modeling method.

[0063] In step one, the cell parameters a and c of ice Ih structure need to be accurate to The methane hydrate cell parameter a needs to be accurate to Ensure that the deviation of the initial structure from the real crystal structure is within 1%.

[0064] By strictly limiting the accuracy of the cell parameters, it is ensured that the initial model is highly consistent with the atomic arrangement, bond length and bond angle of the real crystal structure, avoiding the deviation of the subsequent simulation results from the actual situation due to the distortion of the basic structure.

[0065] In step two, the bubble introduction process assumes that the bubble is spherical and only applies to ideal situations with low content and uniform interfacial tension. When the bubble content is high or coalescence occurs, a non-spherical bubble model needs to be used and the volume calculation method needs to be corrected. The bubble radius formula needs to introduce a shape factor to ensure that the error of bubble volume ratio calculation is not more than 0.5%.

[0066] The non-spherical bubble model and shape factor are introduced to address the volume calculation deviation caused by the traditional spherical assumption, reduce the error in calculating the bubble volume ratio, ensure that the geometric parameters of the model match the actual system in high-pressure and high-gas scenarios, and improve the reliability of multiphase flow simulation.

[0067] In step four, the convergence criterion is set as the interatomic force being less than Only applicable to initial structures without severe overlap, if the initial structure has an atomic spacing less than the sum of the van der Waals radii (such as O-H spacing ), a rigid constraint optimization is required first, and then the degrees of freedom are gradually released to prevent the energy minimization from falling into a local minimum.

[0068] For initial structures with severe atomic overlap, a rigid constraint optimization is performed first, and then the degrees of freedom are released to effectively avoid the energy minimization from falling into a local minimum, ensuring that the system reaches a global energy minimum state.

[0069] In step five, the temperature is set to 273K, and the phase state of the system needs to be determined. If the model contains liquid water regions, the temperature needs to be adjusted to 298K and the TIP3P force field is used to describe the liquid water molecules. The collision frequency needs to be dynamically adjusted according to the mass density of the system to ensure that the temperature control accuracy is within ±1K. When the system contains liquid water, the hydrogen atom charge of the TIP3P force field is adjusted to +0.417e, and the oxygen atom charge is adjusted to -0.834e.

[0070] The temperature and force field are dynamically adjusted according to the phase state of the system, and the collision frequency is optimized to solve the molecular description contradiction in the liquid water region at the edge of the bubble, improving the temperature control accuracy from ±5K to ±1K, and avoiding the virtual high molecular diffusion coefficient caused by the mismatch between the force field and the phase state.

[0071] In step six, the pressure is set to 1 atm, which is only applicable to surface or shallow layer environments. If deep-sea high-pressure environments are simulated, the Berendsen pressure control method is used, and the pressure coupling constant is shortened to 0.1 ps, while allowing anisotropic adjustment of the cell size in three directions.

[0072] Adjustments are made for deep-sea high-pressure environments to avoid lattice distortion caused by improper pressure control, ensuring the credibility of the crystal symmetry and mechanical property simulation of the model under high pressure.

[0073] In step seven, atomic coordinates and velocities are saved every 100 steps, which needs to be adjusted according to the simulation target. For fast kinetic processes, the interval should be shortened to 10 steps; for long-term structural evolution, it can be extended to 500 steps.

[0074] For fast processes, high-frequency sampling is used to capture transient behaviors that may be missed by traditional fixed intervals, improving the effectiveness of data analysis.

[0075] Wherein, in step eight, RDF (Radial Distribution Function) calculation needs to contain at least 1000 time frames to ensure statistical significance, and the cut-off radius needs to be greater than the first solvation shell. When comparing experimental data, the experimental conditions need to be clear.

[0076] The RDF calculation is required to contain ≥1000 time frames and set a reasonable cut-off radius, so that the statistical error is greatly reduced, and when comparing with experimental data, the conditions are matched to avoid model misjudgment caused by insufficient samples or wrong conditions, and to ensure the scientific rigor of structure verification and dynamics analysis.

[0077] Wherein, in step ten, the results are optimized and the model is corrected. When the deviation between the simulation value and the experimental value of the diffusion coefficient is more than 20%, first, the charge distribution of the methane molecule is calibrated (±0.05e), and then the ε value of the LJ potential is adjusted (±5%). Finally, by extending the simulation to a diffusion coefficient error <10%, the deviation of the optimization direction caused by the error in the description of the basic interaction is avoided.

[0078] The correction path is provided when the deviation of the diffusion coefficient is more than 20%, so as to improve the parameter optimization efficiency, ensure that each round of optimization is based on clear physical mechanisms, and improve the pertinence and effectiveness of model correction.

[0079] It should be noted that, in this article, relational terms such as first and second are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "contain" or any other variant thereof are intended to cover non-exclusive inclusion, so that the process, method, article or equipment including a series of elements not only includes those elements, but also includes other elements not explicitly listed or inherent to such process, method, article or equipment.

[0080] Although the embodiments of the present application have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and variations can be made to the embodiments without departing from the principles and spirit of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents, characterized in that: The specific steps are as follows: Step 1: Initial structure construction The crystal structures of ice and methane hydrate are selected as the basis of the model. The most common hexagonal ice Ih structure in nature is adopted, and its unit cell parameters are Each unit cell contains 4 water molecules, which are arranged in a hexagonal pattern through hydrogen bonds. The methane hydrate selects structure I (sI) type, and its unit cell parameters are Each unit cell contains 46 water molecules and 8 methane molecules. The water molecules form two cage structures, a pentagonal dodecahedron and a tetradecahedron, and the methane molecules fill the cages. Molecular modeling software was used to construct the initial unit cell to ensure that the atomic coordinates were accurate and the positions of the water and methane molecules conformed to the crystal structure characteristics. Step 2: Bubble introduction and content control The bubble content is expressed as a volume fraction φ, and the calculation formula is φ = (V _bubble / V _total ) × 100%, where V _bubble is the bubble volume, V _total is the total volume of the unit cell. For the sI-type methane hydrate unit cell, V _total = a 3 ; assuming the bubble is spherical, the radius r = (3V _bubble / (4π))^(1 / 3). The bubble radius r needs to satisfy 2r < min(Lx, Ly, Lz), where L is the side length of the unit cell. Otherwise, the Wigner-Seitz primitive cell is used to correct the periodic boundary effect, and it is necessary to ensure that the bubble does not exceed the unit cell boundary; random coordinates (x, y, z) uniformly distributed within the unit cell are generated by a random number generator within the range of the unit cell size, and it is checked whether the selected position is outside the existing bubble region to avoid overlap; The water molecules at that location are removed to form a vacuum area, and the bubbles can be selectively filled with gaseous methane molecules to simulate the actual environment. This process is repeated until the target number of bubbles is reached; Step 3: Force field selection and parameter setting Ice uses the TIP4P / Ice force field, which accurately describes the hydrogen bonding interactions between water molecules. The parameters include the Lennard-Jones potential parameters of oxygen atoms. ε = 78.0 kJ / mol, and the charge distribution of hydrogen atoms; methane uses the OPLS-AA force field, and the carbon-oxygen interaction parameters are obtained from literature data or quantum chemical calculations to ensure that the van der Waals and electrostatic interactions between methane and water molecules are consistent with reality; for the ice-hydrate interface region, the Lorentz-Berthelot mixing rule is used to optimize the LJ parameters of heterogeneous atoms, and the interface charge distribution is calibrated through quantum chemical calculations; Step 4: Energy Minimization The conjugate gradient method is used to minimize the energy of the initial structure, eliminate unreasonable overlaps between atoms and high stress areas, and set the convergence standard to the interatomic force less than 0.01kJ / During the energy minimization process, the unit cell size is fixed and only the atomic coordinates are optimized. The number of iterations is 10,000 steps or until the total potential energy of the system no longer decreases significantly. If a local energy minimum occurs, the optimization can be switched to the steepest descent method to ensure that the system reaches the global energy minimum. Step 5: NVT balance simulation The energy-minimized structure was placed in the canonical ensemble (NVT) for equilibrium simulation, with the temperature controlled at 273 K and the time step set to 1 fs; an Andersen thermostat was used with a collision frequency of 1 ps. -1 To maintain temperature stability, the simulation time is 100 ps. During this period, the potential energy, temperature and density changes of the system are monitored to ensure that they tend to be stable. If the temperature fluctuates greatly, the simulation time can be extended or the temperature control parameters can be adjusted; After the equilibrium is completed, the atomic coordinates and velocities of the system are saved as the initial state of the subsequent simulation; Step 6: NPT equilibrium simulation The NVT-equilibrated structure is converted to an isothermal isobaric ensemble (NPT) while controlling both temperature and pressure, maintaining the temperature at 273 K and the pressure at 1 atm. The Parrinel lo-Rahman pressure-controlled method is used with a pressure coupling constant of 1 ps, allowing the unit cell size to be freely adjusted in three directions. The simulation duration is at least three times the system's diffusion characteristic time, during which changes in the system's pressure, volume, and density are monitored to ensure that it reaches a stable state. During the equilibration process, the unit cell volume may vary due to different bubble content, and the final unit cell parameters must be recorded for production simulations. Step 7: Production simulation and data collection A production simulation is performed based on the structure after NPT equilibrium, with a duration of 1-2 ns and a time step of 1 fs. The Andersen temperature controller and the Parrinel lo-Rahman pressure control method are continued to be used to ensure the stability of the system state. During the simulation, the atomic coordinates and velocities are saved every 100 steps for subsequent analysis. At the same time, the system parameters such as energy, temperature, pressure, and density are monitored in real time, and log files are generated. If any abnormal parameters are found, the initial structure needs to be rechecked or the simulation parameters need to be adjusted. Step 8: Data Analysis and Model Validation The production simulation trajectory is analyzed and RDF (Radial Distribution Function) is calculated to study the interaction between molecules; the OO RDF of ice should be The main peak appears at , which corresponds to the hydrogen bond distance; the CO RDF of methane hydrate should be Peaks appear on the left and right, indicating that the methane molecules are located in the cage. In addition, the MSD (mean square displacement) and diffusion coefficient are calculated to analyze the motion characteristics of water and methane molecules. The simulation results are compared with experimental data to ensure the accuracy of the model. Step 9: Construction and comparison of models with different bubble contents Repeat steps 2 to 8 to construct models with different bubble contents. Models with different bubble contents must be constructed starting from the same initial structure after energy minimization to ensure variable uniqueness. During the bubble introduction process, the bubble content is precisely controlled by adjusting the number of water molecules removed and the volume of methane molecules filled. The same simulation process and data analysis are performed on each model to compare the structural stability, molecular dynamics behavior, and mechanical properties of ice and methane hydrates at different bubble contents, and analyze the mechanism by which bubble content affects model performance. Step 10: Result optimization and model modification Based on the data analysis results, the model is optimized and modified. If the structure is found to be unstable or deviates significantly from the experimental data, the force field parameters, initial structure, or simulation settings need to be checked. If the RDF peak position is shifted, the hydrogen bond parameters of the force field need to be adjusted. If the diffusion coefficient is abnormal, the equilibrium simulation time needs to be extended or the temperature control method needs to be adjusted. Through repeated iterative optimization, it is ensured that the model can accurately reflect the molecular dynamics characteristics of ice and methane hydrate at different bubble contents. Step 11: Visualization and Documentation Use molecular visualization software to visualize the simulation results and generate images and animations of unit cell structure, bubble distribution, and molecular motion trajectory. The bubble distribution image must be annotated with Voronoi grid density, and the frame rate of the molecular motion trajectory animation must be ≥1ps / frame. Organize the simulation data, including energy change curves, RDF graphs, MSD curves, etc., to form a detailed technical report. The report should include the model construction method, simulation parameter settings, data analysis results, and conclusions to ensure the repeatability of the method and the credibility of the results.

2. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: The unit cell parameters a and c of the ice Ih structure described in step 1 must be accurate to The methane hydrate unit cell parameter a must be accurate to The deviation between the initial structure and the true crystal structure was ensured to be within 1%.

3. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: The bubble introduction process described in step 2 assumes that the bubbles are spherical, which is only applicable to the ideal situation with low bubble content and uniform interfacial tension. When the bubble content is high or coalescence occurs, a non-spherical bubble model must be used and the volume calculation method must be modified. In this case, the bubble radius formula needs to include a shape factor to ensure that the calculation error of the bubble volume fraction does not exceed 0.5%.

4. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: The convergence criterion in step 4 is set to an interatomic force less than 0.01 kJ / It is only applicable when the initial structure has no serious overlap. If the initial structure has interatomic distances less than the sum of the van der Waals radii (such as OH distance < ), it is necessary to first perform rigid constraint optimization and then gradually release the degrees of freedom to prevent energy minimization from falling into a local minimum.

5. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: When the temperature is set to 273K in step 5, the phase state of the system needs to be clarified. If the model contains liquid water areas, the temperature needs to be adjusted to 298K and the TIP3P force field needs to be used to describe the liquid water molecules. At the same time, the collision frequency needs to be dynamically adjusted according to the mass density of the system to ensure that the temperature control accuracy is within ±1K. When the system contains liquid water, the hydrogen atom charge of the TIP3P force field is adjusted to +0.417e and the oxygen atom charge is adjusted to -0.834e.

6. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: The pressure of 1 atm in step 6 is only applicable to surface or shallow surface environments. If the deep-sea high-pressure environment is to be simulated, the Berendsen pressure-controlled method should be adopted and the pressure coupling constant should be shortened to 0.1 ps, while allowing anisotropic adjustment of the unit cell size in three directions.

7. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: As described in step 7, atomic coordinates and velocities are saved every 100 steps. This interval needs to be adjusted according to the simulation goal. When studying fast dynamic processes, the interval should be shortened to 10 steps; when studying long-term structural evolution, it can be extended to 500 steps.

8. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: The RDF (Radial Distribution Function) calculation described in step 8 must include at least 1000 time frames to ensure statistical significance, and the cutoff radius must be larger than the first solvation shell. The experimental conditions must be clearly defined when comparing experimental data.

9. The method for constructing a molecular dynamics model of ice and methane hydrate with different bubble contents according to claim 1, characterized in that: In step 10, the results were optimized and the model was modified. When the simulated diffusion coefficient deviated from the experimental value by more than 20%, the charge distribution of the methane molecules was first calibrated (±0.05e), followed by adjusting the ε value of the LJ potential (±5%). Finally, the simulation was extended until the diffusion coefficient error was <10% to avoid the optimization direction deviation caused by the incorrect description of the basic interaction.

Images