Hydrate-containing sediment depressurization exploitation numerical simulation method, medium, equipment and product

Abstract

The invention provides a hydrate-containing sediment depressurization exploitation numerical simulation method, medium, equipment and product, and relates to the technical field of natural gas hydrate exploitation numerical simulation, the method comprises the following steps: gridding a hydrate-containing sediment depressurization exploitation area based on a voronoi grid, constructing a CVT grid, and carrying out boundary treatment; establishing a control equation under the natural gas hydrate depressurization exploitation condition; establishing the relationship between the density of each phase and the temperature and pressure; and substituting the relationship between the density of each phase and the temperature and the pressure into a control equation, establishing a numerical model by using a finite volume theory according to the divided CVT grid, carrying out discretization, and solving the control equation. According to the method, discrete errors caused by traditional grid centroid deviation are effectively improved, and the calculation efficiency and the solving stability are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for natural gas hydrate extraction, and in particular to numerical simulation methods, media, equipment, and products for depressurization extraction of hydrate-bearing sediments. Background Technology

[0002] Natural gas hydrates, as a clean fossil energy source of significant strategic importance, involve complex phase change dynamics, multiphase flow, and thermodynamic coupling, among other multi-physical field synergies, and have become a cutting-edge topic in the field of energy development research.

[0003] Due to the high cost and long cycle of field tests for natural gas hydrates, indoor simulation experiments and numerical simulations have become the main means to study the exploitation response of hydrate reservoirs. Numerical simulation, as a core technology for studying the exploitation laws of marine natural gas hydrates, hinges on establishing an accurate multi-field coupled mathematical model and achieving precise and efficient solutions. Current mainstream simulators (such as TOUGH+HYDRATE) employ a fully coupled framework and the Voronoi mesh finite volume method. While offering significant advantages in computational accuracy, their core algorithms are limited by commercial software dependence and technological barriers, facing the risk of technological blockade. Furthermore, existing mesh generation techniques generally lack support for centroidal Voronoi meshes (CVT). The centroidal deviation of traditional Voronoi meshes leads to discretization errors, while CVT meshes possess unique geometric advantages in improving the computational accuracy of the finite volume method—aligning the control volume vertices with the centroids of the generation points. Summary of the Invention

[0004] The purpose of this invention is to address the problems of traditional grid centroid deviation, simplified numerical simulation models, and weak mechanical coupling by proposing a numerical simulation method for depressurization mining of hydrate-bearing sediments, comprising the following steps: S1. Based on the Voronoi grid, the depressurization mining area of ​​hydrate-bearing sediments is gridded, a CVT grid is constructed, and boundary processing is performed; S2. Establish the governing equations under depressurization extraction conditions for natural gas hydrates; establish the relationship between the density of each phase and temperature and pressure; S3. Substitute the relationship between the density of each phase and temperature and pressure into the control equation. Based on the divided CVT mesh, establish a numerical model using finite volume theory and discretize it to solve the control equation.

[0005] Furthermore, the specific steps for constructing the CVT mesh are as follows: By abstracting the construction of CVT meshes into solving an energy functional problem, and based on the theory of energy functional minimization, this module formalizes the construction problem of centroid Voronoi mosaics into a nonlinear optimization problem:

[0006] in, express, For seeds, For the mesh to be partitioned, Define the initial region. For planar regions, For Voronoi tessellation units, Let be the density function. y is the seed point; The Lloyd's algorithm is used to solve the nonlinear optimization problem, and through continuous iteration, the following is achieved:

[0007] The solution process is as follows: (1) Initialization: Define the initial distribution of generated points within the computational domain; (2) Voronoi mosaicking: Constructing a Voronoi graph of the generated points; (3) Centroid calculation: Calculate the centroid of each Voronoi element; (4) Point position update: Move each generated point to the centroid of its corresponding Voronoi element; (5) Termination condition: Iterate through (2)–(4) until the convergence criterion is met; A three-dimensional CVT is generated using an axial layer discretization method.

[0008] Furthermore, the boundary processing specifically involves: Whether a seed point is within a region is implicitly defined using the following formula:

