Gas-liquid two-phase flow numerical simulation method considering coal bed gas desorption lag

Abstract

The invention discloses a gas-liquid two-phase flow numerical simulation method considering coal bed gas desorption lag. The method comprises the following steps: acquiring productivity calculation parameters by using logging and indoor experiment means; defining a desorption equation and a diffusion equation of the matrix system in combination with a desorption hysteresis effect; defining a gas-water two-phase seepage equation of the fracture system, and forming a dual-medium gas-water two-phase coupling flow equation set in combination with a desorption equation and a diffusion equation of the matrix system; the porosity and permeability changes of a matrix system and a fracture system in the production process are calculated respectively, and a physical property change equation set of the double media is deduced; and solving the coupling equation set by using a numerical solver of COMSOL and adopting an iteration method, and calculating to obtain gas production and water production. The method is easy and convenient to operate, the influence of the coal bed gas desorption lag effect on the gas well productivity can be evaluated more scientifically and reasonably, meanwhile, the gas well numerical simulation calculation cost and time are greatly reduced, and important inspiration is also provided for development of other unconventional gas reservoirs influenced by desorption lag.

Description

Technical Field

[0001] The present invention relates to the technical field of coalbed methane development, and in particular to a gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis. Background Art

[0002] my country boasts abundant coalbed methane (CBM) resources and broad development prospects. CBM is primarily stored in an adsorbed state within the coal matrix. When pressure decreases, it desorbs, diffuses into fracture systems, and migrates to production wells via seepage. Compared to conventional natural gas reservoirs, CBM reservoirs exhibit distinctly low porosity and permeability. Their occurrence and migration are controlled by coal seam heterogeneity, stress sensitivity, and the desorption hysteresis effect. Desorption hysteresis refers to the hysteresis between the desorption and adsorption processes during CBM pressure fluctuations. Specifically, when CBM pressure decreases, the desorption process fails to remain completely symmetrical with the adsorption process, and the desorption curve exhibits a lag relative to the adsorption curve, meaning that the desorption process experiences a longer delay than the adsorption process. This hysteresis indicates that the desorption process is affected by pressure, preventing CBM from immediately reaching a new equilibrium desorption state when pressure changes. Accurately assessing the impact of desorption hysteresis on gas-liquid two-phase flow and optimizing CBM development strategies is a key research priority.

[0003] Currently, engineers typically use numerical simulation methods for production capacity prediction and optimization analysis in coalbed methane (CBM) reservoir development research. Mainstream commercial numerical simulation software (such as CMG and ECLIPSE) can simulate CBM seepage characteristics by building complex geological models and creating fine grids. However, these software applications require detailed modeling of the dual-media characteristics and dynamic seepage mechanisms of coalbed methane. This computational complexity, high computer hardware requirements, and lengthy simulation processes make it difficult to meet the demand for efficient decision-making during CBM development. Furthermore, traditional numerical simulation methods are primarily based on the equilibrium desorption assumption and fail to fully consider the dynamic effects of desorption hysteresis on gas-liquid two-phase flow, leading to biased predictions and, in turn, hindering the optimization of CBM reservoir development plans. Therefore, it is necessary to develop a numerical simulation method for gas-liquid two-phase flow that accounts for the hysteresis effect of CBM desorption and the complex seepage mechanisms of gas-liquid two-phase flow. This approach can improve the accuracy of CBM reservoir production capacity predictions and provide theoretical and technical support for the scientific development of CBM. Summary of the Invention

[0004] To solve the above problems, the present invention proposes a numerical simulation method for gas-liquid two-phase flow considering coalbed methane desorption hysteresis. This method takes into account the coalbed methane desorption hysteresis effect and the complex seepage mechanism of gas-liquid two-phase flow, and can reasonably evaluate coalbed methane development.

[0005] In order to achieve the above technical objectives, the present invention adopts the following technical solutions:

[0006] A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis, characterized by comprising the following steps:

[0007] S1: Use well logging and indoor experiments to obtain production capacity prediction parameters;

