A method, apparatus, system, and storage medium for developing a numerical simulator.

Abstract

This invention discloses a method, apparatus, system, and storage medium for developing a numerical simulator. The development method includes: using mesh generation software to mesh the target study area and outputting a mesh information file, which is then converted into a node coordinate information matrix and a node number information matrix for the target study area; calculating the pressure, temperature, and three-phase saturation distribution of natural gas hydrates during the extraction process based on a derived incremental analytical formula, while simultaneously calculating the displacement distribution within the area using the Galerkin finite element method; constructing control volumes around the mesh vertices when using the incremental analytical formula; and decoupling between different physical fields in an explicit recursive manner. This method eliminates the need to solve linear equations when solving the seepage heat transfer system using an incremental approach, avoiding the generation of large coefficient matrices and saving memory. This approach also demonstrates good adaptability when considering the elastoplastic characteristics of natural gas hydrate deposits.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for the exploitation of natural gas hydrate resources, and in particular to a method, apparatus, system and storage medium for developing a numerical simulator. Background Technology

[0002] Natural gas hydrate (commonly known as combustible ice) is a clean and efficient energy source. It is widely distributed, has huge reserves, and high energy density, making it a future energy source for humankind.

[0003] The extraction of natural gas hydrates (hereinafter referred to as hydrates) involves multiple physical fields, multiple phases and multiple components. The evolution process of each physical field is very complex. Related experimental research and field experiments are very time-consuming and labor-intensive, and it is difficult to fully demonstrate the evolution process of each physical field.

[0004] Numerical simulation technology, by solving the governing equations, can fully display the distribution cloud map of various field parameters throughout the entire time domain, facilitating the development of reasonable and efficient mining schemes. Therefore, numerical simulation technology is an indispensable, cost-effective, and efficient method for studying this complex evolutionary process.

[0005] Numerical simulation technology uses mathematical modeling to model the problem under study, resulting in a set of governing equations (usually a system of partial differential equations) and algebraic equations, which constitute a differential-algebraic system. Under certain assumptions, this system can describe the evolution of each physical field in terms of time and space with relatively high accuracy. After the system is established, the region under study needs to be meshed, and numerical methods are used to calculate the values ​​of each physical quantity on the mesh, thereby obtaining the spatiotemporal evolution characteristics of each field parameter.

[0006] However, most existing studies are based on commercial software, which tends to generate large coefficient matrices when solving linear equations, consuming a certain amount of memory, and is not well adapted to considering the elastoplastic characteristics of hydrate deposits. Summary of the Invention

[0007] This invention provides a method, apparatus, system, and storage medium for developing a numerical simulator to solve the technical problem of multi-physics coupling encountered in the hydrate mining process in the prior art.

[0008] To address the aforementioned problems, the primary objective of this invention is to provide a method for developing a numerical simulator for multiphysics coupling of hydrate sediments, the method comprising:

[0009] S 100 The target study area is meshed using mesh generation software, and the mesh information file is output.

[0010] S 200Import the grid information file and convert it into a matrix of coordinate information for each node in the target study area, a matrix of node numbers in each grid, and a matrix of grid numbers and node numbers in the boundary area.

[0011] S 300 Based on the initial and boundary conditions of the temperature field, seepage field, mechanical field and phase transition field of the hydrate, the temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation and porosity are initialized.

[0012] S 400 Construct a control volume centered on a grid vertex, calculate the flux of each boundary of the control volume, and accumulate them to obtain the flux difference of each control volume;

[0013] S 500 Substitute the flux difference of each control volume into the incremental formula, and calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation and hydrate saturation of each grid vertex at the next moment based on the field parameters at the current moment.

[0014] S 600 For the mechanical governing equations of the multiphysics coupled mathematical model of the hydrate-bearing sediments, the Galerkin finite element method is used to solve the purely elastic problem, and the incremental variable stiffness method is used to solve the elastoplastic problem.

[0015] S 700 The reservoir porosity is updated based on the calculated volumetric strain, temperature, and pore pressure, and the absolute permeability is updated based on the porosity.

[0016] S 800 Determine whether the calculation results meet the stability conditions. If they do, continue the loop calculation. If they do not meet the stability conditions, reduce the time step and recalculate until the stability conditions are met. Output the calculation results of each physical field at the current time in the form of a field diagram.

[0017] Furthermore, in step S 400 In this context, the calculation of fluxes at each boundary of the control volume specifically includes water flux and gas flux, wherein:

[0018] The flux of water about a certain cross section

[0019] (1)

[0020] Gas flux about a certain cross section

[0021] (2)

[0022] in: It is the permeability of the aqueous phase. It is the permeability of the gas phase. and These are the pore water pressure and pore gas pressure, respectively, assembled from the finite element basis functions and nodal function values. , , These are finite element basis functions. It is the number of grid vertices. It is the unit normal vector of that cross section. This represents the differential of the current cross section.

[0023] Furthermore, in step S 500 The specific expression for the pore gas pressure increment is as follows:

[0024] (3)

[0025] The specific expression for the increment of pore water pressure is as follows:

[0026] (4)

[0027] The specific expression for the gas phase saturation increment is as follows:

[0028] (5)

[0029] The specific expression for the increment of water phase saturation is as follows:

[0030] (6)

[0031] The specific expression for the hydrate saturation increment is as follows:

[0032] (7)

[0033] In the formula: , and ,

[0034] in: ( ) represents the divergence of water flux. ( ) represents the divergence of gas flux. Represents gas phase saturation. Represents the water phase saturation. Represents hydrate saturation. Represents effective water saturation. Represents the stiffness of the gas phase. Represents the stiffness of the water phase. Represents the stiffness of hydrates. It is the gas phase pore pressure. It is the pore pressure of the aqueous phase. It is porosity. It is volumetric strain. It is the time step. It is the molar mass of the gas phase. It is the molar mass of the aqueous phase. It is the molar mass of the hydrate. It is water and numbers. It is the gas phase density. It is the density of the aqueous phase. It is the density of the hydrate. It is the coefficient of thermal expansion of the gas phase. It is the coefficient of thermal expansion of water phase. It is the coefficient of thermal expansion of hydrates. It is the coefficient of thermal expansion of sediment particles. It is the number of moles of hydrate decomposition. It's temperature. It is capillary pressure. It is pore air pressure. It is pore water pressure. It is gas saturation. It's water saturation. It is the increase in hydrate saturation;