[0009] in, This represents the distance between the seed point and the region boundary. The signed distance function is represented as follows: -1 indicates that the seed point is within the region, and +1 indicates that it is not within the region. This represents the region boundary, where y is the seed point and m is the boundary point. For the mesh to be partitioned, Define the initial region. It is a planar region; For seed points near the boundary, a boundary seed point set is formed by performing a boundary reflection operation:

[0010] in, Seed points are located near the boundary points. This represents a distance function that is proportional to the cell size. , For the boundary seed point set, Represents the gradient operator; CVT meshes are formed by seed points within the region.

[0011] Furthermore, the governing equations include: Hydrate mass conservation equation:

[0012]

[0013]

[0014] in, Indicates porosity; , and These represent the densities of gases, water, and hydrates, respectively. , and These represent the saturation levels of the gas phase, aqueous phase, and hydrate phase, respectively. , and These are Darcy velocities relative to the representative unit cell for gas, water, and hydrate, respectively. Let represent the rates of gas generation, water generation, and hydrate dissociation within the control body, respectively, and t represent time. Represents the gradient operator; Hydrate reaction kinetic equation:

[0015] in, This represents the rate at which gas is produced by the decomposition of hydrates per unit time. Represents the reaction rate constant. R represents the activation energy, R represents the ideal gas constant, and T represents the temperature. Indicates the molar mass of the gas. Indicates the specific surface area of ​​the reaction. This represents the phase equilibrium pressure of natural gas hydrates. Indicates gas pressure; Energy conservation equation for hydrates:

[0016] in, Indicates the total heat capacity. , and These represent the relative velocities of the skeleton for gas, water, and hydrate, respectively. Indicates the total thermal conductivity. This represents the latent heat of phase transition during the dissociation of hydrates. , , and These represent the specific heat capacities of the solid matrix, gas, water, and hydrate, respectively. , , and These represent the thermal conductivity of the solid matrix, gas, water, and hydrate, respectively.

[0017] Furthermore, the relationships between the density of each phase and temperature and pressure include: the relationships between the density of gases, water, and hydrates and temperature and pressure, respectively expressed as:

[0018]

[0019]

[0020] in, , and These represent the densities of gases, water, and hydrates, respectively. , and These represent the pressures of gases, water, and hydrates, respectively. Z represents the molar mass of the gas; R represents the correction factor; T represents the ideal gas constant; and T represents the temperature. This indicates the temperature difference from the standard temperature. Indicates the compression factor; Indicates the coefficient of thermal expansion; , , and Indicates a constant coefficient; , The molar mass of the hydrate is Indicates the number of hydrates.

[0021] Furthermore, substituting the obtained relationship between density and temperature and pressure into the governing equations and using finite volume discretization, we obtain the following equation:

[0022]

[0023] Where E represents the unit integration domain; P represents pressure; Indicates the direction of the face; and These represent auxiliary variables related to gas and water, respectively.

[0024] Furthermore, a sequential solution algorithm and Newton's method are used to solve the governing equations, dividing the solutions for temperature and pressure, and saturation, into two subsystems: (1) Update the pressure field in the (k+1)th Newton iteration. and temperature field ,in , The final pressure solution at time n+1 is obtained after iterative convergence. and temperature solution ; (2) Determine the saturation of the hydrate, liquid phase and gas phase at time n+1. , , .

[0025] The present invention also proposes a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described numerical simulation method for depressurization mining of hydrate-bearing sediments.

[0026] The present invention also proposes an electronic device, including a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program including computer-readable instructions, and the processor is configured to invoke the computer-readable instructions to execute the above-described numerical simulation method for depressurization mining of hydrate-bearing sediments.

[0027] The present invention also proposes a computer program product, including a computer program / instruction, which, when executed by a processor, implements the steps of the numerical simulation method for depressurization mining of hydrate-bearing sediments described above.

[0028] The beneficial effects of the technical solution provided by this invention are: This invention first constructs a CVT mesh for depressurization extraction of natural gas hydrates and performs boundary treatment to effectively improve the discretization error caused by centroid deviation of traditional meshes. By conserving hydrate mass, hydrate reaction kinetics, and hydrate energy, the governing equations under depressurization extraction conditions of natural gas hydrates are established. The relationship between the density of each phase of gas, water, and hydrate and temperature and pressure is established. A numerical model is established using finite volume theory and discretized. The governing equations are solved using a sequential solution algorithm and Newton's method, which significantly improves the computational efficiency and solution stability. Attached Figure Description