[0008] S2: Establish the computational domain of the CBM dual-medium two-phase model in COMSOL, and define the desorption equation and diffusion equation of the matrix system in combination with the desorption hysteresis effect;

[0009] S3: Using the finite element method of COMSOL, the gas-water two-phase seepage equation of the fracture system is defined, and combined with the desorption equation and diffusion equation of the matrix system to form a dual-medium two-phase coupled flow equation system;

[0010] S4: Considering the desorption hysteresis effect and the gas-water two-phase seepage characteristics, the porosity and permeability changes of the matrix system and the fracture system are calculated respectively, and the physical property change equations of the dual medium are derived;

[0011] S5: Using the numerical solver of COMSOL, set the corresponding initial conditions, and use the iterative method to solve the coupled equations in steps S2-S4, and obtain the gas and water production through calculation.

[0012] Furthermore, steps S1 to S5 are subject to the following assumptions:

[0013] 1) The coal body is divided into two systems: matrix and fracture. Adsorbed gas mainly exists in the matrix system, fractures are the main seepage channels for fluids, and water only exists in the fracture system.

[0014] 2) The strain of the coal body is considered to be infinitesimal;

[0015] 3) Adsorbed gas desorbs from matrix pores and enters the fracture system, subsequently flowing into the wellbore;

[0016] 4) Ignore the effects of gravity and temperature;

[0017] 5) The gas is equivalent to an ideal gas, and its viscosity remains unchanged during the flow process. Water is an incompressible fluid;

[0018] 6) The adsorption of methane in the matrix is saturated adsorption.

[0019] Furthermore, the capacity forecast parameters in step S1 include:

[0020] 1) Coal seam parameters

[0021] Coal seam length, coal seam width, coal seam thickness, coal seam density, matrix initial porosity, matrix initial permeability, fracture initial porosity, fracture initial permeability, coal seam temperature, coal seam initial pressure, coal seam Young's modulus, initial water saturation, irreducible water saturation;

[0022] 2) Adsorption and desorption parameters

[0023] Langmuir pressure, Langmuir volume, desorption hysteresis coefficient;

[0024] 3) Fluid parameters

[0025] Methane viscosity, methane density, water viscosity, water density, bottom hole flowing pressure, standard temperature.

[0026] Furthermore, the desorption equation of the matrix system in step S2 is:

[0027] W de =W ad +a·exp(b·p m )

[0028] Where a is the hysteresis gas content coefficient, m 3 / kg; b is the hysteresis gas pressure coefficient, MPa -1 , coefficients a and b can be fitted in actual adsorption and desorption experiments. m is the matrix system pressure, MPa; W de is the gas content of coal rock during desorption, m 3 / kg,W ad is the adsorbed gas content in the coal seam, m 3 / kg, which can be defined by the following formula:

[0029]

[0030] Where V L is the Langmuir volume constant of the adsorption process, m 3 / kg; P L is the Langmuir pressure constant of the adsorption process, MPa; V L-de is the Langmuir volume constant of the desorption process, m 3 / kg; P L-de is the Langmuir pressure constant of the desorption process, MPa; W re is the residual adsorption capacity, m 3 / kg.

[0031] The desorption hysteresis pressure is a threshold pressure in the desorption state, which is defined as p q , MPa, can be calculated according to the following formula:

[0032]

[0033] Where ρ c is the density of coal rock, kg / m 3 ; VM is the molar volume of methane gas, m 3 / mol; R is the gas molar constant, J / (mol·K); T is the temperature, K

[0034] According to Fick's diffusion law, the mass conservation equation of gas in the matrix is:

[0035]

[0036] where Q m =Dσ c (c m -c f ), represents the mass transfer rate of methane between the matrix and the fracture system, kg / (m 3 ·s);m m is the mass of methane per unit volume of coal matrix, kg / m 3 σ c is the coal matrix shape factor, m -2 , defined as L is the crack spacing, m; D is the diffusion coefficient, m 2 / s.

[0037] c m and c f are the concentrations of methane in the matrix and in the fracture system, kg / m 3 According to the ideal gas law:

[0038]

[0039] Where M g is the molar mass of methane, kg / mol; p fg is the gas pressure in the fracture system, MPa. Considering the desorption hysteresis effect, the mass transfer rate of methane between the matrix and the fracture system is obtained as:

[0040]

[0041] The adsorbed gas in the matrix is represented by the Langmuir equation, while the free gas satisfies the ideal gas state equation. The total gas content in the matrix system can be calculated by the following formula:

[0042]

[0043] Where φ m is the matrix porosity, dimensionless.

[0044] Furthermore, the diffusion equation of the matrix system in step S2 is:

[0045]

[0046] Furthermore, the gas-water two-phase seepage equation of the fracture system in step S3 can be written as:

[0047]

[0048] Where m g =ρ g φ f S g , represents the mass of free gas in cracks per volume of coal, kg / m 3 ;m w =ρ w φ f S w , represents the mass of water in the cracks per volume of coal, kg / m 3 ρ g is the gas density, kg / m 3 ρ w is the density of water, kg / m 3 ;S g and S w are gas and water saturations, dimensionless; φ f is the porosity of the fracture, dimensionless; k f is the fracture permeability, mD; μ g and μ w are the dynamic viscosities of methane and water, Pa·s; b1 is the Klinkenberg coefficient, Pa; k rg and k rw are the relative permeabilities of the gas phase and the water phase, dimensionless; p fw is the partial pressure of water in the fracture system.

[0049] The relative permeability of gas and water in a fracture system can be given by the following formula:

[0050]

[0051] Where m is the relative permeability fitting coefficient, dimensionless, and in this study, m = 0.65; S e is the effective water saturation, dimensionless, and can be defined by the following formula:

[0052]

[0053] Where S wr is the residual water saturation, S gr is the residual gas saturation.

[0054] Furthermore, the dual-medium two-phase coupled flow equations in step S3 are:

[0055]

[0056] Furthermore, considering the desorption hysteresis in step S4, the porosity change of the matrix system can be expressed as follows:

[0057]

[0058] Where S = ε V +(p m / K s )-ε sq , ε V is the volume strain of coal, dimensionless; K s is the bulk modulus of the coal skeleton, GPa; ε sq represents the strain during gas desorption, ε sq =α sg W de , dimensionless; α sg is the matrix expansion coefficient, kg / m 3 , where the subscript '0' represents the initial value of the corresponding variable.

[0059] Furthermore, considering the desorption hysteresis in step S4, the porosity change of the fracture system can be expressed as follows:

[0060]

[0061] Where M is the axial constraint modulus, MPa, and M = E(1-v) / [(1+v)(1-2v)]. Considering the cubic law between the permeability and porosity of porous media, the permeability of the matrix and fracture system can be expressed as:

[0062]

[0063] Furthermore, the formula for calculating the gas and water production in step S5 is:

[0064]

[0065] where Q g and Q w are daily gas and water production, m 3 ; and are the average pressures of gas phase and water phase, MPa; p wf is the bottom hole pressure, MPa; r w is the wellbore radius, m; r e is the drainage radius, m; S is the skin coefficient.

[0066] This paper provides a numerical simulation method for gas-liquid two-phase flow that takes into account coalbed methane desorption hysteresis. Combining well logging and experimental data, the method defines the desorption and diffusion equations for the matrix system and the gas-water two-phase flow equations for the fracture system. This system forms a set of coupled gas-water flow equations for the dual-medium system. The coupled equations are iteratively solved using COMSOL software to calculate gas and water production. This method is simple to operate, and actual testing has shown good fit with field data. It enables a more scientific and rational assessment of the impact of coalbed methane desorption hysteresis on gas well productivity, while significantly reducing the computational cost and time of gas well numerical simulations. This method also provides important insights for the development of other unconventional gas reservoirs affected by desorption hysteresis.

[0067] Beneficial effects:

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