[0035] The divergence of the water flux ( The result is obtained from equation (8):

[0036] (8)

[0037] The divergence of the gas flux ( The result is obtained from equation (9):

[0038] (9)

[0039] In the formula: Represents the boundary controlling the volume. Represents the density of water. Represents gas density; Represents the seepage velocity of the aqueous phase. Represents the seepage velocity in the gas phase. It is the area of ​​the current control body. Represents the control body. Represents the boundary of the control volume. It is the unit vector of the outward normal to the boundary of the control volume. The differential represents the boundary of the control volume.

[0040] Furthermore, in step S 500 In the above, the expression for calculating the temperature increment of each grid vertex at the next moment is:

[0041] (10)

[0042] in: It is an increase in temperature. It is the weighted average of the phase heats. It is air compared to heat, Water is hotter than water. It is the specific heat of hydrates. It is the thermal conductivity tensor. It is the enthalpy change caused by the decomposition of hydrates. It is the effective stress. It is plastic strain. ( ) represents the divergence of heat flux with respect to the control volume.

[0043] Furthermore, in step S 600 The method of using the Galerkin finite element method to solve purely elastic problems specifically includes:

[0044] According to step S 500 The calculated hydrate saturation is used to update the elastic constants;

[0045] According to step S 500 The calculated temperature and pore pressure were used to calculate the displacement considering the thermal-fluid-structure interaction using a purely elastic finite element program.

[0046] The stress-strain and volumetric strain were calculated.

[0047] Furthermore, in step S 600 The method of using incremental variable stiffness to solve elastoplastic problems specifically includes:

[0048] Gradually apply the load and determine whether the stress in the mesh reaches the yield state;

[0049] Based on whether the yield state has been reached, the correction coefficient matrix is ​​modified, and the displacement solution of the target study area is finally obtained by iterative process.

[0050] The strain and stress of the target study area were calculated.

[0051] A second objective of this invention is to provide a development apparatus for a numerical simulator, comprising:

[0052] Mesh partitioning cells are used to partition the target study area into meshes and output mesh information files;

[0053] The conversion unit is used to convert the grid information file into a matrix of coordinate information of each node in the target study area, a matrix of node numbers contained in each grid, and a matrix of grid numbers and node numbers in the boundary area.

[0054] The initialization unit is used to initialize temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation, and porosity.

[0055] The flux calculation unit is used to calculate the flux at each boundary of the control volume and accumulate them to obtain the flux difference for each control volume.

[0056] The incremental calculation unit is used to calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation, and hydrate saturation at each grid vertex in the next time step, based on the field parameters at the current time step.

[0057] The solution unit is used to solve the purely elastic problem using the Galerkin finite element method and the elastoplastic problem using the incremental variable stiffness method for the mechanical governing equations of the multiphysics coupled mathematical model of the hydrate-bearing sediments.

[0058] The update unit is used to update the reservoir porosity based on the calculated volumetric strain, temperature and pore pressure, and update the absolute permeability based on the porosity.

[0059] The decision unit is used to determine whether the calculation result meets the stability condition.

[0060] A third objective of this invention is to provide a numerical simulator development system, comprising a lower-level simulator and a higher-level computer that interacts with the lower-level simulator, wherein the lower-level simulator is used to implement the numerical simulator development method described above, and the higher-level computer is used to receive and parse the updated parsed data generated by the lower-level simulator.

[0061] A fourth objective of this invention is to provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the development method of the numerical simulator as described above.

[0062] Compared with the prior art, the present invention has significant advantages and beneficial effects, specifically reflected in the following aspects:

[0063] This invention proposes a simulator construction method based on incremental schemes and the finite element method. This method calculates the pressure, temperature, and three-phase saturation distribution of natural gas hydrates during extraction based on derived incremental analytical formulas. When solving for pore water pressure, pore gas pressure, water phase saturation, gas saturation, hydrate saturation, and temperature across the entire domain, a control volume is constructed centered on the grid vertices. The flux difference with respect to this control volume is calculated, and then the incremental analytical formulas are used to calculate the time increment of each physical quantity at each grid node. Simultaneously, the parameters of the mechanical field are solved using the classical finite element method for both purely elastic and elastoplastic problems. In this way, the basic information calculated by the solver is stored at the grid vertices, avoiding interpolation errors between the grid center data and the grid vertex data. Decoupling between different physical fields is performed in an explicit recursive manner. This method eliminates the need to solve linear equations when using an incremental approach to solve the seepage heat transfer system, thus avoiding the generation of large coefficient matrices and saving memory. It also has good adaptability when considering the elastoplastic characteristics of hydrate deposits. Through this method, the water and gas production and deformation laws of hydrate reservoirs during the mining process can be obtained, and the evolution laws of various physical fields can be understood in depth. Attached Figure Description

[0064] Figure 1 This is a flowchart of the development method of the numerical simulator in this embodiment of the invention;

[0065] Figure 2 This is a schematic diagram of the phase diagram fitting curve of methane hydrate in an embodiment of the present invention;

[0066] Figure 3 This is a comparison diagram of phase equilibrium model curves in embodiments of the present invention;

[0067] Figure 4 This is a schematic diagram of multiphysics coupling of hydrate-containing sediments in an embodiment of the present invention;

[0068] Figure 5 This is a control volume diagram centered on the grid vertices in this embodiment of the invention;

[0069] Figure 6 This is a boundary flux map inside a grid in an embodiment of the present invention;

[0070] Figure 7 This is a schematic diagram of the solution of a multi-field coupling model for depressurization mining of hydrate-bearing sediments in an embodiment of the present invention. Detailed Implementation

[0071] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0072] Please see Figure 1As shown, this embodiment of the invention provides a method for developing a numerical simulator for multiphysics coupling of hydrate sediments. The development method includes:

[0073] S 100 The target study area is meshed using mesh generation software, and the mesh information file is output.

[0074] In this step, the mesh generation software used is, for example, the open-source software Gmsh, to generate meshes for the target study area, dividing it into general triangular, quadrilateral, tetrahedral, and hexahedral meshes. It is important to note that triangular and quadrilateral meshes correspond to two-dimensional regions, while tetrahedral and hexahedral meshes correspond to three-dimensional regions. This allows for the generation of appropriate meshes based on the structural characteristics of the actual target study area.

[0075] It should be noted that by using existing preprocessing software such as Hypermesh, Gambit, and Gmsh for mesh generation, the simulator only needs to convert the mesh file into the data required for computation. Similarly, by using existing post-processing software such as Tecplot and Paraview for post-processing, the simulator only needs to output the corresponding files.

[0076] S 200 Import the grid information file and convert it into a matrix of coordinate information for each node in the target study area, a matrix of node numbers in each grid, and a matrix of grid numbers and node numbers in the boundary area.

[0077] Specifically, in this step, after the target study area is divided into grids, the grids are divided by the connection of nodes. Each node is assigned a number and coordinate information. In this way, the grids in the entire target study area can be converted into a node coordinate information matrix, an information matrix of the node numbers contained in each grid, and a grid number and node number information matrix of the grids in the boundary area.

[0078] S 300 Based on the initial and boundary conditions of the temperature field, seepage field, mechanical field and phase transition field of the hydrate, the temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation and porosity are initialized.

[0079] A multiphysics coupled mathematical model was established for hydrate-bearing sediments. The model establishment process is as follows:

[0080] (1) Conservation of mass:

[0081] During the decomposition of hydrates, the gas, water, and hydrate phases must satisfy the mass conservation equation during their evolution.

[0082] (11)

[0083] (12)

[0084] (13)

[0085] Of the three formulas above, Water encapsulated within hydrates Represents the gas encapsulated within the hydrate, and has ; The vector representing the flow of water. A vector representing the flow of gas; Represents the rate of gas production from hydrate decomposition. This represents the rate at which hydrates decompose and produce water.

[0086] (2) Energy conservation equation

[0087] (14)

[0088] In the above formula, It is the effective specific heat, determined by the following formula:

[0089] (15)

[0090] For an effective thermal conductivity tensor, the volume average of different components can be taken; It is the enthalpy change during the endothermic decomposition of hydrates. ); The strain energy consumption is due to plastic strain. In this model, it is assumed that all the work done due to plastic strain is converted into heat energy.

[0091] (3) Changes in pore volume

[0092] Meanwhile, the changes in pore volume are as follows:

[0093] (16)

[0094] in It is volumetric strain.

[0095] The change in pore volume can be expressed as the sum of volumetric strain and temperature strain. Here it is assumed that the solid (hydrates and sand particles) will not be compressed or expanded, but the soil skeleton composed of particles can still be compressed.

[0096] (4) Kinetic equations for hydrate decomposition

[0097] The decomposition and regeneration of hydrates are governed by first-order kinetics, meaning the phase transition rate is proportional to the product of the specific surface area and the driving force (the difference between ambient pressure and phase equilibrium pressure), as shown in the following equation:

[0098] (17)

[0099] (18)

[0100] (19)

[0101] in: Let be the hydrate decomposition constant, and Calculated by the following formula: , ; It is the specific surface area of ​​hydrates; This is the value of the hydrate formation constant.

[0102] (5) The following formula can be used to calculate water and air flux:

[0103] Van Genuchten (1980) unsaturated flow model

[0104] (20)

[0105] (twenty one)

[0106] (twenty two)

[0107] (twenty three)

[0108] Capillary pressure:

[0109] (twenty four)

[0110] (25)

[0111] in: , , , They are all coefficients. This is called effective water saturation, and , It is capillary pressure. It is the absolute permeability without hydrates. It is the absolute permeability of the hydrate.

[0112] (6) Hydrate phase transition

[0113] To describe the reaction behavior of hydrate deposits in a model, it is essential to first determine when the hydrates decompose and when they stabilize. Since most natural gas hydrates in nature are methane hydrates, the phase equilibrium conditions of methane hydrates are the primary focus of this study.

[0114] In previous studies, different scholars have proposed different phase equilibrium models for methane hydrate. (Kamath, 1987) proposed the following phase equilibrium model:

[0115] (26)

[0116] in: It is the phase equilibrium pressure, which is the critical pressure at which natural gas hydrates decompose.

[0117] like Figure 2 , 3 As shown, Moridis et al. presented a phase diagram of methane hydrate and a comparison of two phase equilibrium models, which showed that there are certain differences between the two models. Here, we prioritize Moridis' phase equilibrium model.

[0118] Decomposition kinetic model

[0119] The most commonly used kinetic model for hydrate decomposition is the Kim model. In 1987, Kim's research found that hydrate decomposition and regeneration are governed by first-order kinetics, meaning the phase transition rate is proportional to the product of the specific surface area and the driving force (the difference between ambient pressure and phase equilibrium pressure). The specific phase transition rate is shown in the following equation:

[0120] (27)

[0121] (28)

[0122] (29)

[0123] in: It is the hydrate decomposition constant, and Calculated by the following formula:

[0124] (30)

[0125] in, ; It is the specific surface area of ​​the hydrate, which can be taken as 0.375 for a unit volume of hydrate; It is the hydrate formation constant. .

[0126] (7) Mechanical model

[0127] In the theory of elasticity, the equilibrium equations can be expressed in tensor form as follows:

[0128] (31)

[0129] Expanding the above equation yields:

[0130] (32)

[0131] The tensor form of the geometric equations is then expressed by the following equation:

[0132] (33)

[0133] In indicator form:

[0134] (34)

[0135] The commonly used form in engineering is:

[0136] (35)

[0137] Effective stress formula

[0138] Bishop's formula, derived from Terzaghi's effective stress principle:

[0139] (36)

[0140] mechanical constitutive equations

[0141] The elasticity model can be expanded into the following matrix form:

[0142] (37)

[0143] in:

[0144] (38)

[0145] Many scholars have conducted experiments to explore the relationship between hydrate saturation and elastic modulus E and Poisson's ratio v. The main view is that the elastic modulus increases with the increase of hydrate saturation, while Poisson's ratio is less affected by hydrate saturation.

[0146] One relatively simple form is the one obtained by Ebinuma et al. through triaxial experimental fitting:

[0147] (39)

[0148] Liu Lele et al. proposed the following linear model:

[0149] (40)

[0150] In the above formula: This is the Young's modulus without hydrates. The saturation of hydrate is Young's modulus.

[0151] Please see Figure 4 As shown, the governing equations and algebraic formulas described above describe the characteristics of each physical field, and there are relationships between the physical fields as follows: Figure 4 The coupling relationship is shown.

[0152] S 400 Construct a control volume centered on a grid vertex, calculate the flux of each boundary of the control volume, and accumulate them to obtain the flux difference of each control volume.

[0153] Specifically, in the embodiments of the present invention, in step S 400 In this context, the fluxes at each boundary of the computational control volume specifically include water flux and gas flux.

[0154] in:

[0155] The flux of water about a certain cross section

[0156] (1)

[0157] Gas flux about a certain cross section

[0158] (2)

[0159] in: It is the permeability of the aqueous phase. It is the permeability of the gas phase. It is the unit normal vector of that cross section. and These are the pore water pressure and pore gas pressure, respectively, assembled from the finite element basis functions and nodal function values. , , These are finite element basis functions. These are the pore air pressure and pore water pressure at the grid vertices, respectively.