[0029] Figure 1This is a flowchart of a numerical simulation method for depressurization mining of hydrate-bearing sediments according to an embodiment of the present invention; Figure 2 This is a schematic diagram of a fixed-point iterative algorithm used for CVT mesh generation in an embodiment of the present invention; Figure 3 This is a schematic diagram illustrating the boundary reflection operation for seed points near the boundary points in an embodiment of the present invention. Figure 4 CVT mesh modeling for the indoor experimental area of ​​hydrate-bearing sediment depressurization mining in Comparative Experiment 1 of this invention; Figure 5 A comparison of gas production over time in the Masuda experiment and the GHCUGYSN simulation of sandstone hydrate decomposition experiment. Figure 6 A comparison of pressure versus time results from the Masuda experiment and the GHCUGYSN simulation of sandstone hydrate decomposition. Figure 7 Monitoring points at different locations for the results of Masuda experiment and GHCUGYSN simulation of sandstone hydrate decomposition experiment ( , and A comparison graph of temperature changes over time; Figure 8 This invention provides a regional modeling for the indoor experiment of depressurization mining of hydrate-bearing sediments in Comparative Experiment 2 of this embodiment. Figure 9 A comparison of gas production over time between a Tough+Hydrate depressurization extraction comparative experiment and a GHCUGYSN simulated sandstone hydrate decomposition experiment; Figure 10 A comparison of pressure changes over time at monitoring points (A, B, and C) at different locations, based on the results of a Tough+Hydrate depressurization mining comparative experiment and a GHCUGYSN simulated sandstone hydrate decomposition experiment. Figure 11 A comparison of temperature changes over time at monitoring points (A, B, and C) at different locations, based on the results of a Tough+Hydrate depressurization mining comparative experiment and a GHCUGYSN simulated sandstone hydrate decomposition experiment. Figure 12 This is a block diagram of an electronic device in an exemplary embodiment. Detailed Implementation

[0030] To make the objectives, technical solutions, and advantages of the present invention clearer, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0031] The flowchart of the numerical simulation method for depressurization mining of hydrate-bearing sediments according to an embodiment of the present invention is as follows: Figure 1 Specifically, it includes the following steps: S1. Based on the Voronoi grid, the depressurization mining area of ​​hydrate-bearing sediments is gridded, a CVT grid is constructed, and boundary processing is performed.

[0032] The specific steps for constructing the CVT mesh are as follows: The number of meshes generated is determined in the solution domain. By abstracting the construction of CVT meshes into solving an energy functional problem, and based on the energy functional minimization theory, this module formalizes the construction problem of centroid Voronoi mosaics into a nonlinear optimization problem:

[0033] in, express, For seeds, For the mesh to be partitioned, Define the initial region. For planar regions, For Voronoi tessellation units, Let be the density function. y and y are seed points.

[0034] The functional reaches its minimum if and only if the centroid coincides with the spatial coordinates of the seed point. To solve this optimization problem, a fixed-point iterative algorithm for CVT mesh generation is illustrated in the diagram below. Figure 2 As shown, Lloyd's algorithm is used to solve the nonlinear optimization problem, and through continuous iteration, the following is achieved:

[0035] The solution process is as follows: (1) Initialization: Define the initial distribution of generated points within the computational domain; (2) Voronoi mosaicking: Constructing a Voronoi graph of the generated points; (3) Centroid calculation: Calculate the centroid of each Voronoi element; (4) Point position update: Move each generated point to the centroid of its corresponding Voronoi element; (5) Termination condition: Iterate through (2)–(4) until the convergence criterion is met; For the generation of 3D CVT, an axial layer discretization method is adopted, and parametric stretching along the Z-axis is used to achieve 3D spatial partitioning, which meets the requirements for generating meshes in various dimensions.

[0036] The specific steps for boundary processing are as follows: Whether a seed point is within a region is implicitly defined using the following formula:

[0037] in, This represents the distance between the seed point and the region boundary. The signed distance function is represented as follows: -1 indicates that the seed point is within the region, and +1 indicates that it is not within the region. This represents the region boundary, where y is the seed point and m is the boundary point. For the mesh to be partitioned, Define the initial region. It is a planar region.

[0038] For seed points near the boundary, the boundary reflection operation is illustrated in the diagram below. Figure 3 As shown, the boundary seed point set is formed:

[0039] in, Seed points are located near the boundary points. This represents a distance function that is proportional to the cell size. , For the boundary seed point set, Represents the gradient operator; Finally, only the CVT mesh formed by the seed points within the region is taken, which can automatically trim the boundary. A virtual point set is generated by normal mirror reflection, and after constructing an extended Voronoi diagram, the cells formed by the internal seeds are selected to achieve the boundary fitting solution; for complex boundaries, the boundary is approximated by multi-boundary composite reflection, which effectively ensures the geometric fidelity of complex smooth boundaries.

[0040] S2. Establish the governing equations for natural gas hydrate depressurization extraction under certain conditions; establish the relationship between the density of each phase and temperature and pressure.

[0041] The governing equations include: (1) Mass conservation equation for hydrates:

[0042]

[0043]

[0044] in, Porosity (dimensionless); , and These represent the densities of gases, water, and hydrates, respectively. , and These represent the saturation (dimensionless) of the gas phase, aqueous phase, and hydrate phase, respectively. , and These are Darcy velocities relative to the representative unit cell for gas, water, and hydrate, respectively. These represent the rates of gas generation, water generation, and hydrate dissociation within the control body, respectively, in kg / (m³). 3 •s), where t represents time. Represents the gradient operator; (2) Reaction kinetic equation of hydrate:

[0045] in, This represents the rate at which gas is produced by the decomposition of hydrates per unit time. This represents the reaction rate constant, with units of mol / (m²). 2 •Pa•s), The value represents the activation energy, expressed in J / (mol•K); R represents the ideal gas constant, expressed in J / (mol•K); and T represents the temperature, expressed in K. This indicates the molar mass of a gas, expressed in kg / mol. The specific surface area of ​​the reaction is expressed in units of 1 / m². This represents the phase equilibrium pressure of natural gas hydrates, expressed in Pa. Indicates gas pressure; (3) Energy conservation equation for hydrates:

[0046] in, Indicates the total heat capacity. , and These represent the relative velocities of the skeleton for gas, water, and hydrate, respectively. Indicates the total thermal conductivity. The latent heat of phase transition during the dissociation of hydrates is expressed in J / (kg•s). , , and The values ​​represent the specific heat capacities of the solid matrix (s), gas (g), water (w), and hydrate (h), respectively, in J / (kg•K). , , and These represent the thermal conductivity of the solid matrix, gas, water, and hydrate, respectively, in W / (m•K).

[0047] The relationships between the density of each phase and temperature and pressure include: (1) The relationship between gas density and temperature and pressure

[0048] (2) The relationship between water density and temperature and pressure

[0049] (3) The relationship between the density of hydrates and temperature and pressure

[0050] in, , and These represent the densities of gases, water, and hydrates, respectively. , and These represent the pressures of gases, water, and hydrates, respectively. Z represents the molar mass of the gas; Z represents the correction factor (dimensionless); R represents the ideal gas constant (dimensionless); T represents the temperature. This indicates the temperature difference from the standard temperature. This represents the compressibility factor, with units of Pa. -1 ; This represents the coefficient of thermal expansion, with units of K. -1 ; , , and Indicates a constant coefficient. , , , ; (dimensionless) The molar mass of the hydrate is expressed in kg / mol. Indicates the number of hydrates.

[0051] S3. Substitute the relationship between the density of each phase and temperature and pressure into the control equation. Based on the divided CVT mesh, establish a numerical model using finite volume theory and discretize it to solve the control equation.

[0052] Substituting the relationship between density and temperature and pressure into the governing equations and using finite volume discretization, we obtain the following equation:

[0053]