[0069] COMSOL software was used to rapidly simulate gas-liquid two-phase flow under varying desorption hysteresis characteristics, and the model's accuracy was verified using coalbed methane well production data. This method enables a more scientific and rational assessment of the impact of desorption hysteresis on coalbed methane well productivity, improving numerical simulation accuracy while significantly reducing computational cost and time. This method holds significant application value in the economical and efficient development of coalbed methane and also offers valuable guidance for the development of other unconventional gas reservoirs affected by desorption hysteresis. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 This is a comparison chart of the production volume predicted for a vertical well and production data. DETAILED DESCRIPTION

[0071] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0072] Example:

[0073] A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis, characterized by comprising the following steps:

[0074] S1: Use well logging and indoor experiments to obtain production capacity prediction parameters;

[0075] The capacity calculation parameters include:

[0076] 1) Coal seam parameters

[0077] Coal seam length, coal seam width, coal seam thickness, coal seam density, matrix initial porosity, matrix initial permeability, fracture initial porosity, fracture initial permeability, coal seam temperature, coal seam initial pressure, coal seam Young's modulus, initial water saturation, irreducible water saturation;

[0078] 2) Adsorption and desorption parameters

[0079] Langmuir pressure, Langmuir volume, desorption hysteresis coefficient;

[0080] 3) Fluid parameters

[0081] Methane viscosity, methane density, water viscosity, water density, bottom hole flowing pressure, standard temperature.

[0082] S2: Establish the computational domain of the coalbed methane dual-medium two-phase model in COMSOL, and define the desorption equation and diffusion equation of the matrix system in combination with the desorption hysteresis effect;

[0083] The desorption equation of the matrix system is:

[0084] W de =W ad +a·exp(b·p m )

[0085] Where a is the hysteresis gas content coefficient, m 3 / kg; b is the hysteresis gas pressure coefficient, MPa -1 , coefficients a and b can be fitted in actual adsorption and desorption experiments. m is the matrix system pressure, MPa; W de is the gas content of coal rock during desorption, m 3 / kg,W ad is the adsorbed gas content in the coal seam, m 3 / kg, which can be defined by the following formula:

[0086]

[0087] Where V L is the Langmuir volume constant of the adsorption process, m 3 / kg; P L is the Langmuir pressure constant of the adsorption process, MPa; V L-de is the Langmuir volume constant of the desorption process, m 3 / kg; P L-de is the Langmuir pressure constant of the desorption process, MPa; W re is the residual adsorption capacity, m 3 / kg.

[0088] The desorption hysteresis pressure is a threshold pressure in the desorption state, which is defined as p q , MPa, can be calculated according to the following formula:

[0089]

[0090] Where ρ c is the density of coal rock, kg / m 3 ; V M is the molar volume of methane gas, m 3 / mol; R is the gas molar constant, J / (mol·K); T is the temperature, K

[0091] According to Fick's diffusion law, the mass conservation equation of gas in the matrix is:

[0092]

[0093] where Q m =Dσ c (c m -c f ), represents the mass transfer rate of methane between the matrix and the fracture system, kg / (m 3 ·s);m m is the mass of methane per unit volume of coal matrix, kg / m 3 σ c is the coal matrix shape factor, m -2 , defined as L is the crack spacing, m; D is the diffusion coefficient, m 2 / s.

[0094] c m and c f are the concentrations of methane in the matrix and in the fracture system, kg / m 3 According to the ideal gas law:

[0095]

[0096] Where M g is the molar mass of methane, kg / mol; p fg is the gas pressure in the fracture system, MPa. Considering the desorption hysteresis effect, the mass transfer rate of methane between the matrix and the fracture system is obtained as:

[0097]

[0098] The adsorbed gas in the matrix is represented by the Langmuir equation, while the free gas satisfies the ideal gas state equation. The total gas content in the matrix system can be calculated by the following formula:

[0099]

[0100] Where φ m is the matrix porosity, dimensionless.

[0101] The diffusion equation of the matrix system is:

[0102]

[0103] S3: Using the finite element method of COMSOL, the gas-water two-phase seepage equation of the fracture system is defined, and combined with the desorption equation and diffusion equation of the matrix system to form a dual-medium two-phase coupled flow equation system;

