A method for determining the wettability of rock surface in cap rock of aquifer hydrogen storage reservoir

Abstract

The present invention discloses a method for measuring the surface wettability of the cap rock of an aquifer hydrogen storage reservoir, and belongs to the technical field of wettability measurement. The method includes: establishing a semi-cylindrical water molecule aggregate model on the surface of the cap rock mineral model; dividing the simulation area and boundary conditions; selecting a force field function to describe the interaction forces within and between molecules; structurally optimizing the contact surface model of the hydrogen-water-cap layer three-phase coexistence; determining relevant simulation parameters and ensembles; and performing molecular dynamics simulation on the hydrogen-water-cap layer wettability molecular model to obtain wettability simulation results. The present invention is applied to the measurement of the surface wettability of the cap rock of an aquifer hydrogen storage reservoir, and solves the problem that the existing wettability measurement method cannot accurately reflect the surface wettability of the cap rock of an aquifer hydrogen storage reservoir and has low calculation efficiency. The method has the characteristics of being able to calculate the wetting angle of water molecules in a hydrogen environment on the mineral surface under different mineral materials and having high calculation accuracy.

Description

Technical Field

[0001] The invention belongs to the technical field of wettability determination, and in particular relates to a method for determining the wettability of a rock surface in a cap layer of an aquifer hydrogen storage reservoir. Background Art

[0002] The wettability of the cap rock of an aquifer hydrogen storage reservoir refers to the tendency of formation water in hydrogen to expand or adhere to the surface of the cap rock minerals. It is one of the important properties of the surface characteristics of the cap rock of a hydrogen storage reservoir. Wettability has an important influence on the flow law of hydrogen and water two-phases and the sealing of the cap rock, and it plays an important role in the site selection and operation of aquifer hydrogen storage reservoirs.

[0003] Traditional wettability characterization methods are mainly based on macroscopic qualitative analysis, and there are few quantitative studies on the wettability of caprocks at meso-micro scales, especially wettability in hydrogen. Chinese patent CN112382345A discloses a shale oil reservoir contact angle calculation method and system, which respectively constructs a shale oil reservoir rock model, a shale oil molecule model and a water molecule model, builds a microscopic oil-water-rock wettability molecular model, performs molecular dynamics simulation on the oil-water-rock wettability molecular model, extracts water molecules in the oil-water-rock wettability molecular model from the simulation results, clusters the water molecules to form water clusters, and then digitizes the water clusters to determine the contact angle of the digitized water clusters on the rock surface of the shale oil reservoir at each simulation moment, overcomes the defects of the prior art that the two-dimensional image in the angle measurement method cannot well reflect the morphological characteristics of the water clusters and the subjective factors interfere with the artificial measurement angle, and obtains the molecular coordinates, extracts the parameters required for calculating the contact angle for calculation, and greatly improves the calculation accuracy of the contact angle.

[0004] However, the analysis object of the above method is shale oil reservoir. In the calculation of the wetting angle of shale oil reservoir, the main focus is on the interaction between oil molecules and the rock surface and water. The wetting behavior of oil is relatively simple because its interfacial energy and the interaction between water and rock are relatively clear. However, in the analysis of the cap rock of the aquifer hydrogen reservoir, the wetting angle not only involves the interaction between water and the rock surface, but also needs to consider the competitive adsorption and interface effect between hydrogen and water. Since hydrogen molecules have a smaller volume in the gaseous state, their interaction with water is more complex than that of common hydrocarbon molecules. Compared with hydrocarbon molecules, the interaction between hydrogen and water has significant differences at the physical and chemical level, resulting in more complex wetting behavior. The combined action of water and hydrogen may affect the adsorption behavior of water on the rock wall, especially in the cap rock of the hydrogen reservoir. The presence of hydrogen may change the wettability of the rock surface to water, resulting in changes in the wetting angle. Therefore, the contact angle calculation method disclosed in the above method is not applicable to the determination of the wettability of the cap rock surface of the aquifer hydrogen reservoir, and there is also a problem of low calculation efficiency. Summary of the invention

[0005] In view of the shortcomings of the prior art, the technical problem to be solved by the present invention is to overcome the problem that the existing wettability measurement method cannot accurately reflect the wettability of the surface of the cap rock of the aquifer hydrogen storage reservoir and has low calculation efficiency, and proposes a method for measuring the surface wettability of the cap rock of the aquifer hydrogen storage reservoir, which has the function of calculating the wetting angle of water molecules on the mineral surface under a hydrogen environment under different mineral materials and has high calculation accuracy.

[0006] In order to solve the technical problem, the technical solution adopted by the present invention is:

[0007] The present invention provides a method for determining the surface wettability of a cap rock of an aquifer hydrogen storage reservoir, comprising:

[0008] The modeling steps include establishing a semi-cylindrical water molecule aggregation model on the surface of the cap rock mineral model, constructing a hydrogen-water-cap wettability molecular model, and the plane of the semi-cylindrical water molecule aggregation model in the hydrogen-water-cap wettability molecular model is in linear contact with the cap rock mineral model;