[0054] Where E represents the element integration domain; P represents pressure; and n represents the surface direction. and These represent auxiliary variables related to gas and water, respectively.

[0055] A sequential solution algorithm and Newton's method are used to systematically decouple the strongly nonlinear control equations, while taking into account both numerical stability and physical consistency. The solutions for temperature and pressure, and saturation, are divided into two subsystems, and the thermo-pressure variables and saturation are solved hierarchically.

[0056] (1) Update the pressure field in the (k+1)th Newton iteration. and temperature field ,in , The final pressure solution at time n+1 is obtained after iterative convergence. and temperature solution ; (2) Determine the saturation of the hydrate, liquid phase and gas phase at time n+1. , , .

[0057] The results are visualized by integrating a spatiotemporal visualization engine, enabling dynamic analysis of key processes such as phase change interface evolution, gas production rate, and dynamic evolution of temperature and pressure. This invention, through the above process, develops a numerical simulator for natural gas hydrate depressurization extraction (GHCUGYSN), including a preprocessing module, a solver module, and a post-processing module. The GHCUGYSN preprocessing module generates a CVT mesh and performs boundary treatment; the GHCUGYSN solver module establishes a numerical model using finite volume theory and discretizes it, then solves it using a sequential solution method; the GHCUGYSN post-processing module visualizes the data results. This achieves an integrated workflow from mesh generation to numerical solution to result visualization. The technical solution of this invention will be further described in detail below through experiments.

[0058] Example 1: Comparison of results from the Masuda indoor sandstone depressurization extraction experiment (a landmark classic experiment in the field of natural gas hydrate research): This invention provides an embodiment of a numerical simulation experiment for indoor sandstone depressurization mining in Masuda. The specific steps are as follows: The preprocessing module of GHCUGYSN is used to model the area of ​​the indoor depressurization mining experiment for hydrate-bearing sediments. The CVT mesh modeling of the indoor depressurization mining area for hydrate-bearing sediments in Experiment 1 is compared with that in Experiment 1. Figure 4 As shown, the geometric parameters of the model are shown in Table 1.

[0059] Table 1

[0060] The meshed model was solved using the GHCUGYSN solver module, and experimental parameters were assigned. Specific parameter values ​​are shown in Table 2. Furthermore, an isothermal and isobaric boundary was used at the left end of the model, while no other seepage boundaries were observed.

[0061] Table 2

[0062] Based on this model, this example compares and simulates the system temperature, pore pressure, water vapor production rate, hydrate saturation, and displacement evolution under a 5-hour indoor experiment in Masuda sandstone. The results of the Masuda experiment and the GHCUGYSN simulation of sandstone hydrate decomposition are as follows: Figures 5-7 As shown, where, Figure 5 This is a comparison graph showing the change in gas production over time between the Masuda experiment and the GHCUGYSN simulation. Figure 6 This is a comparison graph of pressure changes over time in the Masuda experiment and the GHCUGYSN simulation. Figure 7 These are monitoring points at different locations, as simulated by the Masuda experiment and GHCUGYSN. , and A comparison chart of temperature changes over time.

[0063] A comparison with the results of Masuda sandstone hydrate decomposition experiments shows that the GHCUGYSN simulator can effectively capture the evolution of key multiphase flow characteristics such as gas production, pressure, and temperature propagation. However, there are deviations in the reconstruction of the far-end pressure field and the characterization of local temperature monitoring point data. The possible reasons are as follows: (1) Missing salt effect coupling: The phase equilibrium condition shift caused by the salinity (1%) of the pore fluid in the experimental system has not been fully incorporated into the current version of the phase equation, resulting in a deviation in the prediction of the hydrate decomposition interface.

[0064] (2) Non-uniform heat transfer mechanism: Abnormal fluctuations in the pressure field at the far end (first decrease and then increase) may be related to the lack of boundary condition parameters for heat conduction of the plug in the experimental device. Overestimation of the heat flow diffusion rate in the numerical model may cause a non-physical "pre-decomposition" phenomenon, thus causing such errors.