[0104] The gas-water two-phase flow equation in the fracture system can be written as:

[0105]

[0106] Where m g =ρ g φ f S g , represents the mass of free gas in cracks per volume of coal, kg / m 3 ;m w =ρ w φ f S w , represents the mass of water in the cracks per volume of coal, kg / m 3 ρ g is the gas density, kg / m 3 ρ w is the density of water, kg / m 3 ;S g and S w are gas and water saturations, dimensionless; φ f is the porosity of the fracture, dimensionless; k f is the fracture permeability, mD; μ g and μ w are the dynamic viscosities of methane and water, Pa·s; b1 is the Klinkenberg coefficient, Pa; k rg and k rw are the relative permeabilities of the gas phase and the water phase, dimensionless; p fw is the partial pressure of water in the fracture system.

[0107] The relative permeability of gas and water in a fracture system can be given by the following formula:

[0108]

[0109] Where m is the relative permeability fitting coefficient, dimensionless, and in this study, m = 0.65; S eis the effective water saturation, dimensionless, and can be defined by the following formula:

[0110]

[0111] Where S wr is the residual water saturation, S gr is the residual gas saturation.

[0112] The dual medium two-phase coupled flow equations are:

[0113]

[0114] S4: Considering the desorption hysteresis effect and the gas-water two-phase seepage characteristics, the porosity and permeability changes of the matrix system and the fracture system are calculated respectively, and the physical property change equations of the dual medium are derived;

[0115] Considering the desorption hysteresis, the porosity change of the matrix system can be expressed as follows:

[0116]

[0117] Where S = ε V +(p m / K s )-ε sq , ε V is the volume strain of coal, dimensionless; K s is the bulk modulus of the coal skeleton, GPa; ε sq represents the strain during gas desorption, ε sq =α sg W de , dimensionless; α sg is the matrix expansion coefficient, kg / m 3 , where the subscript '0' represents the initial value of the corresponding variable.

[0118] Considering the desorption hysteresis, the porosity change of the fracture system can be expressed as follows:

[0119]

[0120] Where M is the axial constraint modulus, MPa, and M = E(1-v) / [(1+v)(1-2v)]. Considering the cubic law between the permeability and porosity of porous media, the permeability of the matrix and fracture system can be expressed as:

[0121]

[0122] S5: Using the numerical solver of COMSOL, set the corresponding initial conditions, and use the iterative method to solve the coupled equations in steps S2-S4, and obtain the gas and water production through calculation.

[0123] The formula for calculating gas and water production is:

[0124]

[0125] where Q g and Q w are daily gas and water production, m 3 ; and are the average pressures of gas phase and water phase, MPa; p wf is the bottom hole pressure, MPa; r w is the wellbore radius, m; r e is the drainage radius, m; S is the skin coefficient.

[0126] Steps S1 to S5 have the following assumptions:

[0127] 1) The coal body is divided into two systems: matrix and fracture. Adsorbed gas mainly exists in the matrix system, fractures are the main seepage channels for fluids, and water only exists in the fracture system.

[0128] 2) The strain of the coal body is considered to be infinitesimal;

[0129] 3) Adsorbed gas desorbs from matrix pores and enters the fracture system, subsequently flowing into the wellbore;

[0130] 4) Ignore the effects of gravity and temperature;

[0131] 5) The gas is equivalent to an ideal gas, and its viscosity remains unchanged during the flow process. Water is an incompressible fluid;

[0132] 6) The adsorption of methane in the matrix is saturated adsorption.

[0133] Example 1:

[0134] The data used in this example comes from Well No. 1 in the southern area of the Fanzhuang block in the Qinshui Basin. The calculation parameters are shown in the following table.

[0135] Table 1 Calculation parameters

[0136]

[0137] According to the present invention, the calculation data in Table 1 is used, the time step is set to 1 day, and the total time length is 1000 days, to calculate the daily gas production and daily water production of the production well. Figure 1 The comparison between the daily production of horizontal wells calculated by the method of the present invention and the actual production is shown. Figure 1 It can be seen that the results calculated using the method of the present invention are highly consistent with the actual data, which verifies the accuracy of the method of the present invention.