[0009] In the analysis of caprock rocks of aquifer hydrogen reservoirs, the wetting angle not only involves the interaction between water and the rock surface, but also needs to consider the competitive adsorption and interface effects between hydrogen and water. Since hydrogen molecules are smaller in gaseous state, their interaction with water is more complex than that of common hydrocarbon molecules, with stronger polarity and unique molecular dynamics behavior. Compared with hydrocarbon molecules, the interaction between hydrogen and water is significantly different at the physical and chemical level, resulting in more complex wetting behavior. The combined action of water and hydrogen may affect the adsorption behavior of water on the rock wall, especially in the caprock of hydrogen reservoirs. The presence of hydrogen may change the wettability of the rock surface to water, resulting in changes in the wetting angle. The above modeling steps establish a semi-cylindrical water molecule cluster structure on the surface of the mineral model, saturate the remaining space with hydrogen molecules, and establish a contact surface model of hydrogen-water-mineral three-phase coexistence, which fully considers the influence of hydrogen. At the same time, the present invention uses semi-cylindrical water molecule aggregates to study wettability. Compared with the use of hemispherical water molecule aggregates to calculate the wetting angle, the contact line between the semi-cylindrical water molecule aggregate and the mineral surface is a straight line, which can reduce the influence of line tension caused by the curvature of the contact line and improve the calculation accuracy of the wetting angle.

[0010] The simulation steps include dividing the simulation area and boundary conditions; selecting force field functions to describe the interaction forces within and between molecules; optimizing the structure of the contact surface model of hydrogen-water-cap layer coexistence; determining relevant simulation parameters and ensembles; and performing molecular dynamics simulation on the hydrogen-water-cap layer wettability molecular model to obtain wettability simulation results.

[0011] The calculation steps include:

[0012] The water molecule aggregate is evenly divided into multiple water molecule layers along the y-axis direction. Then, based on the Delaunay triangulation algorithm, with the help of finer grid division, the boundary contour points of each water molecule layer are calculated to obtain the boundary contour points of the water molecule aggregate, which enhances the processing ability of complex surfaces and droplet morphology. The cross product characteristic between two vectors is introduced to screen and remove the boundary points of the cover surface and the points with weak direct correlation with the hydrogen-water interface, accurately locate the interface contour points where the droplet and hydrogen are in contact, and accurately describe the geometric shape of the droplet surface, effectively avoiding the simple geometric approximation error caused by the traditional method that relies on the calculation of the highest point of the droplet. In the case of complex droplet surface morphology, it can provide higher accuracy;

[0013] In the calculation of wetting angle, the contact angle between the droplet and the solid surface usually presents a nonlinear relationship. The Levenberg-Marquardt algorithm is used to select multiple function models for optimization fitting to minimize the error of the calculation model. R 2 and MSE Evaluate the optimal fitting function f(x) to parameterize the droplet surface shape, thereby more accurately characterizing the droplet surface morphology and calculating the droplet contact angle, thereby improving calculation accuracy.

[0014] For the existing calculation method, the above calculation steps use finer grid division to calculate the boundary contour points of each water molecule layer, obtain the boundary contour points of the water molecule aggregate, introduce the cross product characteristics between the two vectors, screen and remove the boundary points on the surface of the cover layer and the points with weak direct correlation with the hydrogen-water interface, accurately locate the interface contour points where the droplet and hydrogen are in contact, and accurately describe the geometric shape of the droplet surface; use the Levenberg-Marquardt algorithm, select a variety of function models for optimization fitting, minimize the error of the calculation model, and according to R 2 and MSE Evaluate the optimal fitting function f(x) to parameterize the droplet surface shape.

[0015] According to the boundary points of the outer convex hull, nonlinear curve fitting is performed based on the Levenberg-Marquardt iterative algorithm, and the contour is further optimized to accurately characterize the complex boundary morphology of water molecule aggregates; in this process, it is necessary to select multiple models to fit the curve, and by comparing the fitting effects of different models, select the function form that can most accurately characterize the complex boundary morphology of water molecule aggregates. Then, the Levenberg-Marquardt iterative algorithm is used to perform nonlinear least squares optimization on the selected fitting function. The algorithm adjusts the fitting parameters repeatedly and iteratively to minimize the error between the actual boundary data and the fitting curve, and gradually converges to the optimal solution, thereby ensuring the high accuracy and reliability of the fitting results.

[0016] The above method is based on the microscopic scale. By conducting molecular dynamics simulation research on the wetting behavior of the three-phase coexistence of hydrogen, water and minerals, the wettability characteristics of hydrogen, water and minerals under different cap minerals, temperatures and pressure conditions can be determined, which can provide support for characterizing the wetting heterogeneity of the cap layer of the aquifer hydrogen storage reservoir and lay the foundation for evaluating the cap layer sealing of the aquifer hydrogen storage reservoir. In addition, the present invention uses a model in which a semi-cylindrical water molecule aggregate is wrapped by hydrogen molecules, which is closer to the actual state when the two-phase fluid initially contacts the solid surface. The surface wettability of the minerals in the aquifer can be determined, and the wettability of different components of the cap layer minerals can be determined, which is conducive to the analysis of the sealing of the cap layer of the hydrogen storage reservoir.