[0160] Please see Figure 5 As shown, unlike the traditional finite volume method, this embodiment of the invention borrows the idea of ​​the control volume finite element method. Instead of establishing a control volume for the mesh, it establishes a control volume centered on the mesh vertices, and then performs calculations using incremental formulas. The advantage of this approach is that all solution information is concentrated at the mesh vertices, facilitating coupled solution development.

[0161] S500 Substitute the flux difference of each control volume into the incremental formula, and calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation, and hydrate saturation of each grid vertex at the next moment based on the field parameters at the current moment.

[0162] For solving the multi-field coupling model of hydrates, due to the complexity of the problem, the current decoupling method is sequential decoupling, that is, first solve one physical field, and then solve the next physical field based on the updated physical field, and continue to cycle until the maximum time step is reached.

[0163] The solution of the above physical field can be divided into the solution of the seepage subsystem and the geomechanical subsystem, that is, the incremental analytical formula is used to calculate the increment of pore gas pressure, the increment of pore water pressure, the increment of gas phase saturation, the increment of water phase saturation and the increment of hydrate saturation.

[0164] The specific expression for the pore gas pressure increment is as follows:

[0165] (3)

[0166] The specific expression for the increment of pore water pressure is as follows:

[0167] (4)

[0168] The specific expression for the gas phase saturation increment is as follows:

[0169] (5)

[0170] The specific expression for the increment of water phase saturation is as follows:

[0171] (6)

[0172] The specific expression for the hydrate saturation increment is as follows:

[0173] (7)

[0174] As can be seen, the above formula consists of four parts.

[0175] The first part represents the effect of fluid flow, the second part represents the effect of mechanical deformation, the third part represents the effect of hydrate decomposition or regeneration, and the fourth part represents the thermal effect.

[0176] In this way, the calculation of each physical variable will take into account the influence of other physical fields, and naturally the effect of coupling will be taken into account.

[0177] in: , and ,

[0178] in: The divergence representing water flux The divergence representing gas flux Represents gas phase saturation. Represents the water phase saturation. Represents hydrate saturation. Represents effective water saturation. , Represents the stiffness of the gas phase. Represents the stiffness of the water phase. Represents the stiffness of hydrates. It is the gas phase pore pressure. It is the pore pressure of the aqueous phase. It is porosity. It is volumetric strain. It is the time step. It is the molar mass of the gas phase. It is the molar mass of the aqueous phase. It is the molar mass of the hydrate. It is water and numbers. It is the gas phase density. It is the density of the aqueous phase. It is the density of the hydrate. It is the coefficient of thermal expansion of the gas phase. It is the coefficient of thermal expansion of water phase. It is the coefficient of thermal expansion of hydrates. It is the coefficient of thermal expansion of sediment particles. It is the number of moles of hydrate decomposition. It's temperature. It is capillary pressure. , It is pore air pressure. It is pore water pressure. It is gas saturation. It's water saturation. It is the increase in hydrate saturation.

[0179] The calculation of the above incremental analytical formula requires prior calculation of the divergence of the water flux. and the divergence of gas flux The calculations for these two items are performed according to the following formula:

[0180] The divergence of the water flux Calculated from equation (8):

[0181] (8)

[0182] The divergence of the gas flux Calculated from equation (9):

[0183] (9)

[0184] In the formula: The number representing the component that controls the volume in each surrounding grid. Represents the density of water. Represents gas density; Represents the seepage velocity of the aqueous phase. Represents the seepage velocity in the gas phase. It is the area of ​​the current control body. Represents the control body. Represents the boundary of the control volume. It is the unit vector of the outward normal to the boundary of the control volume. The differential represents the control volume boundary.

[0185] This represents the portion of the control volume to which each grid cell belongs. It is the unit vector of the outward normal of the control volume boundary in this part of the mesh. This represents the portion of the water in each grid cell that belongs to the control volume. This represents the portion of the gas in each grid that belongs to the control volume. Represents the density of water. Represents the seepage velocity of the aqueous phase. Represents the density of water. Represents the seepage velocity of the aqueous phase. Represents gas density, Represents the seepage velocity in the gas phase. Represents the density of water. Represents the seepage velocity of the aqueous phase. It represents the area of ​​the current controlled volume; for a 3D controlled volume, this is the volume. Represents the boundary of the control volume. Represents the control body. It is the unit vector of the outward normal to the boundary of the control volume. The differential represents the control volume boundary.

[0186] Specifically, in step S 500 In the above, the expression for calculating the temperature increment of each grid vertex at the next moment is:

[0187] (10)

[0188] in: It is an increase in temperature. It is the weighted average of the phase heats. It is air compared to heat, Water is hotter than water. It is the specific heat of hydrates. It is the thermal conductivity tensor. It is the enthalpy change caused by the decomposition of hydrates. It is the effective stress. It is plastic strain. ( ) represents the divergence of temperature flux with respect to the control volume, which represents the energy change in the region due to heat conduction.

