Aquifer hydrogen storage cover rock surface wettability determination method

Abstract

The invention discloses a method for measuring the surface wettability of overlying rock of a hydrogen storage reservoir of an aquifer, and belongs to the technical field of wettability measurement. Comprising the following steps: establishing a semi-cylindrical water molecule aggregate model on the surface of a cover rock mineral model; dividing simulation areas and boundary conditions; selecting a force field function, and describing interaction forces inside molecules and among the molecules; performing structure optimization on the hydrogen-water-cover layer three-phase coexisting contact surface model; related simulation parameters and ensembles are determined; and performing molecular dynamics simulation on the hydrogen-water-cover layer wettability molecular model to obtain a wettability simulation result. The method is applied to the aspect of measuring the surface wettability of the overlying rock of the hydrogen storage reservoir of the aquifer, and solves the problems that an existing wettability measuring method cannot accurately reflect the surface wettability of the overlying rock of the hydrogen storage reservoir of the aquifer and is low in calculation efficiency. The method has the characteristics of being capable of calculating the wetting angles of the water molecules under the mineral surface hydrogen environment of different mineral materials and high in calculation precision.

Description

Technical Field

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

[0002] The wettability of the caprock of an aquifer hydrogen storage reservoir refers to the tendency of formation water in hydrogen to spread or adhere to the mineral surface of the caprock. It is one of the most important properties of the hydrogen storage caprock surface. Wettability has a significant impact on the flow of hydrogen and water phases and the sealing performance of the caprock, playing a crucial role in the site selection and operation of aquifer hydrogen storage reservoirs.

[0003] Traditional wettability characterization methods primarily rely on macroscopic qualitative analysis, while quantitative studies of caprock wettability at the meso-micro scale, particularly in hydrogen, are limited. Chinese patent CN112382345A discloses a method and system for calculating contact angles in shale oil reservoirs. This system constructs a shale oil reservoir rock model, a shale oil molecule model, and a water molecule model, building a microscopic oil-water-rock wettability molecular model. Molecular dynamics simulations are then performed on the oil-water-rock wettability molecular model. Water molecules in the model are extracted from the simulation results and clustered to form water clusters. These clusters are then digitized to determine the contact angles of the digitized water clusters on the shale oil reservoir rock surface at each simulation moment. This method overcomes the shortcomings of prior art methods, such as the inability of two-dimensional images to accurately reflect the morphological characteristics of water clusters and the subjective interference of manual angle measurement. By obtaining molecular coordinates and extracting the parameters required for contact angle calculation, the accuracy of contact angle calculation is significantly improved.

[0004] However, the analysis object of the above method is shale oil reservoirs. In the calculation of the wetting angle of shale oil reservoirs, 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 aquifer hydrogen storage reservoir, the wetting angle not only involves the interaction between water and the rock surface, but also needs to consider the competitive adsorption and interfacial effects between hydrogen and water. Because 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 storage reservoir. The presence of hydrogen may change the wettability of the rock surface to water, resulting in a change 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 storage reservoir, and there is also the problem of low calculation efficiency. Summary of the Invention

[0005] In view of the shortcomings of the existing technology, 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. A method for measuring the surface wettability of the cap rock of the aquifer hydrogen storage reservoir is proposed, which can calculate the wetting angle of water molecules on the mineral surface under 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: The present invention provides a method for measuring the surface wettability of a cap rock of an aquifer hydrogen storage reservoir, comprising: The modeling steps include establishing a semi-cylindrical water molecule aggregate model on the surface of the cap rock mineral model, constructing a hydrogen-water-cap layer wettability molecular model, wherein the plane of the semi-cylindrical water molecule aggregate model in the hydrogen-water-cap layer wettability molecular model is in linear contact with the cap rock mineral model; In the analysis of caprock rocks in aquifer hydrogen reservoirs, the wetting angle involves not only the interaction between water and the rock surface, but also the competitive adsorption and interfacial effects between hydrogen and water. Because hydrogen molecules are smaller in their gaseous state, their interactions with water are more complex than those of common hydrocarbon molecules, with stronger polarity and unique molecular dynamics. Compared to hydrocarbon molecules, the interactions between hydrogen and water differ significantly 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. In particular, in the caprock of hydrogen reservoirs, the presence of hydrogen may alter the wettability of the rock surface to water, leading to 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, fully considering 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.