[0065] Example 2: Comparative Experiment of Tough + Hydrate Decompression Mining: The preprocessing module of GHCUGYSN was used to model the area of ​​the indoor experiment on depressurization mining of hydrate-bearing sediments, and the modeling was compared with the area modeling reference of the indoor experiment on depressurization mining of hydrate-bearing sediments in Experiment 2. Figure 6 The geometric parameters of the model are shown in Table 3.

[0066] Table 3

[0067] The meshed model was solved using the GHCUGYSN solver module, and experimental parameters were assigned. Specific parameter values ​​are shown in Table 4. Furthermore, an isothermal and isobaric boundary was used at the left end of the model, while no other seepage boundaries were observed.

[0068] Table 4

[0069] Figures 9-11 These are the results of a comparative experiment on Tough+Hydrate depressurization mining and a simulated sandstone hydrate decomposition experiment using GHCUGYSN. Figure 9 This is a comparison graph showing the change in gas production over time between Tough+Hydrate and GHCUGYSN simulations. Figure 10 This is a comparison graph showing the pressure changes over time at monitoring points (A, B, and C) at different locations simulated by Tough+Hydrate and GHCUGYSN. Figure 11 This is a comparison graph showing the temperature changes over time at monitoring points (A, B, and C) at different locations simulated by Tough+Hydrate and GHCUGYSN.

[0070] By comparing the prediction results with those of the mainstream simulator TOUGH+HYDRATE using the same initial conditions, it was found that the GHCUGYSN of this invention exhibits excellent multi-field coupling calculation performance, and its cumulative gas production and the evolution of physical property parameters at key monitoring points are in high agreement with its calculation results.

[0071] This invention has the following significant advantages: (1) In terms of mesh generation, based on the adaptive algorithm of reflection boundary and various optimizer solution strategies, the bottleneck of traditional CVT mesh generation under complex geometric models has been broken. (2) Through Masuda experimental calibration and cross-validation with TOUGH+HYDRATE software, a three-field coupled solution architecture covering hydrate phase change-seepage-heat transfer was constructed, providing an important foundation for large-scale site-level multiphysics simulation. The modular design of this platform lays the algorithmic foundation for subsequent extensions of cutting-edge models such as heat flow solidification coupling.

[0072] In one exemplary embodiment, a computer-readable storage medium is included, which stores a computer program that, when executed by a processor, implements the numerical simulation method for depressurization mining of hydrate-bearing sediments described above.

[0073] Please see Figure 12 In one exemplary embodiment, the device further includes an electronic device including at least one processor, at least one memory, and at least one communication bus.

[0074] The memory stores a computer program, which includes computer-readable instructions. The processor calls the computer-readable instructions stored in the memory through the communication bus to execute the numerical simulation method for depressurization mining of hydrate-bearing sediments.

[0075] In one exemplary embodiment, a computer program product is proposed, comprising a computer program / instructions that, when executed by a processor, implement the steps of the numerical simulation method for depressurization mining of hydrate-bearing sediments described above.

[0076] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A numerical simulation method for gas production from a hydrate deposit under pressure reduction, characterized in that, The method comprises the following steps: S1, gridding a hydrate deposit depressurization area based on a voronoi grid, constructing a CVT grid, and performing boundary processing; S2, establishing a control equation under a natural gas hydrate depressurization condition, and a relationship between the density of each phase and temperature and pressure; S3, substituting the relationship between the density of each phase and temperature and pressure into the control equation, establishing a numerical model and performing discretization according to the divided CVT grid by using a finite volume theory, and solving the control equation.

2. The numerical simulation method for pressure depletion of a hydrate deposit according to claim 1, wherein, The construction of the CVT grid is specifically as follows: By abstracting the construction of the CVT grid as solving an energy functional problem, based on the energy functional minimization theory, the module formulates the construction problem of the centroid Voronoi tessellation into a nonlinear optimization problem: wherein, denotes, is a seed, is a grid domain to be subdivided, is an initially defined region, is a planar region, is a Voronoi tessellation cell, is a density function, and y is a seed point; Lloyd's algorithm is used to solve the nonlinear optimization problem, and by constantly iterating: The solving process is as follows: (1) initialization: defining an initial generated point distribution in a calculation domain; (2) Voronoi tessellation: constructing a Voronoi diagram of the generated points; (3) centroid calculation: calculating the centroid of each Voronoi unit; (4) point position updating: moving each generated point to the centroid position of its corresponding Voronoi unit; (5) termination condition: iteratively performing (2)-(4) until the convergence criterion is met; A three-dimensional CVT is generated by using an axial layering discretization method.