[0189] Therefore, the two-phase pressure, three-phase saturation and temperature of each grid vertex can be calculated according to the above algorithm. At the same time, the displacement of each node can be calculated according to the finite element method. In this way, all the basic information is stored on the grid vertex.

[0190] S 600 For the mechanical governing equations of the multiphysics coupled mathematical model of the hydrate-bearing sediments, the Galerkin finite element method is used to solve the purely elastic problem, and the incremental variable stiffness method is used to solve the elastoplastic problem.

[0191] Specifically, regarding step S of the present invention 600 The method of using the Galerkin finite element method to solve purely elastic problems specifically includes:

[0192] Update the elastic constants based on the hydrate saturation calculated in the previous step;

[0193] Based on the temperature and pore pressure calculated in the previous step, the displacement considering the thermal-fluid-structure interaction is calculated using a purely elastic finite element program.

[0194] The stress-strain and volumetric strain are obtained.

[0195] Specifically, the steps S of this invention 600 The method of using incremental variable stiffness to solve elastoplastic problems specifically includes:

[0196] Gradually apply the load and determine whether the stress in the mesh reaches the yield state;

[0197] Based on whether the yield state has been reached, the correction coefficient matrix is ​​modified, and the displacement solution of the target study area is finally obtained by iterative process.

[0198] The strain and stress of the target study area were calculated.

[0199] Please see Figure 7 As shown, according to Figure 7 By continuously calculating in the manner shown, the parameters of each physical field in the spatiotemporal domain can be solved.

[0200] S 700The reservoir porosity is updated based on the calculated volumetric strain, temperature, and pore pressure, and the absolute permeability is updated based on the porosity.

[0201] S 800 Determine whether the calculation results meet the stability conditions. If they do, continue the loop calculation. If they do not meet the stability conditions, reduce the time step and recalculate until the stability conditions are met. Output the calculation results of each physical field at the current time in the form of a field diagram.

[0202] Since the incremental explicit scheme is used in this embodiment of the invention, the time step cannot be too large. For simplicity, this embodiment of the invention directly adopts an incremental saturation control method. Specifically, an upper limit value is set for the increment of saturation and pressure. Once the calculation result exceeds this value, the time step is reduced to ensure the stability and reliability of the simulator.

[0203] Therefore, this embodiment of the invention conducts simulations based on in-situ experiments and existing benchmark problems for hydrate extraction, eliminating various program issues, ensuring good comparison with similar simulators, accurately solving differential-algebraic systems, and showing good comparability with indoor experimental results. When solving the seepage heat transfer system using an incremental approach, this method does not require solving linear equations, avoids the generation of large coefficient matrices, and saves memory. In addition, it has good adaptability when considering the elastoplastic characteristics of hydrate deposits.

[0204] In addition, the above simulator program design adopts a program framework based on a general structure, and ensures that it can adapt to regions with arbitrary geometric shapes and common two-dimensional and three-dimensional mesh types.

[0205] To enable the development of the aforementioned algorithm simulator, common scientific computing languages ​​such as Fortran, MATLAB, C, and C++ are used in this embodiment of the invention.

[0206] Using common scientific computing languages ​​such as Fortran, MATLAB, C, and C++, we will develop a simulator based on the aforementioned algorithms. Once the simulator is developed, we will encapsulate it and develop a user interface for users. We will continuously update and improve the numerical simulator based on user feedback.

[0207] The simulator uses existing preprocessing software such as Hypermesh, Gambit, and GMSH for mesh generation, and only needs to convert the mesh file into the data required for computation. Existing post-processing software such as Tecplot and Paraview is used for post-processing, and the simulator only needs to output the corresponding files.

[0208] Another embodiment of the present invention also provides a development apparatus for a numerical simulator, the development apparatus comprising:

[0209] Mesh partitioning cells are used to partition the target study area into meshes and output mesh information files;

[0210] The conversion unit is used to convert the grid information file into a matrix of coordinate information of each node in the target study area, a matrix of node numbers contained in each grid, and a matrix of grid numbers and node numbers in the boundary area.

[0211] The initialization unit is used to initialize temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation, and porosity.

[0212] The flux calculation unit is used to calculate the flux at each boundary of the control volume and accumulate them to obtain the flux difference for each control volume.

[0213] The incremental calculation unit is used to calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation, and hydrate saturation at each grid vertex in the next time step, based on the field parameters at the current time step.

[0214] The solution unit is used to solve the purely elastic problem using the Galerkin finite element method and the elastoplastic problem using the incremental variable stiffness method for the mechanical governing equations of the multiphysics coupled mathematical model of the hydrate-bearing sediments.

[0215] The update unit is used to update the reservoir porosity based on the calculated volumetric strain, temperature and pore pressure, and update the absolute permeability based on the porosity.

[0216] The decision unit is used to determine whether the calculation result meets the stability condition.