[0007] The simulation steps include dividing the simulation area and boundary conditions; selecting force field functions to describe the interaction forces within and between molecules; structural optimization of the interface 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 calculation steps include: The water molecule aggregate is evenly divided into multiple water molecule layers along the y-axis. Then, based on the Delaunay triangulation algorithm and with the help of a 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. This enhances the ability to handle complex surfaces and droplet morphologies. The cross product property between two vectors is introduced to filter out the boundary points of the cover surface and points with weak direct correlation with the hydrogen-water interface. The interface contour points where the droplet and hydrogen are in contact are accurately located, and the geometric shape of the droplet surface is accurately described. This effectively avoids the simple geometric approximation error caused by the traditional method that relies on the calculation of the droplet's highest point. In the case of complex droplet surface morphology, it can provide higher accuracy. In terms of wetting angle calculation, 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 and 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.

[0008] 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 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 contact, and accurately describe the geometric shape of the droplet surface; use the Levenberg-Marquardt algorithm, select multiple 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.

[0009] Based on the perimeter convex hull boundary points, a nonlinear curve fitting method based on the Levenberg-Marquardt iterative algorithm is used to further optimize the contour to accurately characterize the complex boundary morphology of water molecule aggregates. During this process, multiple models are selected for curve fitting. By comparing the fitting results of different models, the function form that most accurately represents the complex boundary morphology of water molecule aggregates is selected. The selected fitting function is then optimized using a nonlinear least-squares method using the Levenberg-Marquardt iterative algorithm. This algorithm iteratively adjusts the fitting parameters to minimize the error between the actual boundary data and the fitted curve, gradually converging to the optimal solution, thereby ensuring high accuracy and reliability of the fitting results.

[0010] The above method is based on the microscopic scale. By conducting molecular dynamics simulations on the wetting behavior of the coexisting hydrogen-water-mineral three-phase, it can determine the wettability characteristics of hydrogen-water-minerals under different cap minerals, temperatures, and pressure conditions. It can provide support for characterizing the wetting heterogeneity of the cap layer of aquifer hydrogen storage reservoirs and lay the foundation for evaluating the cap layer sealing of aquifer hydrogen storage reservoirs. In addition, the present invention uses a model in which semi-cylindrical water molecule aggregates are wrapped by hydrogen molecules, which is closer to the actual state when the two-phase fluids initially contact the solid surface. It can determine the wettability of the mineral surface of the aquifer and the wettability of different components of the cap layer minerals, which is beneficial to the analysis of the sealing of the hydrogen storage reservoir cap layer.

[0011] 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.

[0012] Preferably, the modeling step specifically includes: S1: Use Packmol software to construct water molecule models, hydrogen molecule models and cap rock mineral models. The cap rock mineral models include graphene model, montmorillonite model and quartz model. S2: Based on the water molecule, hydrogen molecule and cap rock mineral models, a semi-cylindrical water molecule aggregate structure is established on the mineral model surface to construct a hydrogen-water-cap rock wettability molecular model. Among them, the fixed base is composed of a caprock mineral model. To ensure the quantification of water molecule aggregates under different reservoir conditions, the number of water molecules is set to 1000. Based on the water molecule model, a semi-cylindrical water molecule aggregate initial model is constructed, and by filling different numbers of hydrogen molecules, semi-cylindrical water molecule aggregate models under different pressures (5 MPa, 15 MPa, 20 MPa, 30 MPa, and 40 MPa) in a hydrogen environment are obtained; graphene is selected as the organic mineral to construct a single-layer or multi-layer graphene base; this is spliced ​​with the water molecule aggregate model to obtain a 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 a 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 a hydrogen-water-quartz wetting initial model.