[0017] Furthermore, a semi-cylindrical water molecule aggregate model was established by the following method: constructing an initial model of a semi-cylindrical water molecule aggregate, and filling it with different numbers of hydrogen molecules to obtain a semi-cylindrical water molecule aggregate model under different pressure conditions in a hydrogen environment.

[0018] Preferably, the modeling step specifically includes:

[0019] S1: Use Packmol software to construct water molecule model, hydrogen molecule model and cap rock mineral model. The cap rock mineral model includes graphene model, montmorillonite model and quartz model. S2: According to the water molecule, hydrogen molecule and cap rock mineral model, a semi-cylindrical water molecule aggregate structure is established on the surface of the mineral model to construct the hydrogen-water-cap rock wettability molecular model. Among them, the fixed base is composed of the cap rock mineral model. In order to ensure the quantification of water molecule aggregates under different reservoir conditions, the number of water molecules is set to 1000. According to the water molecule model, a semi-cylindrical water molecule aggregate initial model is constructed, and different numbers of hydrogen molecules are filled to obtain the semi-cylindrical water molecule aggregate model under different pressures (5 MPa, 15 MPa, 20 MPa, 30 MPa, 40 MPa) in the hydrogen environment; graphene is selected as the organic mineral to construct a single-layer or multi-layer graphene base; it is spliced ​​with the water molecule aggregate model to obtain the hydrogen-water-graphene wetting initial model; or montmorillonite crystals are selected as clay minerals to construct a montmorillonite base, which is spliced ​​with the water molecule aggregate model to obtain the hydrogen-water-montmorillonite wetting initial model; or quartz crystals are selected as siliceous minerals to construct a quartz base, which is spliced ​​with the water molecule aggregate model to obtain the hydrogen-water-quartz wetting initial model.

[0020] Preferably, the simulation step specifically includes:

[0021] S3: Divide the simulation area and boundary conditions: the base molecule is set as a fixed atomic layer to prevent the loss of atoms due to deformation caused by heat. The model's X and Y directions are set as periodic boundaries, and a soft repulsive wall is set in the Z direction.

[0022] S4: Select force field functions to describe the interaction forces within and between molecules, including: using CVFF force field to describe the interaction between water molecules, hydrogen molecules and graphene surface, using LJ parameters to calculate the interaction between hydrogen molecules, using SHAKE algorithm in LAMMPS to constrain bond angles and bond lengths, and using SPC / E model to model water molecules, using CLAYFF force field to calculate montmorillonite crystals and quartz crystals, some force field parameters are shown in Table 1:

[0023]

[0024]

[0025] The calculation of LJ potential energy , Using the Lorentz-Berthelot mixing rule, the formula is:

[0026]

[0027] Table 1 Some force field parameters

[0028]

[0029] S5: Structural optimization of the contact surface model of hydrogen-water-cap layer coexistence, including: fixing the substrate molecules to maintain stability, minimizing the energy of water molecule aggregates and hydrogen molecules, and obtaining a stable hydrogen-water-cap layer wetting model.

[0030] S6: Determine relevant simulation parameters and ensemble, including: simulation time step of 1.0 fs, cutoff radius of 10.0 Å, initial temperatures of 292 K, 313 K, 323 K, and 343 K, and pressures of 5 MPa, 15 MPa, 20 MPa, 30 MPa, and 40 MPa, respectively; set the simulation system to use the NVT ensemble, use the Nose-Hover heat bath method to control the temperature, and use the number of hydrogen molecules to control the pressure.

[0031] S7: Molecular dynamics simulation of the wetting model was performed to obtain wettability simulation results, including: molecular simulation using LAMMPS software, firstly performing a 1 ns relaxation equilibrium on the optimized structural model, then performing a 5 ns simulation calculation, and exporting the molecular space trajectory coordinate file at each moment, calculating the interaction energy between the water molecule aggregate and the cap layer wall, and the interaction energy between hydrogen and the cap layer wall during the simulation time.

[0032] Preferably, the calculation step specifically includes:

[0033] S8: Simulate and visualize the two-dimensional density distribution cloud map of water molecule aggregates in the hydrogen-water-cap layer wetting model, and intuitively show the wetting state of water molecule aggregates, including:

[0034] S81: Import the trajectory file obtained by simulation into OVITO to visualize the simulated morphology of each time period; due to the irregular motion of molecules, the system is usually in a dynamic equilibrium state, and the morphology of water molecule aggregates fluctuates over time, resulting in the wetting angle at a certain moment cannot be used as the final calculation result. Therefore, the coordinates of the stabilized droplets are averaged to obtain an average coordinate representing the stable morphology of the droplets, and the hydrogen-water-cap layer wettability is calculated.

