A Method for Establishing a Numerical Model for Hydrate Extraction Based on Start-up Pressure Gradient
By establishing a numerical model for hydrate exploitation based on the starting pressure gradient, the problem of slowed seepage velocity in low-permeability hydrate reservoirs in the South China Sea was solved, production capacity was increased, pressure and temperature distribution were optimized, bottom water coning was suppressed, and the gas-water ratio was improved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-05-10
- Publication Date
- 2026-03-10
AI Technical Summary
Existing technologies fail to effectively consider the characteristics of the initiation pressure gradient in low-permeability hydrate reservoirs in the South China Sea, resulting in a slow change in seepage velocity within low-permeability reservoirs, affecting well productivity, and causing poor injection-production connectivity.
A numerical model for hydrate exploitation based on the initiation pressure gradient was established. The quantitative relationship between the initiation pressure gradient and reservoir parameters was obtained through experiments on hydrate reservoir samples from the South China Sea. The non-Darcy flow section was characterized by a pseudo-initiation pressure gradient. The mass conservation equation in the hydrate exploitation simulation control module was improved. A numerical model was built by combining formation physical property data to simulate the hydrate exploitation process.
It improved the productivity of silty mudstone hydrate reservoirs in the South China Sea, increased gas production, reduced water production, optimized pressure transmission rate and temperature distribution, suppressed bottom water coning, improved wellbore water saturation, and enhanced the gas-water ratio.
Smart Images