[0138] This paper provides a numerical simulation method for gas-liquid two-phase flow that considers coalbed methane desorption hysteresis. Combining well logging and experimental data, the desorption hysteresis effect defines the desorption and diffusion equations for the matrix system and the gas-water two-phase flow equations for the fracture system, forming a set of coupled flow equations for the dual-medium two-phase flow. The coupled equations are iteratively solved using COMSOL software, and the gas and water production rates are calculated. This method is simple to operate and can more scientifically and rationally assess the impact of coalbed methane desorption hysteresis on gas well productivity. It also significantly reduces the computational cost and time of gas well numerical simulations, providing important insights for the development of other unconventional gas reservoirs affected by desorption hysteresis.

[0139] The above description does not limit the present invention in any form. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A numerical simulation method for gas-liquid two-phase flow considering coalbed methane desorption hysteresis, characterized in that: The following steps are involved: S1: Use well logging and indoor experiments to obtain production capacity prediction parameters; S2: Establish the computational domain of the CBM dual-medium two-phase model in COMSOL, and define the desorption equation and diffusion equation of the matrix system in combination with the desorption hysteresis effect; S3: Using the finite element method of COMSOL, the gas-water two-phase seepage equation of the fracture system is defined, and combined with the desorption equation and diffusion equation of the matrix system to form a dual-medium two-phase coupled flow equation system; S4: Considering the desorption hysteresis effect and the gas-water two-phase seepage characteristics, the porosity and permeability changes of the matrix system and the fracture system are calculated respectively, and the physical property change equations of the dual medium are derived; S5: Using the numerical solver of COMSOL, set the corresponding initial conditions, and use the iterative method to solve the coupled equations in steps S2-S4, and obtain the gas and water production through calculation; Steps S1 to S5 have the following assumptions: 1) The coal body is divided into two systems: matrix and fracture. Adsorbed gas mainly exists in the matrix system, fractures are the main seepage channels for fluids, and water only exists in the fracture system. 2) The strain of the coal body is considered to be infinitesimal; 3) Adsorbed gas desorbs from matrix pores and enters the fracture system, subsequently flowing into the wellbore; 4) Ignore the effects of gravity and temperature; 5) The gas is equivalent to an ideal gas, and its viscosity remains unchanged during the flow process. Water is an incompressible fluid; 6) The adsorption of methane in the matrix is saturated adsorption.

2. A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis as claimed in claim 1, characterized in that: The capacity prediction parameters in step S1 include: 1) Coal seam parameters Coal seam length, coal seam width, coal seam thickness, coal seam density, matrix initial porosity, matrix initial permeability, fracture initial porosity, fracture initial permeability, coal seam temperature, coal seam initial pressure, coal seam Young's modulus, initial water saturation, irreducible water saturation; 2) Adsorption and desorption parameters Langmuir pressure, Langmuir volume, desorption hysteresis coefficient; 3) Fluid parameters Methane viscosity, methane density, water viscosity, water density, bottom hole flowing pressure, standard temperature.