[0035] S82: Based on the simulation results, draw the interaction energy curve between the water molecule aggregates and the cap layer wall, find the dynamically stable time period, extract the spatial coordinates of the water molecule aggregates in this time period, and perform statistical averaging to obtain a spatial coordinate file that characterizes the stable form of the water molecule aggregates; in this step, in order to ensure that the water molecule aggregates are in a stable state, it is necessary to count multiple stable state time periods and select multiple stable state segments on different time scales. By statistically analyzing the data of time periods of different lengths, the short-term fluctuations of the system can be effectively eliminated to ensure that the results accurately reflect the true stable form of the water molecule aggregates.

[0036] S83: Based on the spatial coordinate file characterizing the water molecule aggregates, calculate the density distribution of the water molecule aggregates in the hydrogen-water-cap layer wetting model on the xz plane, and draw a two-dimensional density distribution cloud map to intuitively show the wetting results of the water molecule aggregates on the cap layer wall.

[0037] S9: Based on Delaunay triangulation and combined with the cross product characteristics of vectors, the boundary points of water molecule aggregates are calculated;

[0038] S10: The wettability of the rock surface of the cap rock of the aquifer hydrogen storage reservoir was determined by iteratively fitting the boundary contour curve using the Levenberg-Marquardt algorithm;

[0039] The above S9 and S10 specifically include:

[0040] Firstly, the height and contact length of water molecule aggregates are obtained according to the density distribution cloud map of water molecule aggregates, and the approximate range of the wetting angle of water molecule aggregates on the surface of the cap rock is estimated using the following formula.

[0041]

[0042] Furthermore, according to the spatial coordinate file characterizing the water molecule aggregate, it is evenly divided into multiple water molecule layers along the y-axis direction, and the boundary coordinates of each layer are calculated. The specific process is as follows:

[0043] According to the spatial coordinate file characterizing the stable morphology of water molecule aggregates, the aggregates are equally divided into multiple water molecule layers along the y-axis. For a single water molecule layer, the y-axis coordinate information is ignored, and efficient meshing is performed based on the Delaunay triangulation method according to the coordinate information of the x-axis and z-axis to obtain the boundary contour points of the layer, and the boundary point coordinate file of all water molecule layers is obtained;

[0044] According to the boundary point coordinate file, the cross product characteristics between two vectors are introduced to screen and remove the boundary points on the surface of the cap layer and the points with weak direct correlation with the hydrogen-water interface, obtain the boundary points of the interface between the water molecule aggregate and hydrogen, and obtain the accurate geometric contour of the hydrogen-water interface;

[0045] According to the boundary points of the hydrogen-water interface, with the help of the Levenberg-Marquardt iterative algorithm, multiple fitting functions are selected to implement nonlinear curve fitting, and the minimized fitting parameters are obtained. Then the relevant contours are further optimized, and several parameterized interface curves are obtained by parameterizing the discrete points on the interface.

[0046] According to several interface curves, calculate R 2 , MSE In order to determine the curve with the best fitting effect, the accurate parametric characterization of the complex boundary morphology of water molecule aggregates can be achieved.

[0047] Based on the parameterized curves that characterize the complex boundary morphology of water molecule aggregates f(x) , calculate the slope of the curve at the point (x0, y0) where the water molecule aggregate intersects the fitting curve at the water-cap layer interface k=fˊ(x 0 ,y 0 ) , accurately solve the wetting angle of water molecule aggregates on the surface of caprock;

[0048] Optionally, the contact surface of the droplet is determined based on the radial distribution function between the cap layer molecules and the water molecules.

[0049] The above calculation steps simplify the calculation process by statistically simulating the stabilized droplet coordinates and taking their average value as the basis for calculating the wetting angle, thereby avoiding complex calculations at each moment and significantly improving the calculation efficiency. At the same time, using stabilized data for wetting angle calculation reduces the impact of droplet morphology fluctuations, can more accurately reflect the final wettability, and makes the calculation results more stable and reliable.

[0050] Compared with the prior art, the present invention has the following beneficial effects:

[0051] (1) The present invention provides a molecular simulation method for determining the wettability of mineral surfaces in a hydrogen environment, which can calculate the wetting angle of water molecules under different mineral materials;

[0052] (2) A semi-cylindrical geometric structure is used, and its straight contact line between water and substrate can eliminate the influence of the finite contact line of spherical water droplets and improve the calculation accuracy;

[0053] (3) The Delaunay triangulation algorithm is combined with the LM (Levenberg-Marquardt) algorithm to achieve accurate fitting of the boundaries of water molecule aggregates. The Delaunay triangulation algorithm provides an efficient grid division scheme and introduces the cross product of vectors to ensure the geometric accuracy of the boundary contour. The LM algorithm further accurately fits the complex boundary morphology of water molecule aggregates through nonlinear optimization and more accurately characterizes the interface shape of water molecule aggregates, thereby improving the accuracy of the calculation of the wetting angle of water molecule aggregates on the cap layer surface.