[0013] Preferably, the simulation step specifically includes: S3: Divide the simulation area and set boundary conditions: the base molecules are set to a fixed atomic layer to prevent atomic loss due to deformation caused by heat. The model is set to periodic boundaries in the X and Y directions, and a soft repulsive wall is set in the Z direction. S4: Select force field functions to describe the intra-molecular and intermolecular interactions. Specifically, the CVFF force field is used to describe the interactions between water molecules, hydrogen molecules, and the graphene surface. The LJ parameters are used to calculate the intermolecular forces of hydrogen. The SHAKE algorithm is used in LAMMPS to constrain bond angles and bond lengths, and the SPC / E model is used to model water molecules. The CLAYFF force field is used for calculations of montmorillonite and quartz crystals. Some force field parameters are shown in Table 1:

[0014]

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

[0016] Table 1 Some force field parameters

[0017] S5: Structural optimization of the hydrogen-water-cap layer three-phase coexistence interface model, specifically 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.

[0018] 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 hot bath method to control temperature, and use the number of hydrogen molecules to control pressure.

[0019] S7: Molecular dynamics simulation of the wetting model was performed to obtain wettability simulation results, specifically including: molecular simulation using LAMMPS software, first 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 aggregates and the cap layer wall during the simulation time, and the interaction energy between hydrogen and the cap layer wall.

[0020] Preferably, the calculation step specifically includes: S8: Simulate and visualize the two-dimensional density distribution of water molecule aggregates in the hydrogen-water-cap layer wetting model, and intuitively show the wetting state of water molecule aggregates, including: 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.

[0021] S82: Based on the simulation results, draw the interaction energy curve between the water molecule aggregates and the cover 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 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, ensuring that the results accurately reflect the true stable form of the water molecule aggregates.

[0022] S83: Based on the spatial coordinate file representing the water molecule aggregates, calculate the density distribution of 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 water molecule aggregates on the cap layer wall.

[0023] S9: Based on Delaunay triangulation and combined with the cross product property of vectors, the boundary points of water molecule aggregates are calculated; S10: The wettability of the rock surface of the caprock of the aquifer hydrogen reservoir is determined by iteratively fitting the boundary contour curve using the Levenberg-Marquardt algorithm; The above S9 and S10 specifically include: First, the height and contact length of the water molecule aggregates are obtained based on the density distribution cloud map of the water molecule aggregates, and the approximate range of the wetting angle of the water molecule aggregates on the cover rock surface is estimated using the following formula.

[0024]

[0025] Furthermore, according to the spatial coordinate file representing 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: According to the spatial coordinate file representing 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 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; Based on the boundary point coordinate file, the cross product property between two vectors is introduced to screen out boundary points on the cap layer surface and points with weak direct correlation with the hydrogen-water interface, thereby obtaining the boundary points of the phase interface between water molecule aggregates and hydrogen and obtaining the precise geometric contour of the hydrogen-water interface. Based on 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 to obtain the minimized fitting parameters, and then the relevant contours are further optimized. By parameterizing the discrete points on the interface, several parameterized interface curves are obtained.

[0026] 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.

[0027] According to the parameterized curves characterizing 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; Optionally, the contact surface of the droplet is determined based on a radial distribution function between the cap layer molecules and the water molecules.

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