[0217] Another embodiment of the present invention provides a numerical simulator development system, including a lower-level simulator and a host computer that interacts with the lower-level simulator. The lower-level simulator is used to implement the numerical simulator development method described above, and the host computer is used to receive and parse the updated parsed data generated by the lower-level simulator.

[0218] Another embodiment of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the development method of the numerical simulator as described above.

Claims

1. A method for developing a numerical simulator for multiphysics coupling calculations of hydrate sediments, characterized in that, include: S 100 The target study area is meshed using mesh generation software, and the mesh information file is output. S 200 Import the grid information file and convert it into a matrix of coordinate information for each node in the target study area, a matrix of node numbers in each grid, and a matrix of grid numbers and node numbers in the boundary area. S 300 Based on the initial and boundary conditions of the temperature field, seepage field, mechanical field and phase transition field of the hydrate, the temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation and porosity are initialized. S 400 Construct a control volume centered on a grid vertex, calculate the flux of each boundary of the control volume, and accumulate them to obtain the flux difference of each control volume; The calculation of fluxes at each boundary of the control volume specifically includes water flux and gas flux, wherein: The flux of water about a certain cross section （1） Gas flux about a certain cross section （2） in: It is the permeability of the aqueous phase. It is the permeability of the gas phase. and These are the pore water pressure and pore gas pressure, respectively, assembled from the finite element basis functions and nodal function values. , , These are finite element basis functions. It is the number of grid vertices. It is the unit normal vector of that cross section. Represents the differential of the current cross section; S 500 Substitute the flux difference of each control volume into the incremental formula, and calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation and hydrate saturation of each grid vertex at the next moment based on the field parameters at the current moment. The specific expression for the pore gas pressure increment is as follows: （3） The specific expression for the increment of pore water pressure is as follows: （4） The specific expression for the gas phase saturation increment is as follows: （5） The specific expression for the increment of water phase saturation is as follows: （6） The specific expression for the hydrate saturation increment is as follows: （7） In the formula: , and , in: ( ) represents the divergence of water flux. ( ) represents the divergence of gas flux. Represents gas phase saturation. Represents the water phase saturation. Represents hydrate saturation. Represents effective water saturation. Represents the stiffness of the gas phase. Represents the stiffness of the water phase. Represents the stiffness of hydrates. It is the gas phase pore pressure. It is the pore pressure of the aqueous phase. It is porosity. It is volumetric strain. It is the time step. It is the molar mass of the gas phase. It is the molar mass of the aqueous phase. It is the molar mass of the hydrate. It is water and numbers. It is the gas phase density. It is the density of the aqueous phase. It is the density of the hydrate. It is the coefficient of thermal expansion of the gas phase. It is the coefficient of thermal expansion of water phase. It is the coefficient of thermal expansion of natural gas hydrate. It is the coefficient of thermal expansion of sediment particles. It is the number of moles of hydrate decomposition. It's temperature. It is capillary pressure. It is pore air pressure. It is pore water pressure. It is the gas saturation. It's water saturation. It is the increase in hydrate saturation; The divergence of the water flux ( The result is obtained from equation (8): （8） The divergence of the gas flux ( The result is obtained from equation (9): （9） In the formula: The number representing the component that controls the volume in each surrounding grid. This represents the portion of the control volume to which each grid cell belongs. Represents the density of water. Represents gas density, Represents the seepage velocity of the aqueous phase. Represents the seepage velocity in the gas phase. It is the area of ​​the current control body. Represents the control body. Represents the boundary of the control volume. It is the unit vector of the outward normal to the boundary of the control volume; The expression for calculating the temperature increment of each grid vertex at the next moment is as follows: （10） in: It is an increase in temperature. It is the weighted average of the phase heats. It is air compared to heat, Water is hotter than water. It is the specific heat of hydrates. It is the thermal conductivity tensor. It is the enthalpy change caused by the decomposition of hydrates. It is the effective stress. It is plastic strain. ( ) represents the divergence of heat flux with respect to the control volume; S 600 For the mechanical governing equations of the multi-physics coupled mathematical model of hydrate-bearing sediments, the Galerkin finite element method is used to solve the purely elastic problem, and the incremental variable stiffness method is used to solve the elastoplastic problem. S 700 The reservoir porosity is updated based on the calculated volumetric strain, temperature, and pore pressure, and the absolute permeability is updated based on the porosity. S 800 Determine whether the calculation results meet the stability conditions. If they do, continue the loop calculation. If they do not meet the conditions, reduce the time step and recalculate until the stability conditions are met. Output the calculation results of each physical field at the current time in the form of a field diagram.

2. The development method of the numerical simulator according to claim 1, characterized in that, In step S 600 The method of using the Galerkin finite element method to solve purely elastic problems specifically includes: According to step S 500 The calculated hydrate saturation is used to update the elastic constants; According to step S 500 The calculated temperature and pore pressure were used to calculate the displacement considering the thermal-fluid-structure interaction using a purely elastic finite element program. The stress-strain and volumetric strain were calculated.

3. The development method of the numerical simulator according to claim 1, characterized in that, In step S 600 The method of using incremental variable stiffness to solve elastoplastic problems specifically includes: Gradually apply the load and determine whether the stress in the mesh reaches the yield state; Based on whether the yield state has been reached, the correction coefficient matrix is ​​modified, and the displacement solution of the target study area is finally obtained by iterative process. The strain and stress of the target study area were calculated.

4. A development apparatus for a numerical simulator, characterized in that, include: Mesh partitioning cells are used to partition the target study area into meshes and output mesh information files; The conversion unit is used to convert the grid information file into a matrix of coordinate information of each node in the target study area, a matrix of node numbers contained in each grid, and a matrix of grid numbers and node numbers in the boundary area. The initialization unit is used to initialize temperature, pore pressure, displacement, stress, water saturation, gas saturation, hydrate saturation, and porosity. The flux calculation unit is used to calculate the flux at each boundary of the control volume and accumulate them to obtain the flux difference of each control volume. The calculation of the flux at each boundary of the control volume specifically includes water flux and gas flux. Among them, the flux of water about a certain cross section （1） Gas flux about a certain cross section （2） in: It is the permeability of the aqueous phase. It is the permeability of the gas phase. and These are the pore water pressure and pore gas pressure, respectively, assembled from the finite element basis functions and nodal function values. , , These are finite element basis functions. It is the number of grid vertices. It is the unit normal vector of that cross section. Represents the differential of the current cross section; The incremental calculation unit is used to calculate the increments of temperature, pore gas pressure, pore water pressure, gas phase saturation, water phase saturation, and hydrate saturation at each grid vertex in the next time step, based on the field parameters at the current time step. The specific expression for the pore gas pressure increment is as follows: （3） The specific expression for the increment of pore water pressure is as follows: （4） The specific expression for the gas phase saturation increment is as follows: （5） The specific expression for the increment of water phase saturation is as follows: （6） The specific expression for the hydrate saturation increment is as follows: （7） In the formula: , and , in: ( ) represents the divergence of water flux. ( ) represents the divergence of gas flux. Represents gas phase saturation. Represents the water phase saturation. Represents hydrate saturation. Represents effective water saturation. Represents the stiffness of the gas phase. Represents the stiffness of the water phase. Represents the stiffness of hydrates. It is the gas phase pore pressure. It is the pore pressure of the aqueous phase. It is porosity. It is volumetric strain. It is the time step. It is the molar mass of the gas phase. It is the molar mass of the aqueous phase. It is the molar mass of the hydrate. It is water and numbers. It is the gas phase density. It is the density of the aqueous phase. It is the density of the hydrate. It is the coefficient of thermal expansion of the gas phase. It is the coefficient of thermal expansion of water phase. It is the coefficient of thermal expansion of natural gas hydrate. It is the coefficient of thermal expansion of sediment particles. It is the number of moles of hydrate decomposition. It's temperature. It is capillary pressure. It is pore air pressure. It is pore water pressure. It is the gas saturation. It's water saturation. It is the increase in hydrate saturation; The divergence of the water flux ( The result is obtained from equation (8): （8） The divergence of the gas flux ( The result is obtained from equation (9): （9） In the formula: The number representing the component that controls the volume in each surrounding grid. This represents the portion of the control volume to which each grid cell belongs. Represents the density of water. Represents gas density, Represents the seepage velocity of the aqueous phase. Represents the seepage velocity in the gas phase. It is the area of ​​the current control body. Represents the control body. Represents the boundary of the control volume. It is the unit vector of the outward normal to the boundary of the control volume; The expression for calculating the temperature increment of each grid vertex at the next moment is as follows: （10） in: It is an increase in temperature. It is the weighted average of the phase heats. It is air compared to heat, Water is hotter than water. It is the specific heat of hydrates. It is the thermal conductivity tensor. It is the enthalpy change caused by the decomposition of hydrates. It is the effective stress. It is plastic strain. ( ) represents the divergence of heat flux with respect to the control volume; The solution element is used to solve the mechanical governing equations of a multiphysics coupled mathematical model for hydrate-bearing sediments. The Galerkin finite element method is used to solve the purely elastic problem, and the incremental variable stiffness method is used to solve the elastoplastic problem. The update unit is used to update the reservoir porosity based on the calculated volumetric strain, temperature and pore pressure, and update the absolute permeability based on the porosity. The decision unit is used to determine whether the calculation result meets the stability condition.

5. A development system for a numerical simulator, characterized in that, The system includes a lower-level simulator and a higher-level computer that interacts with the lower-level simulator. The lower-level simulator is used to implement the development method of the numerical simulator as described in any one of claims 1 to 3, and the higher-level computer is used to receive and parse the updated parsed data generated by the lower-level simulator.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor of the computer, it implements the development method of the numerical simulator as described in any one of claims 1 to 3.

Images