[0054] (4) The statistical averaging method is used to calculate the wetting angle based on the average coordinate file of the stable morphology of water molecule aggregates, thereby avoiding the complex process of frequently calculating the contact angle in each time step. This method not only improves the calculation efficiency and reduces the calculation time, but also effectively eliminates the influence of the morphology fluctuations of water molecule aggregates, thereby improving the calculation accuracy. Through statistical and averaging processing of the stable state, the wettability of the hydrogen-water-cap layer interface can be more accurately characterized. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 A flow chart for determining the wettability of a mineral surface in a hydrogen environment provided by the present invention

[0056] Figure 2 Schematic diagram of the water molecule structure;

[0057] Figure 3 This is the initial model of hydrogen-water-graphene;

[0058] Figure 4 Optimizing models for hydrogen-water-graphene;

[0059] Figure 5 A single-layer front view of the graphene substrate model used in the present invention;

[0060] Figure 6 A multi-layer front view of the graphene substrate model used in the present invention;

[0061] Figure 7 A top view of the graphene substrate model used in the present invention;

[0062] Figure 8A front view of a montmorillonite substrate model used in the present invention;

[0063] Fig. 9 A top view of the montmorillonite substrate model used in the present invention;

[0064] Fig.10 A front view of a quartz substrate model used in the present invention;

[0065] Fig.11 A top view of a quartz substrate model used in the present invention;

[0066] Fig.12 It is a schematic diagram of the final moment of the simulation snapshot of the hydrogen-water-graphene wetting model of the present invention;

[0067] Fig.13 The basis for selecting the statistical average time step for the present invention is the interaction energy during the simulation process;

[0068] Fig.14 It is a density cloud diagram of water molecule aggregates of the present invention;

[0069] Fig.15 A schematic diagram and parameter indication for estimating the wetting angle of the present invention;

[0070] Fig.16 The process of calculating the boundary contour points of water molecule aggregates of the present invention, (a) schematic diagram of water molecule aggregates, wherein (a1) is the layering method of water molecule aggregates, (a2) is the Delaunay triangulation and boundary points of a single layer of water molecules, (b) is the boundary contour points calculated by Delaunay, and (c) is the boundary points after vector cross product optimization;

[0071] Fig.17 The best interface curve fitted by LM;

[0072] Fig.18 Comparison between the two-dimensional density cloud map and the fitted interface curve, as well as a schematic diagram of the wetting angle calculation. DETAILED DESCRIPTION

[0073] The technical scheme in the specific embodiment of the present invention will be described in detail and completely in conjunction with the accompanying drawings. Obviously, the described embodiments are only some specific implementation methods of the overall technical scheme of the present invention, rather than all implementation methods. Based on the overall concept of the present invention, all other embodiments obtained by ordinary technicians in this field fall within the scope of protection of the present invention.

[0074] The present invention provides a molecular simulation method for determining the wettability of mineral surfaces in a hydrogen environment to calculate the wetting angle of water molecules under different mineral materials. Figure 1 As shown, the method includes:

[0075] S1: construct water molecule, hydrogen molecule and cap rock mineral models respectively: the cap rock mineral models include graphene model, montmorillonite model and quartz model;

[0076] S2: Based on the water molecules, hydrogen molecules and cap rock mineral models, a semi-cylindrical water molecule cluster structure is established on the surface of the mineral model to construct a hydrogen-water-cap rock wettability molecular model;

[0077] S3: Divide the simulation area and boundary conditions: the base molecule is set as a fixed atomic layer to prevent the loss of atoms due to deformation caused by heat. The model's X and Y directions are set as periodic boundaries, and a soft repulsive wall is set in the Z direction.

[0078] S4: Select force field functions to describe the interaction forces within and between molecules;

[0079] S5: Structural optimization of the contact surface model of hydrogen-water-cap layer three-phase coexistence;

[0080] S6: Determine relevant simulation parameters and ensemble;

[0081] S7: Perform molecular dynamics simulation on the wetting model to obtain wettability simulation results;

[0082] S8: Simulation visualization, calculating the two-dimensional density distribution cloud map of water molecule aggregates in the hydrogen-water-cap layer wetting model, intuitively showing the wetting state of water molecule aggregates;

[0083] S9: Based on Delaunay triangulation and combined with the cross product characteristics of vectors, the boundary points of water molecule aggregates are calculated;

[0084] S10: The Levenberg-Marquardt algorithm is used to iteratively fit the boundary contour curve to determine the wettability of the rock surface of the cap rock of the aquifer hydrogen storage reservoir.

[0085] The specific implementation is as follows:

[0086] The cap layer of the aquifer hydrogen storage reservoir is composed of organic matter, clay minerals, siliceous minerals, etc., with graphene as a representative of organic matter, montmorillonite as a representative of clay minerals, and quartz as a representative of siliceous minerals. In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0087] The graphene model was established in Materials Studio. The single-layer graphene unit cell was obtained and transformed from an oblique state to an orthogonal state and expanded in the x and y directions to construct a single-layer graphene substrate model with side lengths of 10.45 nm and 56 nm, respectively. Figure 5 , 7As shown;

[0088] The single-layer graphene unit cell is expanded along the z-axis by two or three layers with an interlayer spacing of 0.335 nm to obtain a multilayer graphene model, such as Figure 6 , 7 As shown;

[0089] Montmorillonite crystals are used as the clay mineral substrate. The position parameters of the initial unit cell are obtained from AMSCD (American Mineralogists Crystal Structure Database) and imported into Materials Studio for cell expansion to obtain the montmorillonite crystals used in the simulation, such as Figure 8 , 9 As shown;

[0090] The quartz crystal is used as the siliceous mineral substrate. The position parameters of the initial unit cell are obtained from AMSCD (American Mineralogists Crystal Structure Database), imported into Materials Studio for cell expansion, and the surface of the quartz crystal is hydroxylated to obtain the quartz crystal model used in the simulation, such as Fig.10 , 11 As shown;

[0091] Step S1 specifically includes:

[0092] Selecting graphene as organic matter to construct a graphene substrate;

[0093] Montmorillonite crystals were selected as clay minerals to construct a montmorillonite matrix;

[0094] Selecting quartz crystal as the siliceous mineral to construct the quartz substrate;

[0095] Packmol was used to construct water molecule models and hydrogen molecule models.

[0096] Step S2 specifically includes:

[0097] Based on the density of water being 1.0 g / cm 3 , 1000 such Figure 2 The water molecules shown are placed in a semi-cylindrical space with a radius of 1.8 nm and a length of 5.6 nm by Packmol software to establish an initial model of water molecule aggregates;

[0098] According to the selected simulation pressure of 5 MPa and simulation temperature of 292 K, referring to the NIST (National Institute of Standards and Technology) database, see Table 2, 565 hydrogen molecules are filled above the semi-cylindrical water molecule aggregate to construct a hydrogen-water model, as shown in Figure 3 As shown;

[0099] Table 2 Number of hydrogen molecules under some working conditions (pieces)

[0100]

[0101] The hydrogen-water model and the substrate (graphene substrate, montmorillonite substrate, quartz substrate) were combined in Packmol software to obtain the hydrogen-water-cap layer initial wetting model with side lengths of 10.45 nm, 56 nm, and 10.6 nm, respectively. Figure 3 As shown;

[0102] Steps S3-S7, performing molecular simulation in LAMMPS to obtain wettability simulation results, specifically include:

[0103] The "hydrogen-water-graphene" system (hydrogen-water-graphene wetting model) uses the CVFF force field to describe the intermolecular interaction force. The hydrogen molecules use the LJ parameters, the water molecules use the SPC / E model, and the FIX SHAKE is used to constrain the bond length and bond angle of the water molecules. The simulation calculation process is as follows: first initialize the parameters, read the model file, assign force field parameters and group them; fix the base atoms, and minimize the energy of the initial model, remove atomic overlap, etc., to obtain the energy-optimized model, such as Figure 4 As shown; the simulation step size is set to 1.0 fs, the cutoff radius is 10.0Å, the initial temperature is given to 292 K, the NVT ensemble is used for hydrogen and water molecules, the Nose-Hoover heat bath method is used for temperature control, the pressure is controlled by the number of hydrogen molecules input into the model, periodic boundaries in the x and y directions, fixed boundaries in the z direction, a soft repulsive wall is set, and the system is subjected to 1 ns equilibrium relaxation; finally, the same settings are used to perform a 5 ns molecular simulation of the system, and the spatial trajectory coordinate files within the simulation time are collected, the interaction energy between the cap layer and water and the cap layer and hydrogen during the simulation process is calculated, and the two-dimensional data of the volume density of water molecule aggregates is calculated;

[0104] "Hydrogen-water-montmorillonite" system (hydrogen-water-montmorillonite wetting model), the calculation of montmorillonite crystals uses the CLAYFF force field and is simulated using LAMMPS. The simulation parameter settings are the same as those of the "hydrogen-water-graphene" system;

[0105] In the "hydrogen-water-quartz" system (hydrogen-water-quartz wetting model), the calculation of the quartz crystal adopts the CLAYFF force field and is simulated in LAMMPS. The simulation parameters are set the same as those of the "hydrogen-water-quartz" system.

[0106] During the simulation, water molecule aggregates on the cap layer substrate will drift left and right, deviate from the center position, and even cross the periodic boundary to cause disconnection, affecting the visual analysis of the simulation and the calculation of the contact angle. In LAMMPS, the fix recenter command is used to fix the coordinates of the water molecule aggregates so that the water molecule aggregates are automatically fixed to the center position of the substrate surface during the simulation.