[0029] Compared with the prior art, the present invention has the following beneficial effects: (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; (2) A semi-cylindrical geometric structure is used, and its straight contact line between water and substrate can eliminate the finite contact line effect of spherical water droplets and improve the calculation accuracy; (3) The Delaunay triangulation algorithm and the LM (Levenberg-Marquardt) algorithm are combined to achieve accurate fitting of the boundaries of water molecule aggregates. The Delaunay triangulation 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, more accurately characterizes the interface shape of water molecule aggregates, and thus improves the accuracy of the calculation of the wetting angle of water molecule aggregates on the cover layer surface. (4) The statistical averaging method is used to calculate the wetting angle by using the average coordinate file based on 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 morphological 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

[0030] Figure 1 A flow chart for determining the wettability of a mineral surface in a hydrogen environment provided by the present invention Figure 2 Schematic diagram of the water molecule structure; Figure 3 This is the initial model of hydrogen-water-graphene; Figure 4 Optimizing models for hydrogen-water-graphene; Figure 5 A single-layer front view of the graphene substrate model used in the present invention; Figure 6 A multi-layer front view of the graphene substrate model used in the present invention; Figure 7 A top view of the graphene substrate model used in the present invention; Figure 8 A front view of a montmorillonite substrate model used in the present invention; Figure 9 A top view of the montmorillonite substrate model used in the present invention; Figure 10 A front view of the quartz substrate model used in the present invention; Figure 11 A top view of the quartz substrate model used in the present invention; Figure 12 Schematic diagram of the final moment of the simulation snapshot of the hydrogen-water-graphene wetting model of the present invention; Figure 13 The basis for selecting the statistical average time step in the present invention is the interaction energy during the simulation process; Figure 14 It is a density cloud diagram of water molecule aggregates of the present invention; Figure 15 A schematic diagram and parameter indication for estimating the wetting angle according to the present invention; Figure 16 Figure 1 shows the process for calculating the boundary contour points of water molecule aggregates in the present invention. (a) Schematic diagram of water molecule aggregates, where (a1) is the layering method for 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. Figure 17 The best interface curve fitted by LM; Figure 18 Comparison of the two-dimensional density cloud map and the fitted interface curve, as well as a schematic diagram of the wetting angle calculation. DETAILED DESCRIPTION

[0031] The following is a detailed and complete description of the technical solutions in specific embodiments of the present invention, with reference to the accompanying drawings. It should be understood that the described embodiments are merely partial implementations of the overall technical solution of the present invention, and are not exhaustive. All other embodiments derived by those skilled in the art based on the overall concept of the present invention are intended to fall within the scope of protection of the present invention.

[0032] 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: 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; S2: Based on the water molecule, hydrogen molecule and cap rock mineral model, a semi-cylindrical water molecule cluster structure is established on the mineral model surface to construct a hydrogen-water-cap rock wettability molecular model; S3: Divide the simulation area and set boundary conditions: the base molecules are set to a fixed atomic layer to prevent atomic loss due to deformation caused by heat. The model is set to periodic boundaries in the X and Y directions, and a soft repulsive wall is set in the Z direction. S4: Select force field functions to describe the intra-molecular and inter-molecular interaction forces; S5: Structural optimization of the interface 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; 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; S9: Based on Delaunay triangulation and combined with the cross product property of vectors, the boundary points of water molecule aggregates are calculated; S10: The wettability of the rock surface of the caprock of the aquifer hydrogen storage reservoir is determined by iteratively fitting the boundary contour curve using the Levenberg-Marquardt algorithm.

[0033] The specific implementation is as follows: 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. To make the above-mentioned objects, features, and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0034] The graphene model was established in Materials Studio. A single-layer graphene unit cell was obtained and transformed from an oblique state to an orthorhombic state. The cell was 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 、 7 As shown; 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, as shown in Figure 6 、 7 As shown; 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; Quartz crystal is used as the silica 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 Figure 10 、 11 As shown; Step S1 specifically includes: Graphene is selected as the organic matter to construct a graphene substrate; Montmorillonite crystals were selected as clay minerals to construct the montmorillonite matrix; Quartz crystals were selected as silica minerals to construct the quartz substrate; Packmol was used to construct water molecule models and hydrogen molecule models.

[0035] Step S2 specifically includes: Based on the density of water being 1.0 g / cm 3 , 1000 such Figure 2 The water molecules shown were placed in a semi-cylindrical space with a radius of 1.8 nm and a length of 5.6 nm using Packmol software to establish an initial model of water molecule aggregates; 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; Table 2 Number of hydrogen molecules under some working conditions (pieces)