3. The numerical simulation method for pressure depletion of a hydrate deposit according to claim 1, wherein, The boundary processing is specifically as follows: Whether a seed point is in the region is implicitly defined by the following formula: wherein, denotes the distance of a seed point to the region boundary, denotes the signed distance function: -1 for a seed point inside the region, +1 for outside the region, denotes the region boundary, y is a seed point, m denotes a boundary point, is a grid to be partitioned domain, is an initial defined region, is a planar region; For seed points close to the boundary, a boundary seed point set is formed by reflection operation with respect to the boundary: wherein, is a seed point close to the boundary point, denotes a distance function proportional to the cell size, , is a set of boundary seed points, denotes a gradient operator; The CVT grid formed by the seed points in the region is taken.

4. The numerical simulation method for pressure depletion of a hydrate deposit according to claim 1, wherein, The control equation comprises: A hydrate mass conservation equation: wherein, denotes porosity; , and denote the densities of gas, water and hydrate, respectively; , and are the saturations of gas, water and hydrate phases, respectively; , and are the Darcy velocities of gas, water and hydrate with respect to a representative elementary volume, respectively; and denote the control of the gas generation rate, the water generation rate and the hydrate dissociation rate in the body, respectively, t denotes time, denotes the gradient operator; A hydrate reaction kinetics equation: wherein, represents a rate of gas production per unit time by hydrate dissociation, represents a reaction rate constant, represents an activation energy, R represents an ideal gas constant, and T represents a temperature, represents a molar mass of a gas, represents a specific surface area of a reaction, represents a phase equilibrium pressure of a natural gas hydrate, represents a gas pressure; A hydrate energy conservation equation: wherein, represents the total heat capacity, , and represent the relative skeleton velocities of the gas, water and hydrate, respectively, represents the total thermal conductivity, represents the latent heat of phase transition during hydrate dissociation, , , and represent the specific heat capacities of the solid matrix, gas, water and hydrate, respectively, , , and represent the thermal conductivities of the solid matrix, gas, water and hydrate, respectively.

5. A numerical simulation method for pressure depletion of a hydrate deposit according to claim 4, wherein, The relationship between the density of each phase and temperature and pressure comprises the relationship between the density of gas, water and hydrate and temperature and pressure, which is respectively represented as: wherein , and represent the density of gas, water and hydrate, respectively; , and represent the pressure of gas, water and hydrate, respectively; represents the molar mass of gas; Z represents a correction factor; R represents the ideal gas constant; T represents the temperature; represents the difference to the temperature at standard conditions; represents the compressibility factor; represents the expansion factor; , , and represent the constant factor; , is the molar mass of the hydrate, represents the number of hydrates.

6. The numerical simulation method for hydrate deposit depressurization according to claim 5, characterized in that, The relationship between the density and temperature and pressure is substituted into the control equation and discretized by using the finite volume to obtain the following equation: where E denotes the element integral domain; P denotes pressure; surface normal direction; and denote gas and water related auxiliary variables, respectively.

7. The numerical simulation method for pressure depletion of a hydrate deposit according to claim 1, wherein, The control equation is solved by using a sequential solving algorithm and Newton's method, and the solving of temperature and pressure and the solving of saturation are divided into two subsystems: (1) Update the pressure field in the (k+1)th Newton iteration and the temperature field where , After the iteration converges, the final pressure solution at time n+1 and temperature solution are obtained; (2) solving the hydrate, liquid phase and gas phase saturation at time n+1 , , .

8. A computer readable storage medium storing a computer program, characterized in that: The computer program is executed by the processor to realize the method of any one of claims 1-7.

9. An electronic device, comprising: The computer program is executed by the processor to realize the method of any one of claims 1-7.

10. A computer program product comprising computer programs / instructions, characterized in that, The computer program / instructions are executed by the processor to realize the steps of the method of any one of claims 1-7.

Images