[0107] Steps S8-S10, visualize the simulation results, calculate the two-dimensional cloud map of the volume density distribution of water molecule aggregates in the hydrogen-water-cap layer wetting model, and calculate the contact angle of water molecule aggregates on the surface of the cap layer mineral, specifically including:

[0108] The trajectory file obtained by simulation is imported into OVITO to visualize the simulated morphology of each time period. Due to the irregular motion of molecules, the system is usually in a dynamic equilibrium state. The morphology of water molecule aggregates fluctuates over time, resulting in the wetting angle at a certain moment not being able to be used as the final calculation result. Therefore, the coordinates of the stabilized droplets are averaged to obtain an average coordinate representing the stable morphology of the droplets, and the wettability of hydrogen-water-cap layer is calculated.

[0109] According to the simulation results, the interaction energy curve between the water molecule aggregate and the cap layer wall is drawn, such as Fig.13 As shown, the dynamically stable time period is found, the spatial coordinates of the water molecule aggregates in this time period are extracted, and statistical averaging is performed to obtain the spatial coordinate file that characterizes the stable morphology of the water molecule aggregates. The stable morphology of the water molecule aggregates is shown in Fig.12 As shown;

[0110] In order to accurately characterize the stable morphology of water molecule aggregates, it is necessary to select time steps of different lengths for statistical averaging in a stable time period to ensure statistical accuracy;

[0111] According to the spatial coordinate file characterizing the water molecule aggregates, the density distribution of the water molecule aggregates in the hydrogen-water-cap layer wetting model on the xz plane is calculated, and a two-dimensional density distribution cloud map is drawn to intuitively show the wetting results of the water molecule aggregates on the cap layer wall, such as Fig.14 As shown;

[0112] First, according to the density distribution cloud map of water molecule aggregates, the height and contact length of water molecule aggregates are obtained, such as Fig.15 As shown, the following formula is used to preliminarily estimate the approximate range of the wetting angle of water molecule aggregates on the surface of the cap rock.

[0113]

[0114] Furthermore, according to the spatial coordinate file characterizing the water molecule aggregate, it is evenly divided into multiple water molecule layers along the y-axis direction, and the boundary coordinates of each layer are calculated. The specific process is as follows: Fig.16 As shown:

[0115] According to the spatial coordinate file that characterizes the stable morphology of water molecule aggregates, it is divided into multiple water molecule layers along the y-axis. For a single water molecule layer, the y-axis coordinate information is ignored, and the boundary contour points of the water molecules in this layer are obtained based on the Delaunay triangulation method according to the coordinate information of the x-axis and z-axis. The boundary point coordinate file of all water molecule layers is obtained, such as Fig.16 As shown in (b);

[0116] According to the boundary point coordinate file, the cross product characteristics between two vectors are introduced to screen and remove the boundary points on the surface of the cap layer and the points with weak direct correlation with the hydrogen-water interface, obtain the boundary points of the interface between the water molecule aggregate and hydrogen, and obtain the precise geometric contour of the hydrogen-water interface, such as Fig.16 As shown in (c);

[0117] According to the boundary points of the hydrogen-water interface, with the help of the Levenberg-Marquardt iterative algorithm, multiple fitting functions are selected to implement nonlinear curve fitting, and the minimized fitting parameters are obtained, and then the relevant profiles are further optimized. By parameterizing the discrete points on the interface, several parameterized interface curves are obtained.

[0118] According to several interface curves, calculate R 2 , MSE To determine the curve with the best fitting effect, such as Fig.17 As shown, accurate parameterized characterization of the complex boundary morphology of water molecule aggregates is achieved.

[0119] Based on the parameterized curves that characterize the complex boundary morphology of water molecule aggregates f(x) , calculate the slope of the curve at the point (x0, y0) where the water molecule aggregate intersects the fitting curve at the water-cap layer interface k=fˊ(x 0 ,y 0 ) , accurately solve the wetting angle of water molecules aggregates on the surface of caprock, such as Fig.18 As shown;

[0120] The contact surface of the droplet was determined based on the radial distribution function between the cap layer molecules and the water molecules, and the results are shown in Table 3.

[0121] Table 3 Wetting angles of some conditions of the “hydrogen-water-graphene” system

[0122]

Claims