Figure CN115292870B_ABST
Abstract
Description
[0001] Field of study
[0002] This invention belongs to the field of natural gas hydrate development, and specifically relates to a method for establishing a numerical model for hydrate extraction based on a starting pressure gradient. Background Technology
[0003] Natural gas hydrate (commonly known as combustible ice) is a cage-like crystalline compound formed by gas and water under certain temperature and pressure conditions. It is widely found in shallow seabed sediments and polar permafrost.
[0004] The hydrate reservoirs discovered in the Shenhu area of the South China Sea are all typical silty mudstone reservoirs. The Guangzhou Marine Geological Survey found at three geophysical stations (SH2, SH3, and SH7) in the Shenhu area that the hydrate saturation exceeded 20%, reaching as high as 47.3%. However, core analysis of the strata showed that the content of large-diameter sand was less than 10%, while the content of silty mudstone was as high as 80-90%, and the clay content of the hydrate producing layer was high.
[0005] The results of seepage experiments on fine-grained sediments with high clay content show that pore water flow in these sediments deviates significantly from the flow law described by the Navier-Stokes equations and does not conform to Darcy's law. The flow rate exhibits a non-linear relationship with the pressure gradient, and a starting pressure gradient exists; pore water only begins to flow when the pressure gradient exceeds the minimum starting pressure gradient. Therefore, the variation of the reservoir starting pressure gradient with formation seepage conditions and fluid properties during the development of silty mudstone hydrate reservoirs in the South China Sea warrants further investigation and clarification, and its application in practical field operations. In the exploitation of low-permeability oil and gas reservoirs, the starting pressure gradient has a significant impact on production capacity. Influenced by the starting pressure gradient, the seepage velocity change within low-permeability reservoirs exhibits a significant delay effect, resulting in a step-like change in seepage velocity. This significantly affects the pressure transmission range within the reservoir, leading to a rapid increase in water production rate in the later stages of well development. If the well spacing is too large or the injection-production pressure difference is too small, the starting pressure gradient prevents the formation of an effective injection-production connection.
[0006] With the steady progress of hydrate development in the South China Sea, numerous numerical simulation studies have been conducted on hydrate reservoirs. Most current research utilizes the TOUGH+HYDRATE hydrate simulator, which simulates the non-isothermal gas release, phase characteristics, fluid flow, and heat changes of natural gas hydrate reservoirs under various complex formation conditions by solving the mass and energy conservation equations. This simulator accurately describes all mechanisms of hydrate decomposition, including calculations of depressurization, heating, and chemical inhibitors. However, current studies have not considered the initiation pressure characteristic of low-permeability hydrate reservoirs in the South China Sea. Summary of the Invention
[0007] This invention provides a method for establishing a numerical model for hydrate mining based on the starting pressure gradient, which solves the problem that the starting pressure characteristics of low-permeability hydrate reservoirs in the South China Sea were not considered in the study of silty mudstone hydrate reservoirs in the South China Sea.
[0008] A method for establishing a numerical model for hydrate extraction based on the initiation pressure gradient includes:
[0009] Step 100: Obtain the quantitative relationship between the initiation pressure gradient of the South China Sea hydrate reservoir and the reservoir parameters through the initiation pressure gradient experiment of the South China Sea hydrate reservoir sample. Based on the quantitative relationship between the initiation pressure gradient of the South China Sea hydrate reservoir and the reservoir parameters, the pseudo-initiation pressure gradient is used to characterize the non-Darcy flow section, and Darcy's law is used to describe the seepage process of the non-Darcy flow section.
[0010] Step 200: Combining the seepage process of hydrates in low-permeability reservoirs in the South China Sea, improve the mass conservation equation in the hydrate mining simulation control module based on the starting pressure gradient, and establish a starting pressure gradient module to control the magnitude of the starting pressure gradient based on the experimental results of the starting pressure gradient of the hydrate reservoir samples. The starting pressure gradient in the hydrate mining simulation control module is controlled by the starting pressure gradient module.
[0011] Step 300: The starting pressure gradient module is coupled with the starting pressure gradient test results of the South China Sea hydrate reservoir samples, thereby improving the mass conservation equation in the hydrate mining simulation control module. The hydrate mining numerical model is built using the stratigraphic property data of the hydrate reservoir at the SH2 station in the South China Sea and the geological model parameters of the hydrate reservoir to simulate the hydrate mining process.
[0012] In a preferred embodiment, the quantitative relationship equation between the South China Sea hydrate reservoir initiation pressure gradient and reservoir parameters in step 100 is as follows:
[0013]
[0014] In the formula, λ is the starting pressure gradient, MPa·m -1 k is the penetration rate, 10 -3 μm 2 μ is the viscosity of the formation water, which is a constant.
[0015] In a preferred embodiment, the equation describing the seepage process in the non-Darcy seepage section using Darcy's law in step 100 is:
[0016]
[0017] In the formula, λ is the starting pressure gradient, MPa·m -1 k is the penetration rate, 10 -3 μm 2μ is the viscosity of the formation water, which is a constant. ρ is the seepage velocity, cm / s; p is the pressure, MPa.
[0018] In a preferred embodiment, the seepage process of South China Sea hydrates in the low-permeability reservoir in step 200 is as follows: the non-Darcy seepage section is characterized by the pseudo-starting pressure gradient, and the seepage process of the non-Darcy seepage section is described by Darcy's law.
[0019] In a preferred embodiment, step 200, which involves perfecting the mass conservation equation in the hydrate extraction simulation control module based on the starting pressure gradient, specifically involves:
[0020] The mobile phase in the liquid phase component conservation equation is modified and improved, with the water component w mass conservation equation modified as follows:
[0021]
[0022] In the formula: A represents the liquid phase; I represents the ice phase; G represents the gas phase; α represents A, I, or G; Porosity; S α ρ represents the saturation of each phase. α X represents the density of each phase; w α q represents the mass fraction of water components in each phase; A q represents the mass of the injected liquid phase. G Q represents the mass of the injected gas phase. w Water is produced by the decomposition of hydrates; k is the absolute permeability; k rA k rG These are the relative permeabilities of the aqueous phase and the gas phase, respectively; μ A μ G These are the liquid phase and gas phase viscosities, respectively. λ represents the pressure gradients of the liquid and gas phases, respectively; g is the gravitational acceleration; D is the depth difference; λ is the starting pressure gradient; and t is the mining time.
[0023] The mass conservation equation for liquid methane component m is modified as follows:
[0024]
[0025] In the formula: X m α Q represents the mass fraction of methane in each phase; m Gas is produced from the decomposition of hydrates;
[0026] The mass conservation equation for water-soluble component i, such as salt and inhibitors, is modified as follows:
[0027]
[0028] In the formula: X iα This represents the mass fraction of methane in each phase.
[0029] In a preferred embodiment, the equation for the initiation pressure gradient within the initiation pressure gradient module in step 200 is:
[0030]
[0031] In the formula, a and b are characteristic parameters. Characteristic parameter a controls the magnitude of the initiation pressure gradient, characteristic parameter b determines the relationship between the initiation pressure and reservoir parameters, and λ is the initiation pressure gradient, MPa·m. -1 k is the penetration rate, 10 -3 μm 2 μ is the viscosity of the formation water, which is a constant.
[0032] In a preferred embodiment, the construction of the hydrate extraction numerical model in step 300 specifically involves: based on the formation property data of the hydrate reservoir at the SH2 station in the South China Sea, simulating a columnar region, using non-uniform grid subdivision, and employing grid densification near the wellbore to establish a two-dimensional axisymmetric columnar grid numerical model for hydrate extraction.
[0033] In a preferred embodiment, the hydrate extraction numerical model includes an extraction well design module.
[0034] The well design module adopts a fixed production pressure differential method to form an injection-production connection.
[0035] In a preferred embodiment, the geological model parameters of the hydrate reservoir in step 300 include a relative permeability model and a van Genuchten model.
[0036] In a preferred embodiment, within the established numerical model for hydrate extraction, a simulation test of sensitivity to the starting pressure gradient is conducted. Specifically, different characteristic parameter 'a' values are set in the starting pressure gradient module to characterize the magnitude of different starting pressure gradient values. Then, the depressurization extraction process is simulated with and without considering the starting pressure gradient. The relationship between the starting pressure gradient value and / or water production and / or gas production and / or the hydrate complete decomposition zone and / or the reservoir three-phase zone and / or the bottom water coning and / or the bottom hole pressure difference is analyzed.
[0037] Compared with the prior art, the present invention has the following advantages:
[0038] (1) There is a starting pressure gradient in the seepage process of the muddy silt hydrate reservoir in the South China Sea. The present invention obtained that the ratio of the starting pressure gradient to the sample permeability-fluid viscosity shows a good power function relationship.
[0039] (2) In this invention, a corresponding numerical model for hydrate mining based on the proposed start-up pressure gradient module is established. This model can be applied to simulation software and the start-up pressure gradient function can be implemented in the simulation software. After using the start-up pressure gradient function, the running time of the simulation software can be increased.
[0040] (3) This invention simulates and calculates the depressurization extraction process with and without considering the starting pressure gradient. When considering the experimental starting pressure gradient, for the SH2 hydrate site in the South China Sea, there is an unexpected increase in production, and the gas production continues to increase over 5 years, while the water production gradually decreases, and the gas-water ratio is good. The pressure transmission speed slows down and is radial, with the pressure wave front almost perpendicular to the horizontal line. Within the pressure wave range, the pressure drop is even greater. The position of the front edge of the hydrate complete decomposition zone does not move forward much, but the range of the hydrate decomposition zone (gas, water, and hydrate three-phase zone) increases significantly, which also leads to an increase in the temperature drop within the reservoir, an expansion of the low-temperature region (<12℃), and even affects the temperature distribution of the underlying layer. Compared with not considering the starting pressure gradient, the water saturation around the production well section is higher, while when considering the starting pressure gradient, the well section is surrounded by a low water saturation zone. The existence of the starting pressure gradient also limits the conical advance of bottom water.
[0041] (4) Through the simulation mining of the present invention, it was found that the existence of the starting pressure gradient can suppress the bottom water coning situation during the vertical well depressurization mining process. The larger the starting pressure gradient, the smaller the volume of the intruded water and the lower the intrusion rate. Attached Figure Description
[0042] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.
[0043] Figure 1 This is a schematic diagram illustrating the process of establishing a numerical model for hydrate extraction based on the starting pressure gradient in this invention.
[0044] Figure 2 This is a schematic diagram of the pressure gradient initiation experiment process in this invention;
[0045] Figure 3 This is a schematic diagram of the physical model for hydrate extraction in this invention;
[0046] Figure 4 This is a schematic diagram illustrating the impact of the start-up pressure gradient on production dynamics in this invention;
[0047] Figure 5This is a schematic diagram illustrating the effect of the starting pressure gradient on the instantaneous gas-water ratio in this invention;
[0048] Figure 6 This is a schematic diagram illustrating the effect of the initiation pressure gradient on reservoir pressure distribution in this invention;
[0049] Figure 7 This is a schematic diagram illustrating the effect of the initiation pressure gradient on the reservoir temperature distribution in this invention;
[0050] Figure 8 This is a schematic diagram illustrating the effect of the starting pressure gradient on the hydrate saturation distribution in this invention;
[0051] Figure 9 This is a schematic diagram illustrating the effect of the starting pressure gradient on the water saturation distribution in this invention;
[0052] Figure 10 This is a schematic diagram illustrating the effect of the starting pressure gradient on the gas saturation distribution in this invention.
[0053] Figure 11 This is a schematic diagram comparing the gas production rates at different start-up pressures in this invention;
[0054] Figure 12 The diagram shows the pressure distribution of the 1800-meter stratum under different starting pressure gradient conditions in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0055] Figure 13 The diagram shows the temperature distribution of the 1800-meter stratum under different starting pressure gradients in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0056] Figure 14 The diagram shows the formation pressure funnel within 100m of the wellbore at a depth of 60m after 1800 days of mining under different starting pressure gradient conditions in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0057] Figure 15 The diagram shows the formation temperature funnel within 100m of the wellbore at a depth of 60m after 1800 days of mining at a constant pressure difference of 10MPa under different starting pressure gradient conditions in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0058] Figure 16 The diagram shows the distribution of hydrate saturation after 1800 days of mining under different starting pressure gradients in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0059] Figure 17 The diagram shows the water saturation distribution after 1800 days of mining under different starting pressure gradients in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4.
[0060] Figure 18 The diagram shows the gas saturation distribution after 1800 days of mining under different starting pressure gradients in this embodiment of the invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4;
[0061] Figure 19 The diagram shows a comparison of the bottom water conic evolution under different starting pressure gradients in the embodiments of the present invention. a) without considering the starting pressure gradient, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4;
[0062] Figure 20 The following are 3D diagrams comparing the evolution of the bottom water conic under different starting pressure gradients in the embodiments of the present invention: a) without considering the starting pressure gradient, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4;
[0063] Figure 21 The diagram shows the advancing front of the hydrate decomposition zone after 1800 days of mining under different starting pressure gradient conditions in the embodiments of the present invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4;
[0064] Figure 22 The diagram shows the advancing front of the hydrate complete decomposition zone after 1800 days of mining under different starting pressure gradient conditions in the embodiments of the present invention. a) Starting pressure gradient is not considered, b) a = 0.07, c) a = 0.14, d) a = 0.28, e) a = 0.7, f) a = 1.4. Detailed Implementation
[0065] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0066] To investigate the characteristics of the initiation pressure gradient in silty mudstone hydrate reservoirs in the South China Sea and its impact on the evolution of reservoir pressure, temperature, gas and water saturation, and hydrate saturation distribution during depressurization exploitation, this invention establishes a mathematical model of seepage based on the initiation pressure gradient, focusing on the actual seepage process of hydrates in low-permeability reservoirs. A method for implementing the initiation pressure gradient function in numerical simulation software is developed, and the initiation pressure gradient test results of target hydrate reservoir samples in the Shenhu area of the South China Sea are coupled. Through production simulation of the SH2 hydrate site in the Shenhu area of the South China Sea, the impact of the initiation pressure gradient on reservoir parameters and production dynamics during hydrate reservoir exploitation is explored, further understanding the characteristics of hydrate exploitation in the South China Sea and promoting the safe and efficient development of natural gas hydrates.
[0067] This invention provides a method for establishing a numerical model for hydrate extraction based on the initiation pressure gradient, such as... Figure 1 As shown, it includes the following steps:
[0068] Step 100: Obtain the quantitative relationship between the initiation pressure gradient of the South China Sea hydrate reservoir and the reservoir parameters through the initiation pressure gradient experiment of the South China Sea hydrate reservoir sample. Based on the quantitative relationship between the initiation pressure gradient of the South China Sea hydrate reservoir and the reservoir parameters, the pseudo-initiation pressure gradient is used to characterize the non-Darcy flow section, and Darcy's law is used to describe the seepage process of the non-Darcy flow section.
[0069] Due to the existence of low-velocity non-Darcy flow in media with high clay content, the initiation pressure gradient can concisely describe this phenomenon. Regarding the initiation pressure gradient experiment on hydrate reservoir samples, specifically: in-situ reservoir cores were used in this example. The selected experimental samples were actual reservoir core samples obtained by drilling from the Guangzhou Marine Geological Survey after hydrate decomposition. Brine was used as the simulated formation injection fluid, with a concentration of 3.5%. The experiment was conducted at room temperature and atmospheric pressure, and the experimental procedure is as follows: Figure 2 As shown, the numerical value of the initiation pressure gradient of the South China Sea hydrochemical reservoir was clarified through experimental data fitting, and a quantitative relationship between the initiation pressure gradient and reservoir parameters was established.
[0070] Experimental results show that the starting pressure gradient and the ratio of permeability to viscosity exhibit a good power function relationship. Using power function regression analysis, the equation relating the proposed starting pressure gradient and the permeability to viscosity ratio in the South China Sea hydrate reservoir is established as follows:
[0071]
[0072] In the formula, λ is the starting pressure gradient, MPa·m -1 k is the penetration rate, 10 -3 μm 2 The viscosity μ of formation water is generally considered to be a constant, and the initiation pressure gradient can be simplified to a simple power function of permeability.
[0073] In constructing the seepage equation, this invention uses the pseudo-starting pressure gradient to characterize the non-Darcy seepage section, that is, it uses Bingham fluid simulation to perform simulation calculations. By using the effective potential gradient method, as shown in equation (2), Darcy's law is used to describe its seepage process.
[0074]
[0075] In the formula ρ is the seepage velocity, cm / s; p is the pressure, MPa.
[0076] Step 200: Combining the seepage process of hydrates in low-permeability reservoirs in the South China Sea, improve the mass conservation equation in the hydrate mining simulation control module based on the starting pressure gradient, and establish a starting pressure gradient module to control the magnitude of the starting pressure gradient based on the experimental results of the starting pressure gradient of the hydrate reservoir samples in the South China Sea. The starting pressure gradient in the hydrate mining simulation control module is controlled by the starting pressure gradient module.
[0077] To establish a more reliable numerical model for hydrate mining in the South China Sea, considering the actual seepage process of hydrates in low-permeability reservoirs and the pressure gradient that forms an effective displacement, a corresponding mathematical model for hydrate mining seepage was established based on the hydrate mining simulation control module. This model was then applied to the simulation software to form a method for implementing the pressure gradient activation function.
[0078] The numerical model for hydrate extraction includes a hydrate extraction simulation control module, a start-up pressure gradient module, and an extraction well design module. In the hydrate extraction simulation control module, the simulation control equations for hydrate extraction have been modified and improved based on the actual seepage process of hydrates in low-permeability reservoirs in the South China Sea.
[0079] The governing equations for hydrate extraction simulation, including the mass and energy conservation of each component, can be expressed as follows:
[0080]
[0081] In the formula, the left side is the cumulative term, k represents different liquid phase components (water, methane, inhibitors, etc.), F is the mobile phase of matter and energy, q is the source-sink term, and n is the surface element dτ. n The normal vector, pointing towards V nt is the mining time, d; Vn is the volume, L 3 M is density, kg / m³ 3 ;τ n The area is the surface area of the region.
[0082] Since the gas phase permeation does not consider the influence of the starting pressure gradient, this invention modifies and improves the mobile phase in the conservation equation for the liquid phase components (water w, methane m, salt or inhibitor i).
[0083] Among them, the mass conservation equation for water component w is modified as shown in equation (4), and the starting pressure gradient λ is added to the liquid mobile phase part on the right.
[0084]
[0085] In the formula: A represents the liquid phase; I represents the ice phase; G represents the gas phase; α represents A, I, or G; Porosity; S α ρ represents the saturation of each phase. α X represents the density of each phase; w α q represents the mass fraction of water components in each phase; A q represents the mass of the injected liquid phase. G Q represents the mass of the injected gas phase. w Water is produced by the decomposition of hydrates; k is the absolute permeability; k rA k rG These are the relative permeabilities of the aqueous phase and the gas phase, respectively; μ A μ G These are the liquid phase and gas phase viscosities, respectively. λ represents the pressure gradients of the liquid and gas phases, respectively; g represents the gravitational acceleration; D represents the depth difference; λ represents the starting pressure gradient; and t represents the mining time.
[0086] The mass conservation equation for liquid methane m component is modified as shown in equation (5), with the starting pressure gradient λ added to the liquid mobile phase part on the right.
[0087]
[0088] In the formula: X m α Q represents the mass fraction of methane in each phase; m It produces gas from the decomposition of hydrates.
[0089] The mass conservation equation for water-soluble components i, such as salts and inhibitors, is modified as shown in equation (6):
[0090]
[0091] In the formula: X i α This represents the mass fraction of methane in each phase.
[0092] In summary, the mass conservation equation considering the initiation of the pressure gradient can be obtained. When considering the initiation of the pressure gradient, the thermal convection portion of the energy conservation equation is modified as follows:
[0093]
[0094] In the formula h A h G ν is the specific enthalpy of the liquid phase and the gas phase; T is the temperature.
[0095] In the simulation equations governing hydrate extraction, λ is determined by the initiation pressure gradient module. Based on the results of seepage experiments in the South China Sea reservoir, the initiation pressure gradient increases with decreasing permeability. At relatively high permeability, the formation water initiation pressure gradient increases slowly with decreasing permeability; at relatively low permeability, it increases sharply with decreasing permeability. The initiation pressure gradient exhibits a good power function relationship with the permeability-viscosity ratio. Therefore, the initiation pressure gradient module uses a power function to characterize the relationship between the initiation pressure gradient and the permeability-viscosity ratio, as follows:
[0096]
[0097] In the formula, a and b are characteristic parameters, with b typically being negative. Parameter a controls the magnitude of the startup pressure, while parameter b determines the relationship between the startup pressure and reservoir parameters. Entering the values of a and b in this module allows you to add a specific startup pressure gradient curve.
[0098] Step 300: The starting pressure gradient module is coupled with the starting pressure gradient test results of the South China Sea hydrate reservoir samples, thereby improving the mass conservation equation in the hydrate mining simulation control module. The hydrate mining numerical model is built using the stratigraphic property data of the hydrate reservoir at the SH2 station in the South China Sea and the geological model parameters of the hydrate reservoir to simulate the hydrate mining process.
[0099] In this embodiment, a model is established based on the stratigraphic property data of the SH2 water level station in the South China Sea to simulate the mining process of hydrates. The experimental values are used in the simulation of the production capacity of specific hydrate reservoirs in the South China Sea, namely a = 0.14 and b = -0.28, which are more in line with the actual situation of the mine.
[0100] The Shenhu sea area is located in the northern continental slope region of the South China Sea, within the Baiyun Depression of the Pearl River Estuary Basin. The Baiyun Depression has a Cenozoic sedimentary layer thickness exceeding 11 km, providing excellent geological conditions for hydrocarbon generation. This invention focuses on the SH2 hydrate reservoir in the Shenhu sea area. According to drilling and logging data from the Guangzhou Marine Geological Survey, the average permeability of the silty mudstone hydrate reservoir at the SH2 well site is 10 mD, the thickness of the hydrate-bearing sedimentary layer is 10-43 m, buried 185-229 m below the seabed, with a mudline depth of 1235 m and a mudline temperature of 3.9℃. The hydrate saturation in the developed hydrate ore bodies reaches a maximum of 48%, with an average saturation of 16.5%. The formation porosity is 40%, and methane accounts for 96.1-99.82% of the gas produced after hydrate decomposition. In-situ observations show a geothermal gradient of 43-67.7℃ / km.
[0101] Based on the above stratigraphic property data, a two-dimensional axisymmetric cylindrical grid numerical model for mining hydrate ore bodies was established. Figure 3 This invention describes a schematic diagram of a simulated natural gas hydrate reservoir model in the Shenhu area of the South China Sea. The simulation region is columnar, with a total of 120 × 110 = 13220 grids. Specifically, the longitudinal (z-direction) grid has 110 grids, a vertical range of 104m, and uses a non-uniform grid subdivision, with a maximum vertical grid size of 10m and a minimum vertical grid size of 0.5m. The transverse (r-direction) grid has 120 grids, a maximum radius of 1000m, and a 0.1m grid refinement near the wellbore.
[0102] The hydrate layer thickness is 44 m, and the thicknesses of the overlying and underlying layers are both 30 m. The bottom interface temperature of the hydrate stability zone is 15.01℃, and the pressure is 15.22 MPa. The top and bottom pressures are 14.5 MPa and 15.47 MPa, respectively, and the temperatures are 11.75℃ and 16.21℃, respectively. The hydrate layer contains two phases: an aqueous phase and a hydrate phase, with an initial hydrate saturation of 16.5%. The overlying and underlying layers are completely saturated with water, and the overall reservoir density is taken as 2600 kg / m³. 3 The initial permeability was set to 10 mD. Relevant formation parameters and physical properties are shown in Table 1.
[0103] Table 1 Physical parameters of the geological model of hydrate reservoirs
[0104]
[0105]
[0106] Where, k rG k rA S represents the relative permeability of gas and water, respectively. A * S G *S represents water saturation and gas saturation, respectively. G and S A S represents the saturation level of gas and water, respectively. irG and S irA P represents the saturation of bound gas and bound water, respectively. cap P0 represents capillary pressure, P0 represents original formation pressure, and v represents flow velocity.
[0107] In the production well design module, the production wellhead radius r w The depth is 0.1m, the perforation section is 34m, and it is located in the middle of the hydrate layer. The extraction adopts a fixed production pressure differential method to form an injection-production connection. Considering the factors of the starting pressure gradient, this invention sets a pressure reduction regime with a constant bottom hole flowing pressure of 4.5MPa. The pressure reduction extraction process was simulated and calculated with and without considering the starting pressure gradient, and the simulation time was 1800 days.
[0108] Figure 4 The paper presents a comparison of gas and water production over 1800 days of production, considering and not considering the start-up pressure gradient. Figure 5 As can be seen from this, when the start-up pressure gradient is not considered, the gas production rate increases from the peak production rate (>3000 sm) at the well opening. 3 / d) Rapidly dropped to 1266sm in 200 days 3 / d, and slowly decreases during the remaining production time, until the simulation ends, when the gas production drops to less than 800sm. 3 / d. When considering the initiation pressure gradient, gas production also drops sharply after well opening, but the decline is less severe than without considering the initiation pressure gradient, with the gas production rate falling to 1580 sm by day 200. 3 / d, but then the gas production began to rise, reaching 2732sm by the end of the simulation. 3 / d. Figure 4 c shows that, without considering the starting pressure gradient, the cumulative gas production over 5 years is only 178.4 × 10⁻⁶. 4 sm 3 Considering the starting pressure gradient, the total gas production reaches 382.1 × 10⁻⁶. 4 sm 3 More than twice that of the former. Similarly, the activation pressure gradient has a significant impact on water production, from... Figure 4 As can be seen from b, without considering the start-up pressure gradient, the change in the permeate rate is not significant, ranging from an initial 320m during the simulation. 3 / d slowly increased to 370m in the final stage 3 / d. However, considering the initiation pressure gradient, the flow of formation water is restricted, and the water production rate continuously decreases. By the end of the simulation, the water production rate had dropped to 68m. 3 / d, with a cumulative water production of only 18.2×10 4 m 3 This is significantly lower than the case without considering the startup pressure gradient (63.4 × 10⁻⁶). 4 m 3 ).
[0109] This result was unexpected, meaning that the start-up pressure gradient of silty clay reservoirs, as indicated by the experimental values, is beneficial for increasing production capacity. This is evident from the instantaneous gas-water ratio during production ( Figure 5 When considering the starting pressure gradient, the gas-water ratio gradually increases with production time, potentially exceeding 40 by the end of the simulation. Without this consideration, the gas-water ratio gradually decreases, approaching 2 by the end of the simulation. An excessively low gas-water ratio makes downhole gas-water separation difficult, which is extremely detrimental to actual production.
[0110] The influence of the initiation pressure gradient on production capacity can be identified by analyzing the evolution of reservoir parameters. To reveal the impact of the initiation pressure gradient on the advancement of the decomposition front during hydrate depressurization mining, pressure, temperature, and saturation distribution maps of the formation at the 50th, 500th, and 1800th days of mining were selected for both cases with and without considering the initiation pressure gradient, and a comparative analysis was conducted.
[0111] The initiation pressure gradient directly affects the pressure evolution and distribution of the South China Sea hydrate reservoir. Figure 6 a, b, and c represent the reservoir pressure distribution on days 50, 500, and 1800, respectively, without considering the initiation pressure gradient. Figure 6 Figures d, e, and f show the reservoir pressure distribution on days 50, 500, and 1800, respectively, considering the initial pressure gradient. As can be seen from the figures, the differences between the two scenarios are significant. Without considering the initial pressure gradient, the pressure diffusion front is approximately horizontal. Although the pressure wave reaches a large area, the pressure drop is low, with the low-pressure region (<10 MPa) front advancing less than 5 m. A notable characteristic of considering the initial pressure gradient is that the pressure drop propagates radially from the perforated section throughout the depressurization process. The pressure advance front is almost vertical. Although the wave reaches only 120 m, the pressure drop in the overlying and overburden layers is significant, and the pressure drop within the wave reaches a large extent. The low-pressure region (<10 MPa) front advances to 40 m, which is highly favorable for hydrate decomposition.
[0112] Hydrate decomposition is a typical endothermic reaction, and changes in the temperature field can often accurately predict the hydrate decomposition process. Figure 7 a, b, and c represent the reservoir temperature distribution on days 50, 500, and 1800, respectively, without considering the startup pressure gradient. Figure 7Figures d, e, and f show the reservoir temperature distribution on days 50, 500, and 1800, respectively, considering the initiation pressure gradient. Without considering the initiation pressure gradient, a low-temperature zone (<12℃) appears within 1-3m of the wellbore in the initial production stage. As production continues, the range of this low-temperature zone gradually decreases and shifts to the upper near-wellbore position of the hydrate reservoir. This indicates a gradual decrease in hydrate decomposition, which corresponds to a decrease in gas production. Figure 6 c shows that the intrusion of bottom water from the overlying strata also led to a rise in temperature at the bottom of the production well and an upward shift of the geothermal gradient line. Considering the starting pressure gradient, the low temperature range (<12℃) in the early stage of production extends to 0-10m and gradually expands in subsequent production. By the end of the simulation, the leading edge of the low temperature region has extended to 50m inside the reservoir, affecting the overlying strata and suppressing the temperature rise caused by the upward intrusion of bottom water.
[0113] Changes in hydrate saturation distribution are the ultimate effect of the initiation pressure gradient. Figure 8 a, b, and c represent the distributions of hydrate saturation on days 50, 500, and 1800, respectively, without considering the starting pressure gradient. Figure 8 Figures d, e, and f show the distribution of hydrate saturation on days 50, 500, and 1800, respectively, considering the start-up pressure gradient. The comparison between not considering and considering the start-up pressure gradient shows that the extent of the hydrate decomposition completion zone (the region containing only gas and water phases) is not significantly different at each production stage. The difference lies in the extent of the hydrate decomposition zone (the region where gas, water, and hydrate phases coexist). Considering the start-up pressure gradient, the leading edge extension rate of the hydrate decomposition zone is much higher than when the start-up pressure gradient is not considered. Figure 8 As can be seen from f, by day 1800, the leading edge of the decomposition zone had reached 80m, and the degree of hydrate decomposition was even greater within the decomposition zone. Figure 8 In section c, due to excessively low temperatures near the wellbore in the upper reservoir of the hydrate, there is a situation where the hydrate cannot decompose. Figure 8 This did not occur in f either.
[0114] The starting pressure gradient also has a significant impact on water and gas saturation. Figure 9 a, b, and c represent the distributions of water saturation on days 50, 500, and 1800, respectively, without considering the starting pressure gradient. Figure 8 d, e, and f represent the distribution of water saturation on days 50, 500, and 1800, respectively, considering the starting pressure gradient. We will focus on these distributions here. Figure 9 c and Figure 9 The comparison of f. As shown in Figure 9c, there is a high water saturation zone around the production well section, causing "water blockage" around the well. However, the starting pressure gradient restricts the migration of formation water or decomposed water far from the well, while the area around the well is a low water saturation zone. Figure 9In figure f, the water intrusion volume of both the upper and lower overburden layers is smaller than that in figure c. Without considering the starting pressure, the intrusion of bottom water into the wellbore also causes an increase in the later water production.
[0115] Figure 10 a, b, and c represent the gas saturation distributions on days 50, 500, and 1800, respectively, without considering the starting pressure gradient. Figure 10 Figures d, e, and f show the gas saturation distribution on days 50, 500, and 1800, respectively, considering the initiation pressure gradient. Under the initiation pressure gradient, the gas enrichment is high around the well section, and an enrichment zone exists at the junction of the reservoir and the overlying strata.
[0116] Within the established and comprehensive numerical model for hydrate extraction, a simulation experiment on the sensitivity of the initiation pressure gradient can be conducted to evaluate the impact of the magnitude of the initiation pressure gradient on the spatiotemporal evolution of multiphase and multifield processes. Specifically, different characteristic parameter 'a' values are set in the initiation pressure gradient module to characterize the magnitude of different initiation pressure gradient values. Then, the depressurization extraction process is simulated with and without considering the initiation pressure gradient. The relationship between the initiation pressure gradient value and / or water production and / or gas production and / or the hydrate complete decomposition zone and / or the reservoir three-phase zone and / or bottom water coning and / or bottom hole pressure difference is analyzed.
[0117] To evaluate the impact of the starting pressure gradient on the spatiotemporal evolution of multiphase and multifield processes, based on the hydrate depressurization mining model established in this invention and parameters a and b obtained from experimental measurement fitting equations, a set of different parameter a values were set to characterize different starting pressure gradient magnitudes. Table 2 shows the starting pressure gradient magnitudes corresponding to different a values.
[0118] Table 2 Model parameters and corresponding starting pressure gradients
[0119] Parameter a <![CDATA[Initial pressure gradient (MPa·m -1 )]]> 0 0 0.07 0.2542 0.14 0.5083 0.28 1.0166 0.7 2.5415 1.4 5.0831
[0120] Figure 11 The figure shows the variation of gas production rate over time under different starting pressure gradients at a constant bottom hole pressure. It can be seen that, without considering the starting pressure gradient, the gas production rate is low and continuously decreases. The optimal production effect is achieved when a = 0.28 and λ = 1.0166 MPa / m, which is twice the experimental value. On day 73, the gas production rate drops to 2129 sm. 3 After / d, it gradually increased, stabilizing at 2700sm after 800 days. 3 / d. When a continues to increase to 0.7, λ=2.5415MPa / m, similar to not considering the starting pressure gradient, the gas production rate shows a continuous decreasing trend, and at the end of the simulation, the gas production rates of the two are close (772sm). 3 / d、664sm 3 / d), but when a=0.7, the initial gas production rate decreases only slightly, and the production rate remains at 1000sm for the first 1000 days. 3 / d or higher. When a continues to increase to 1.4, the production effect is not as good as when the starting pressure gradient is not considered.
[0121] Figures 12 to 17 The distribution characteristics of formation pressure, temperature, hydrate saturation, water saturation, and gas saturation after 1800 days of mining at a constant bottomhole flowing pressure of 4.5 MPa are presented, considering different starting pressure gradients. It can be seen that when the starting pressure gradient is not considered, the pressure difference propagates rapidly into the deeper formation as mining progresses. Figure 12 As can be seen from a, after 1800 days of mining, the pressure difference can be effectively transmitted to a distance of 150m, and simultaneously from... Figure 13 As can be seen from a, due to the lower temperature at the top of the hydrate layer, the hydrate preferentially decomposes from the bottom during the depressurization mining process and continues to the distance.
[0122] When considering the initiation pressure gradient and its variation from small to large, the pressure differential transmission exhibits a unique pattern. First, we analyze the case where the initiation pressure gradient is small, i.e., a ≤ 0.14. This unique pattern manifests in two main aspects: Laterally, the pressure differential transmission differs significantly from that without considering the initiation pressure gradient; under large pressure differential production conditions, it exhibits radial transmission and propagates uniformly outwards. Vertically, considering the initiation pressure gradient, the hydrate layer effectively "blocks" the pressure differential, possessing a focusing effect, causing the pressure differential to propagate only near the wellbore periphery. (Comparison) Figure 12 a and Figure 12 b shows that when the initiation pressure gradient is considered and it is small, the pressure difference is effectively transmitted within a certain range near the wellbore, thus more hydrates participate in decomposition, and the reservoir temperature drops sharply. However, as the initiation pressure gradient gradually increases, the lateral transmission of the pressure difference is gradually blocked, making it difficult for hydrates to decompose, especially when a = 1.4 ( Figure 13 f) When the hydrate has been extracted for 1800 days, the effective decomposition range is only within 20m. Correspondingly, the reservoir temperature also decreases within the same range.
[0123] Figure 14 and Figure 15 Figure 15 presents formation pressure and temperature funnel diagrams within 100m of the wellbore at a depth of 60m after 1800 days of mining under different initiation pressure gradients. The pressure funnel diagram shows that, considering the initiation pressure gradient, the funnel radius decreases with increasing initiation pressure gradient, and the funnel shape also changes. Regarding the temperature funnel diagram shown in Figure 15, although it generally shows a trend of smaller funnel radius with larger initiation pressure gradients, a comparison reveals that... Figure 15For b-15d, the radius of the funnel at temperatures below 11℃ shows a trend of increasing with the increase of the starting pressure gradient, and then the radius of the funnel begins to decrease as the starting pressure gradient further increases.
[0124] Figure 16-18 The distribution characteristics of hydrate saturation, water saturation, and gas saturation after 1800 days of production are presented, considering different initiation pressure gradients. When the initiation pressure gradient is not considered, the decomposition front at the bottom of the hydrate layer extends far into the well, with an effective decomposition zone reaching 25m, and a completely decomposed zone within 15m near the wellbore. When the initiation pressure gradient is considered, the reservoir hydrate region extends significantly, resulting in a large-scale three-phase zone where gas, water, and hydrates coexist. Figure 16 In section b, the hydrate decomposition zone has expanded to 110m from the wellhead, but as the starting pressure gradient gradually increases, the width of the three-phase coexistence zone decreases accordingly. Figure 16 In f, there is almost no three-phase region. Figure 17 As can be seen, without considering the starting pressure gradient, there is a high-saturation water-bearing zone around the production well section. When the starting pressure gradient is small (Fig. 17b, c, d), the high-saturation water-bearing zone around the well disappears and is replaced by a gas-rich zone (Fig. 18b, c, d). However, when the pressure gradient continues to increase, the high-saturation water-bearing zone reappears, but the area near the wellhead is still a gas-rich zone.
[0125] Since the production well in the simulation did not penetrate the hydrate layer, the results showed bottom water coning. Figure 19 To compare the evolution of the bottom water cone under different initiation pressure gradients, Figure 20 This is a 3D rendering. Without considering the initial pressure gradient, the bottom water cone advance is significant; by day 1800, the intrusion height can reach 7m, and the water cone radius exceeds 25m. Figure 19 In case b, the height of the water cone did not decrease significantly, but the radius decreased to 16m. As the starting pressure increased, the height of the intruding water gradually decreased, and the radius of the water cone also decreased rapidly. Figure 19 The height of the e-type has been reduced to less than 2m, and the radius has shrunk to less than 5m. Figure 19 The water cone phenomenon in f has even disappeared.
[0126] Figure 22 The study presents the advancement of the hydrate decomposition zone front after 1800 days of constant bottomhole flowing pressure production under different initiation pressure gradients; the blank area represents the region where hydrate decomposition has begun. With an initiation pressure gradient, the smaller the gradient, the larger the hydrate decomposition zone. Compared to the unconsidered case, the initiation pressure gradient can effectively improve the mobilization of hydrates in the upper reservoir. Figure 22This paper presents the advancement of the hydrate complete decomposition zone front after 1800 days of production at constant bottomhole flowing pressure under different starting pressure gradients, i.e., the hydrate saturation in the blank area is 0. Although the pressure, temperature, and saturation distribution of each phase in the reservoir vary greatly under different starting pressure conditions, the location of the complete decomposition zone front is roughly the same, around 15m. The difference lies only in the hydrate decomposition at the interface between the reservoir and the upper and lower caprocks. When the starting pressure gradient promotes the pressure drop in the pressure sweep zone, the hydrates at the boundary are completely decomposed due to the heat supply from the upper and lower caprocks. Figure 18 Large areas of hydrate decomposition are constrained by energy supply issues, resulting in slow decomposition rates and incomplete decomposition. (Contact) Figure 21 and Figure 22 As the value of 'a' increases, the distance between the leading edge of the fully decomposed region and the leading edge of the decomposed region decreases. Figure 18 f and Figure 19 In terms of f, the two almost overlap.
[0127] In this invention, a sensitivity analysis of the start-up pressure gradient was conducted using a numerical model for hydrate extraction. The start-up pressure gradient is beneficial for increasing production capacity. However, when the start-up pressure gradient is large, water production is further restricted, leading to ineffective pressure differential transmission, failure of hydrates to decompose in distant wells, and ultimately a continuous decrease in gas production. During production, the start-up pressure gradient has little impact on the advancement of the front edge of the fully decomposed hydrate zone, only affecting the area at the boundary between the upper and lower caprocks. However, it has a strong influence on the boundary of the hydrate decomposition zone (three-phase zone). A small start-up pressure gradient is beneficial for the expansion of the three-phase zone, while a gradually increasing start-up pressure gradient restricts the development of the three-phase zone, affecting production. As the start-up pressure gradient increases, the front edge of the three-phase zone gradually approaches the front edge of the fully decomposed hydrate zone.
[0128] The above embodiments are merely exemplary embodiments of this application and are not intended to limit this application. The scope of protection of this application is defined by the claims. Those skilled in the art can make various modifications or equivalent substitutions to this application within its substance and scope of protection, and such modifications or equivalent substitutions should also be considered to fall within the scope of protection of this application.
Claims
1. A method for establishing a hydrate production numerical model based on a starting pressure gradient, characterized in that, The method comprises the following steps: Step 100, obtaining a quantitative relationship between a South China hydrate reservoir threshold pressure gradient and reservoir parameters through a threshold pressure gradient experiment of a South China hydrate reservoir sample, using a pseudo threshold pressure gradient to represent a non-Darcy flow section based on the quantitative relationship between the South China hydrate reservoir threshold pressure gradient and the reservoir parameters, using Darcy's law to describe a seepage process of the non-Darcy flow section, and using a Bingham fluid simulation method to perform simulation calculation; Step 200, combining a seepage process of South China hydrate in a low-permeability reservoir, perfecting a mass conservation equation in a hydrate production simulation control module based on a threshold pressure gradient, and establishing a threshold pressure gradient module to control the size of the threshold pressure gradient according to the threshold pressure gradient experiment results of the South China hydrate reservoir sample, wherein the threshold pressure gradient in the hydrate production simulation control module is controlled in size by the threshold pressure gradient module; Step 300, coupling the threshold pressure gradient test results of the South China hydrate reservoir sample in the threshold pressure gradient module, and then perfecting the mass conservation equation in the hydrate production simulation control module, building the hydrate production numerical model based on the formation physical property data and the hydrate reservoir geological model parameters of the South China SH2 station hydrate reservoir, and simulating a hydrate production process. The quantitative relationship equation between the South China hydrate reservoir threshold pressure gradient and the reservoir parameters in step 100 is: (1); where λ is the threshold pressure gradient, MPa-m −1 ; k is the permeability, 10 −3 m 2 ; μ is the viscosity of the formation water, a constant; In step 200, the mass conservation equation in the hydrate production simulation control module is perfected based on the threshold pressure gradient, and specifically: The flow phase in the liquid phase component conservation equation is modified and perfected, wherein the water component w mass conservation equation is modified as follows, (4); where: A is liquid phase; I is ice phase; G is gas phase; a is A or I or G; is porosity; S α is saturation of each phase; p α is density of each phase; X w α is mass fraction of water component in each phase; q A is mass of injected liquid phase; q G is mass of injected gas phase; Q w is water produced by hydrate decomposition; k is absolute permeability; k rA , k rG are relative permeability of water phase and gas phase respectively; m A , m G are viscosity of liquid phase and gas phase respectively; , are pressure gradient of liquid phase and gas phase respectively; g is acceleration of gravity; D is depth difference; l is threshold pressure gradient; t is production time; The liquid phase methane m component mass conservation equation is modified as follows, (5); wherein: X m α is the mass fraction of the methane component in each phase; Q m is the hydrate dissociation gas; The water-soluble component i component mass conservation equation is modified as follows, wherein the water-soluble components include salt and inhibitor: (6); wherein: X i α is the mass fraction of the methane component in each phase; The hydrate reservoir geological model parameters in step 300 include a relative permeability model and a van Genuchten model.
2. The method of claim 1, wherein, The equation for using Darcy's law to describe the seepage process of the non-Darcy flow section in step 100 is: (2); where λ is the threshold pressure gradient, MPa-m −1 ; k is the permeability, 10 −3 m 2 ; μ is the viscosity of the formation water, which is a constant, is the seepage velocity, cm / s; p is the pressure, MPa; grad is the gradient of the pressure .
3. The method of claim 1, wherein, In step 200, the seepage process of South China hydrate in a low-permeability reservoir is that a pseudo threshold pressure gradient is used to represent a non-Darcy flow section, and Darcy's law is used to describe the seepage process of the non-Darcy flow section.
4. The method of claim 1, wherein, The equation for the threshold pressure gradient in the threshold pressure gradient module in step 200 is: (8); In the formula, a and b are characteristic parameters, the characteristic parameter a controls the size of the threshold pressure gradient, the characteristic parameter b determines the relationship between the threshold pressure and the reservoir parameters, λ is the threshold pressure gradient, MPa·m −1 ; k is the permeability, 10 −3 μm 2 ; μ is the viscosity of the formation water, which is a constant.
5. The method of claim 1, wherein, In step 300, the hydrate production numerical model is built, specifically: based on the formation physical property data of the South China SH2 station hydrate reservoir, the simulation area is columnar, non-uniform grid partitioning is used, grid encryption partitioning is used near the well, and the hydrate production numerical model of the two-dimensional axisymmetric columnar grid is established.
6. The method of claim 1, wherein, The hydrate production numerical model comprises a production well design module, and the production well design module uses a fixed production pressure difference method to form an injection-production connection relationship.
7. A method of establishing a hydrate production numerical model based on the start-up pressure gradient according to any one of claims 1 to 6, characterized in that, In the established perfect hydrate exploitation numerical model, a start-up pressure gradient sensitivity simulation test is carried out. Specifically, different characteristic parameter a values are set in the start-up pressure gradient module to represent different start-up pressure gradient values, and then the depressurization exploitation processes are simulated without considering the start-up pressure gradient and considering the start-up pressure gradient respectively, and the relationship between the start-up pressure gradient value and / or water production and / or gas production and / or hydrate complete decomposition zone and / or reservoir three-phase zone and / or bottom water coning and / or bottom hole differential pressure is analyzed.
Citation Information
Patent Citations
Simulation method and simulation system considering oil reservoir mechanism for shale oil and gas components
CN107133411A
Oil-water two-phase non-darcy seepage numerical simulation method based on discrete fracture model
CN111104766A