3. A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis as claimed in claim 1, characterized in that: The desorption equation of the matrix system in step S2 is: W de =W ad +a·exp(b·p m ) Where a is the hysteresis gas content coefficient, m 3 / kg; b is the hysteresis gas pressure coefficient, MPa -1 , coefficients a and b can be fitted in actual adsorption and desorption experiments; p m is the matrix system pressure, MPa; W de is the gas content of coal rock during desorption, m 3 / kg,W ad is the adsorbed gas content in the coal seam, m 3 / kg, which can be defined by the following formula: Where V L is the Langmuir volume constant of the adsorption process, m 3 / kg; P L is the Langmuir pressure constant of the adsorption process, MPa; V L-de is the Langmuir volume constant of the desorption process, m 3 / kg; P L-de is the Langmuir pressure constant of the desorption process, MPa; W re is the residual adsorption capacity, m 3 / kg; The desorption hysteresis pressure is a threshold pressure in the desorption state, which is defined as p q , MPa, can be calculated according to the following formula: Where ρ c is the density of coal rock, kg / m 3 ; V M is the molar volume of methane gas, m 3 / mol; R is the gas molar constant, J / (mol·K); T is the temperature, K According to Fick's diffusion law, the mass conservation equation of gas in the matrix is: where Q m =Dσ c (c m -c f ), represents the mass transfer rate of methane between the matrix and the fracture system, kg / (m 3 ·s);m m is the mass of methane per unit volume of coal matrix, kg / m 3 ; σ c is the coal matrix shape factor, m -2 , defined as L is the crack spacing, m; D is the diffusion coefficient, m 2 / s; c m and c f are the concentrations of methane in the matrix and in the fracture system, kg / m 3 ; According to the ideal gas law: Where M g is the molar mass of methane, kg / mol; p fg is the gas pressure in the fracture system, MPa; considering the desorption hysteresis effect, the mass transfer rate of methane between the matrix and the fracture system is obtained as: The adsorbed gas in the matrix is represented by the Langmuir equation, while the free gas satisfies the ideal gas state equation. The total gas content in the matrix system can be calculated by the following formula: Where φ m is the matrix porosity, dimensionless; The diffusion equation of the matrix system is: According to this equation, the diffusion behavior of gas in the matrix system can be characterized.

4. A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis as claimed in claim 1, characterized in that: The gas-water two-phase seepage equation of the fracture system in step S3 can be written as: Where m g =ρ g φ f S g , represents the mass of free gas in cracks per volume of coal, kg / m 3 ; m w =ρ w φ f S w , represents the mass of water in the cracks per volume of coal, kg / m 3 ; ρ g is the gas density, kg / m 3 ; ρ w is the density of water, kg / m 3 ; S g and S w are gas and water saturations, dimensionless; φ f is the porosity of the fracture, dimensionless; k f is the fracture permeability, mD; μ g and μ w are the dynamic viscosities of methane and water, Pa·s; b1 is the Klinkenberg coefficient, Pa; k rg and k rw are the relative permeabilities of the gas phase and the water phase, dimensionless; p fw is the partial pressure of water phase in the fracture system; The relative permeability of gas and water in a fracture system can be given by the following formula: Where m is the relative permeability fitting coefficient, dimensionless, and in this study, m = 0.65; S e is the effective water saturation, dimensionless, and can be defined by the following formula: Where S wr is the residual water saturation, S gr is the residual gas saturation; The dual medium two-phase coupled flow equations are: The above formula can characterize the two-phase coupled flow law of dual media.

5. The method for numerical simulation of gas-liquid two-phase flow considering coalbed methane desorption hysteresis as claimed in claim 1, characterized in that: Considering the desorption hysteresis in step S4, the porosity change of the matrix system can be expressed as follows: Where S = ε V +(p m / K s )-ε sq , ε V is the volume strain of coal, dimensionless; K s is the bulk modulus of the coal skeleton, GPa; ε sq represents the strain during gas desorption, ε sq =α sg W de , dimensionless; α sg is the matrix expansion coefficient, kg / m 3 , where the subscript '0' represents the initial value of the corresponding variable; Considering the desorption hysteresis, the porosity change of the fracture system can be expressed as follows: Where M is the axial constraint modulus, MPa, and M = E(1-v) / [(1+v)(1-2v)]. Considering the cubic law between the permeability and porosity of porous media, the permeability of the matrix and fracture system can be expressed as: This set of equations is the calculation formula for the permeability of the matrix and fracture system.

6. A gas-liquid two-phase flow numerical simulation method considering coalbed methane desorption hysteresis as claimed in claim 1, characterized in that: The formula for calculating the gas and water production in step S5 is: where Q g and Q w are daily gas and water production, m 3 ; and are the average pressures of gas phase and water phase, MPa; p wf is the bottom hole pressure, MPa; r w is the wellbore radius, m; r e is the drainage radius, m; S is the skin coefficient.

Images