1. A method for determining the wettability of the rock surface of the cap rock of an aquifer hydrogen storage reservoir, characterized in that: include: The modeling step includes establishing a semi-cylindrical water molecule aggregation model on the surface of the cap rock mineral model and constructing a hydrogen-water-cap wettability molecular model, wherein the plane of the semi-cylindrical water molecule aggregation model in the hydrogen-water-cap wettability molecular model is in linear contact with the cap rock mineral model; The simulation steps include dividing the simulation area and boundary conditions; selecting a force field function to describe the interaction forces within and between molecules; optimizing the structure of the contact surface model of hydrogen-water-cap layer coexistence; determining relevant simulation parameters and ensembles; and performing molecular dynamics simulation on the hydrogen-water-cap layer wettability molecular model to obtain wettability simulation results; The calculation steps include: The water molecule aggregate is evenly divided into multiple water molecule layers along the y-axis direction. Then, based on the Delaunay triangulation algorithm, with the help of finer grid division, the boundary contour points of each water molecule layer are calculated to obtain the boundary contour points of the water molecule aggregate. The cross product characteristic between two vectors is introduced to screen and remove the boundary points on the cover layer surface and the points with weak direct correlation with the hydrogen-water interface, accurately locate the interface contour points where the droplet and hydrogen are in contact, and accurately describe the geometric shape of the droplet surface. The Levenberg-Marquardt algorithm is used to select multiple function models for optimization fitting, minimize the error of the calculation model, and R 2 and MSE Evaluate the optimal fitting function f(x) to parameterize the droplet surface shape; The modeling steps specifically include: S1: using Packmol software to construct a water molecule model, a hydrogen molecule model and a cap rock mineral model, wherein the cap rock mineral model includes a graphene model, a montmorillonite model and a quartz model; S2: According to the water molecules, hydrogen molecules and cap rock mineral model, the semi-cylindrical water molecule aggregation model is established on the surface of the mineral model to construct the hydrogen-water-cap rock wettability molecular model; The simulation steps specifically include: S3: Divide the simulation area and boundary conditions: the base molecule is set as a fixed atomic layer, the model X and Y directions are set as periodic boundaries, and the Z direction is set as a soft repulsive wall; S4: Select force field functions to describe the interaction forces within and between molecules; S5: Structural optimization of the contact surface model of hydrogen-water-cap layer three-phase coexistence; S6: Determine relevant simulation parameters and ensemble; S7: Perform molecular dynamics simulation on the wetting model to obtain wettability simulation results; The calculation steps specifically include: S8: Simulation visualization, calculation of the two-dimensional density distribution cloud map of water molecule aggregates in the hydrogen-water-cap layer wetting model; S9: Based on Delaunay triangulation and combined with the cross product characteristics of vectors, the boundary points of water molecule aggregates are calculated; S10: The Levenberg-Marquardt algorithm is used to iteratively fit the boundary contour curve to determine the wettability of the rock surface of the cap rock of the aquifer hydrogen storage reservoir.

2. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 1, characterized in that: The semi-cylindrical water molecule aggregate model is established by the following method: constructing a semi-cylindrical water molecule aggregate initial model, and filling it with different numbers of hydrogen molecules to obtain a semi-cylindrical water molecule aggregate model under different pressure conditions in a hydrogen environment.

3. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 1, characterized in that: The calculation of the boundary contour points of each water molecule layer specifically includes: calculating the boundary points of each water molecule layer, screening and optimizing all boundary points to obtain contour points that accurately characterize the boundary of the water molecule aggregate; The selecting of multiple function models for optimization fitting specifically includes: selecting multiple models to fit the curve, and selecting a function form that can accurately characterize the complex boundary morphology of water molecule aggregates by comparing the fitting effects of different models.

4. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 1, characterized in that: The S8 specifically includes: S81: Import the trajectory file obtained by simulation into OVITO to visualize the simulation morphology of each time period; S82: Based on the simulation results, draw the interaction energy curve between the water molecule aggregates and the cap layer wall, find the dynamically stable time period, extract the spatial coordinates of the water molecule aggregates in this time period, perform statistical averaging, and obtain the spatial coordinate file that characterizes the stable form of the water molecule aggregates; S83: Based on the spatial coordinate file characterizing the water molecule aggregates, the density distribution of the water molecule aggregates in the hydrogen-water-cap layer wetting model on the xz plane is calculated, and a two-dimensional density distribution cloud map is drawn to intuitively display the wetting results of the water molecule aggregates on the cap layer wall.

5. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 4, characterized in that: The S82 further includes: By counting multiple stable state time periods and selecting multiple stable state segments on different time scales, the data of time periods of different lengths are counted to obtain a spatial coordinate file that accurately reflects the true stable morphology of water molecule aggregates.

6. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 1, characterized in that: The S4 includes: using CVFF force field to describe the interaction between water molecules, hydrogen molecules and graphene surface, using LJ parameters to calculate the force between hydrogen molecules, using SHAKE algorithm in LAMMPS to constrain bond angles and bond lengths, and using SPC / E model to model water molecules, and using CLAYFF force field to calculate montmorillonite crystals and quartz crystals; The S5 includes: fixing the substrate molecules, minimizing the energy of water molecule aggregates and hydrogen molecules, and obtaining a stable hydrogen-water-cap layer wetting model.

7. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 1, characterized in that: The S7 includes: using LAMMPS software to perform molecular simulation, firstly performing a relaxation equilibrium of the optimized structural model for 1 ns, then performing a simulation calculation for 5 ns, and exporting the molecular space trajectory coordinate file at each moment, calculating the interaction energy between the water molecule aggregate and the cap layer wall, and the interaction energy between hydrogen and the cap layer wall during the simulation time.

Images