[0036] The hydrogen-water model and the substrate (graphene substrate, montmorillonite substrate, quartz substrate) were merged 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; Steps S3-S7, performing molecular simulation in LAMMPS to obtain wettability simulation results, specifically include: The "hydrogen-water-graphene" system (hydrogen-water-graphene wetting model) uses the CVFF force field to describe the intermolecular interaction forces. The hydrogen molecules use the LJ parameters, the water molecules use the SPC / E model, and the bond lengths and angles of the water molecules are constrained by FIX SHAKE. The simulation calculation process is as follows: first, the parameters are initialized, the model file is read, the force field parameters are assigned and grouped; the base atoms are fixed, and the initial model is energy minimized, removing atomic overlaps, etc. to obtain the energy-optimized model, such as Figure 4As shown in the figure, 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 temperature control method is the Nose-Hoover heat bath method, the pressure is controlled by the number of hydrogen molecules in the input model, periodic boundaries in the x and y directions and a fixed boundary in the z direction are used, a soft repulsive wall is set, and the system is subjected to equilibrium relaxation for 1 ns; finally, the same settings are used to perform a molecular simulation of the system for 5 ns, and the spatial trajectory coordinate file within the simulation time is 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 density data of the water molecule aggregate is calculated; The "hydrogen-water-montmorillonite" system (hydrogen-water-montmorillonite wetting model) uses the CLAYFF force field to calculate the montmorillonite crystals and is simulated in LAMMPS. The simulation parameters are the same as those for the "hydrogen-water-graphene" system. For the "hydrogen-water-quartz" system (hydrogen-water-quartz wetting model), the quartz crystal is calculated using the CLAYFF force field and simulated in LAMMPS. The simulation parameters are the same as those for the "hydrogen-water-quartz" system. During the simulation, water molecule aggregates on the cap layer substrate may drift left and right, deviating from their center position or even crossing periodic boundaries, causing disconnection. This affects the visual analysis of the simulation and the calculation of contact angles. In LAMMPS, the fix recenter command is used to fix the coordinates of the water molecule aggregates, allowing them to automatically return to the center position of the substrate surface during the simulation. Steps S8-S10, visualizing the simulation results, calculating the two-dimensional cloud map of the volume density distribution of water molecule aggregates in the hydrogen-water-cap layer wetting model, and calculating the contact angle of the water molecule aggregates on the cap layer mineral surface, specifically include: The trajectory file obtained by simulation was 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, and 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 were averaged to obtain an average coordinate representing the stable morphology of the droplets, and the hydrogen-water-cap layer wettability was calculated.

[0037] According to the simulation results, the interaction energy curve between the water molecule aggregates and the cover layer wall is drawn, as shown in Figure 13 As shown, we 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 morphology of the water molecule aggregates. The stable morphology of the water molecule aggregates is shown in Figure 12 As shown; 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 within a stable time period to ensure statistical accuracy; According to the spatial coordinate file representing 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 Figure 14 As shown; 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 Figure 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.

[0038]

[0039] Furthermore, according to the spatial coordinate file representing 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 Figure 16 As shown: 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. According to the coordinate information of the x-axis and z-axis, efficient meshing is performed based on the Delaunay triangulation method to obtain the boundary contour points of the water molecules in this layer, and the boundary point coordinate files of all water molecule layers are obtained, such as Figure 16 As shown in (b); According to the boundary point coordinate file, the cross product characteristics between the two vectors are introduced to screen and remove the boundary points on the cover surface and the points with weak direct correlation with the hydrogen-water interface, and obtain the boundary points of the interface between the water molecule aggregates and hydrogen, and obtain the accurate geometric contour of the hydrogen-water interface, such as Figure 16 As shown in (c); Based on the boundary points of the hydrogen-water interface, the Levenberg-Marquardt iterative algorithm was used to select multiple fitting functions for nonlinear curve fitting, obtaining the minimized fitting parameters and further optimizing the relevant profiles. By parameterizing the discrete points on the interface, several parameterized interface curves were obtained.

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

[0041] According to the parameterized curves characterizing 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 the cap rock, such as Figure 18 As shown; 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.

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

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.

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 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: Based on the water molecules, hydrogen molecules and cap rock mineral model, the 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.

5. 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 simulation steps specifically include: S3: Divide the simulation area and boundary conditions: the substrate molecules are set as fixed atomic layers, 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.

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 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.

7. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 6, 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.

8. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 7, 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.

9. The method for determining the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 5, 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.

10. The method for measuring the surface wettability of the cap rock of an aquifer hydrogen storage reservoir according to claim 5, 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