Sea-land transition phase shale gas productivity evaluation method based on p-EDFM

Abstract

The invention discloses a sea-land transition phase shale gas productivity evaluation method based on p-EDFM, and the method comprises the steps: carrying out the layering processing of a matrix grid according to the longitudinal heterogeneity of sea-land transition phase shale, and introducing physical property parameters and corresponding apparent permeability models for different rocks. Dimension reduction processing is carried out on a fracture surface according to parameters of an artificial fracture, pressure exchange between grids is represented by utilizing a fluid channeling coefficient, flow exchange between different grids is further calculated, and then output of each lithology can be obtained through a horizon tracking mode. The combination of the multi-continuous medium model and the embedded discrete fracture model respectively overcomes the defects of the self model, and can better reflect the heterogeneous formation characteristics of microfractures widely distributed in the shale gas reservoir, large-scale artificial fractures and transition-phase shale gas reservoirs. The method provided by the invention not only can be used for evaluating the productivity of the sea-land transition phase shale gas, but also can be popularized and used for calculating the productivity of unconventional oil and gas reservoirs, and has a relatively high application value.

Description

Technical Field

[0001] The invention relates to the technical field of shale gas exploration and development, and in particular to a method for evaluating the productivity of marine-continental transitional phase shale gas based on p-EDFM. Background Art

[0002] The development of marine-continental transitional shale gas reservoirs faces many challenges, including vertical multi-lithology stacking, poor lateral continuity, thin thickness of favorable layers, and differentiated shale gas enrichment caused by the variability of sedimentary environment. These factors cause large fluctuations in shale gas well production and make it difficult to predict gas well production capacity. At present, multi-stage fracturing horizontal well technology is usually used to develop shale gas reservoirs. After hydraulic fracturing, a complex fracture network structure is formed. In the production capacity evaluation, it is necessary to accurately characterize the fracture morphology and physical properties to ensure the accuracy and reliability of the production capacity evaluation.

[0003] Fracture characterization methods are mainly divided into explicit and implicit methods. The implicit fracture characterization method regards a large number of discrete and disordered reservoir fractures as continuous media. In 1963, Warren and Root proposed a fracture-matrix dual medium model, setting porosity and permeability parameters for the two systems respectively. In 1976, Kazemi modified the Warren-Root model by constructing a multiphase fluid conduction function that considers capillary forces and matrix shape factors. In 1985, Pruess et al. further proposed a multi-medium model, which made a more detailed division of the matrix system in the dual medium model. These models have become classic methods for implicit fracture characterization and simulation of fractured reservoirs. However, for horizontal wells with complex fracture network structures, these methods have large errors from the actual situation. In view of the above situation, researchers tend to use explicit fracture characterization methods. The explicit characterization method models each fracture separately, and the commonly used simulation method is the discrete fracture model (DFN). The DFN model can better characterize the distribution properties of fractures in the reservoir, but there is still a large distortion in the description of the matrix system. To this end, in 2000, Lee et al. proposed the embedded discrete fracture model (EDFM), which constructs two systems, the fracture and the matrix, and uses a connection pair to communicate the two. EDFM inherits the advantages of the DFN model and can effectively construct complex fracture systems. At the same time, it does not require mesh fitting of fracture morphology or local encryption processing, thereby significantly reducing the computational cost. In 2017, Tene et al. further proposed the p-EDFM model, which improved the shortcomings of EDFM in simulating low-permeability fractures.

[0004] At present, these methods are usually used to evaluate the production capacity of the entire gas reservoir, but because there are multiple layers in the marine-continental transitional shale, it is difficult to accurately track the production capacity changes of each layer. Therefore, it is urgent to carry out targeted research to optimize the calculation method of the existing model to achieve accurate tracking of the output of each layer of the marine-continental transitional shale. Summary of the invention

[0005] Aiming at the problem that the existing gas reservoir productivity evaluation methods are difficult to accurately track the productivity changes of each layer of marine-continental transitional phase shale, the present invention provides a marine-continental transitional phase shale gas productivity evaluation method based on p-EDFM.

[0006] The p-EDFM-based method for evaluating the productivity of shale gas in the marine-continental transition phase provided by the present invention comprises the following steps:

[0007] S1: Obtain basic parameters of reservoirs, fractures, and wells;

[0008] The parameters obtained include the geometric parameters of the matrix, fractures and wells, the basic physical properties of the matrix and fractures, and the gas-water relative permeability curves and capillary force curves of the matrix and fractures, as well as the gas-water two-phase PVT data, wellbore parameters and skin factors.

[0009] S2: Use orthogonal grids to divide the matrix and use embedded discrete fracture models to couple natural fractures, artificial fractures and matrix grids.

[0010] The matrix and the crack are coupled by connecting pairs. Specifically, there are three types of connecting pairs:

[0011] The first one is when the fracture grid intersects with the matrix grid, the connection pair between the fracture grid and the matrix grid is calculated as follows:

[0012]

[0013] in,

[0014]

[0015] In the formula, subscript i represents the fracture grid, subscript j and k represent the matrix grid; G mf is the connection pair between the fracture grid and the matrix grid, m 3 ; G m is the connection pair between matrix grids, m 3 ; G f is the connection pair between crack grids, m 3 ; A mf is the projected area of ​​the fracture grid on the matrix grid, m 2 ;k m is the matrix permeability, μm 2 ;k f is the crack permeability, μm 2 ;d mf is the characteristic distance between the matrix grid and the fracture grid, m; is the fracture; A f is the crack area, m 2 ; is the crack opening of crack grid i, m; xn is the normal distance from the volume element to the crack, m; dV is the volume element of the grid block, m 3 ; V is the grid volume, m 3 .

[0016] The second method is to calculate the connection pairs of different cracks when they intersect:

[0017]

[0018] in,

[0019]

[0020] In the formula, G ff is the connection pair when different crack grids intersect, m 3 ; G ffi is the connection pair of the i-th crack grid, m 3 β c is the shape factor, dimensionless; k fi is the permeability of the i-th fracture grid, μm 2 ;ω fi is the opening of the i-th crack, m; L int is the length of the intersection of the two cracks, m; d fi is the average distance from the i-th crack grid to the intersection line, m; dS i is the area element of the i-th crack grid, m 2 ; S i is the area of ​​the i-th crack grid, m 2 ;x n is the normal distance from the area element to the crack intersection line, m.

[0021] The third method is to calculate the adjacent fracture grid connection pairs within the same fracture:

[0022]

[0023] in,

[0024]

[0025] In the formula, G f is the connection pair when adjacent crack grids intersect, m 3 ; G fi is the connection pair of the i-th crack grid, m 3 β c is the shape factor, dimensionless; k fi is the permeability of the i-th fracture grid, μm 2 ; A f is the common area of ​​adjacent cracks, m 2;d fi is the average distance from the i-th crack grid to the intersection line, m.

[0026] S3: The fracture grid where the wellbore and the artificial fracture intersect is set and processed as the wellbore grid;

[0027] When the fracture intersects with the wellbore grid, the fracture grid intersecting with the wellbore trajectory is processed as the wellbore grid, and the calculation method is:

[0028]

[0029] in,

[0030]

[0031] Where WI is the connection coefficient between fracture and wellbore, m 3 ; Δθ is the angle between the wellbore trajectory and the fracture center, dimensionless; k f is the crack permeability, μm 2 ;ω f is the crack opening, m; r e is the effective wellbore radius, m; r w is the actual wellbore radius, m; s is the skin coefficient, dimensionless, L f is the crack length, m; h f is the height of the crack segment, m.

[0032] S4: Calculate the apparent permeability of the matrix, using the corresponding calculation method according to the different matrix types:

[0033] (1) For shale layers, the calculation formula for apparent permeability is as follows:

[0034]

[0035] in,

[0036]

[0037] In the formula, ρ g is the gas density in the porous medium, kg / m 3 ;μ g is the gas viscosity in the porous medium, Pa·s; k 0 is the permeability of the porous medium, μm 2 ; Kn is the Knudsen coefficient, dimensionless; λ is the molecular mean free path, m; D is the diameter of the roar channel, m; M g is the molar mass of gas molecules, g / mol; P is the pressure in the porous medium, MPa; Z is the compression factor, dimensionless; R is the gas constant, whose value is 8.31×10 -3, J / (mol·K); T is the temperature of the porous medium, K; α is the gas rarefaction effect coefficient, dimensionless; b is the slip factor, dimensionless; C g is the gas compressibility coefficient, MPa -1 ;D sur is the gas surface diffusion coefficient, m 2 / s;P sc is the atmospheric pressure under standard conditions, MPa; T sc is the temperature under standard conditions, K; ρ r is the rock density, g / cm 3 ;P L is the Langmuir pressure, MPa; V L is the Langmuir volume, m 3 .

[0038] (2) For sandstone reservoirs, the calculation formula for apparent permeability is as follows:

[0039]

[0040] (3) For coal rock reservoirs, the calculation formula for apparent permeability is as follows:

[0041]

[0042] S5: Based on the continuity equation and the equation of motion, the seepage equations for the matrix, natural fractures and artificial fractures are derived:

[0043] The gas phase percolation equation in the matrix is:

[0044]

[0045] in,

[0046]

[0047] In the formula, k mg is the absolute permeability of the matrix, m 2 ;k app,m is the apparent permeability of the matrix, dimensionless; B mg is the gas phase volume coefficient of the matrix, dimensionless; μ mg is the viscosity of the gas phase, Pa·s; P mg is the pressure of the matrix system, MPa; ρ mg is the matrix gas phase density, kg / m 3 ; g is the acceleration due to gravity, 9.8m / s 2 ;h m represents the height of the matrix system, m; q fm,g is the gas phase flow exchange between the fracture system and the matrix system, m 3 / s;φm is the porosity of the matrix, dimensionless; S mg is the matrix gas phase saturation, dimensionless; ρ r is the matrix rock density, kg / m 3 ;q dec is the matrix analysis flow rate, m 3 / kg·s; V e is the amount of adsorbed gas when adsorption reaches equilibrium, m 3 / kg.

[0048] The gas phase seepage equation in other systems (natural fractures and artificial fractures) is:

[0049]

[0050] Where, subscript s is the system type, f is natural fracture, and F is artificial fracture; q exc,g is the gas phase flow exchange between this system and other systems, m 3 / s;k s,g is the absolute permeability of the system, m 2 ;k s,rg is the gas phase relative permeability of the system, dimensionless; B s,g is the gas phase volume coefficient of the system, dimensionless; μ s,g is the viscosity of the gas phase of the system, Pa·s; P s,g is the pressure of the system, MPa; ρ s,g is the gas phase density of the system, kg / m 3 ;h s is the height of the system, m; φ s is the porosity of the system, dimensionless; S s,g is the gas phase saturation of the system, dimensionless.

[0051] The liquid phase flow equation is:

[0052]

[0053] Where, subscript s is the system type, m is the matrix, f is the natural fracture, and F is the artificial fracture; k s,w is the absolute permeability of the system, m 2 ;k s,rw is the liquid relative permeability of the system, dimensionless; B s,w is the liquid phase system number of the system, dimensionless; μ s,w is the liquid viscosity of the system, Pa·s; P s,w is the pressure of the system, MPa; q exc,w is the liquid phase flow exchange between this system and other systems, m 3 / s;ρ w is the density of formation water, kg / m3 ;h s is the height of the system, m; φ s is the porosity of the system, dimensionless; S s,g is the water phase saturation of the system, dimensionless.

[0054] S6: The seepage equation is discretized using the finite volume method, and the discretization result is solved using the fully implicit Newton-Raphson iteration method to establish the final mathematical model.

[0055] The results of discretizing the seepage equation using the finite volume method are as follows:

[0056]

[0057] in,

[0058]

[0059] ψ i,s =P i,s -ρ i,s gh i,s , ψ j,s =P j,s -ρ j,s gh j,s

[0060] Where, subscripts i and j represent adjacent grids i and j; subscript s represents the system type, m represents the matrix, f represents the natural fracture, and F represents the artificial fracture; subscript l represents the fluid type, w represents the water phase, and g represents the gas phase; T ij,s is the conductivity of grid i and grid j, kg / Pa·s; A i,j is the common area between grid i and grid j, m 2 ;k l,s is the fluid permeability of the system, m 2 ;μ l,s is the fluid viscosity of the system, Pa·s; B l,s is the fluid saturation of the system, dimensionless; is the vector from the center of grid i to the common surface, dimensionless; is the vector from the center of grid j to the common surface, dimensionless; is the normal vector of the common surface, dimensionless; ψ i,s is the potential of grid i, MPa; ψ j,s is the potential of grid j, MPa; n is the current time step, dimensionless; is the mass density in grid i of different systems at n+1 time steps, dimensionless; is the flow rate in grid i of different systems at n+1 time steps, m 3 / s; ΔVi is the volume of grid i; Δt is the time step, s; P i,s is the pressure of grid i in different systems, MPa; ρ i,s is the fluid density in the grid i of different systems, kg / m 3 ;h i,s is the height of grid i in different systems, m.

[0061] Furthermore, the numerical equation obtained by solving is as follows:

[0062]

[0063] In the formula, is the gas phase flow rate at step n+1, m 3 / s; is the liquid flow rate at step n+1, m 3 / s; Respectively represent the coefficients of gas phase pressure change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of liquid phase pressure change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of gas saturation change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of liquid saturation change of the i-th and j-th grids at the current time step; g represents the gas phase, w represents the liquid phase; δP i Indicates the pressure change; δS i Indicates the change in saturation; d g is the coefficient of gas phase flow exchange of the i-th and j-th grids at the current time step; d w is the coefficient of liquid flow exchange of the i-th and j-th grids at the current time step.

[0064] S7: Use the mathematical model established in step S6 to evaluate the productivity of marine-continental transitional shale gas, record the flow between various systems at each time step, and the systems include three systems: natural fractures, artificial fractures and matrix; further calculate the output of each layer.

[0065] The calculation formula for the flow of each system is as follows:

[0066] ∑Q=∑λGΔP ij

[0067] In the formula, Q is the flow rate between different systems, m 3 / s; λ is the fluid velocity, m 2 ·s / kg; G is the connecting pair; ΔP ij is the pressure difference between the two grids, MPa.

[0068] Furthermore, when calculating the production in step S7, the overall flow rate from the hydraulic fracturing cracks to the wellbore is calculated, and the production of the specific layer cannot be accurately determined. Therefore, it is necessary to calculate the flow rates from the matrix and microcracks to the artificial fractures respectively, and determine the specific production layer according to the grid numbers corresponding to the matrix and microcracks.

[0069] The flow rate of matrix to artificial fracture is:

[0070] Q l =G mF λ l (P m -P F )

[0071] Where, subscript l is the fluid type, gas phase is g, liquid phase is w; Q l is the fluid flow rate from the matrix to the artificial fracture, m 3 / s; G mF is the connection pair between matrix and artificial crack, dimensionless; λ l is the fluid velocity, m 2 ·s / kg; P m is the pressure of the matrix grid, MPa; P F is the pressure of the artificial fracture grid, MPa.

[0072] The flow rate from natural fractures to artificial fractures is:

[0073] Q l =G fF λ l (P f -P F )

[0074] In the formula, G fF is the connection pair of natural fractures and artificial fractures, dimensionless; P f is the pressure of the natural fracture grid, MPa.

[0075] Compared with the prior art, the present invention is beneficial in that:

[0076] (1) Based on the p-EDFM model, the present invention is based on the multiple continuous medium model and the embedded discrete fracture model. By dividing the target stratum into matrix grid, microfracture grid and artificial fracture grid, the accurate characterization of the complex fracture system is achieved. In view of the vertical heterogeneity of the marine-continental transitional shale, a layered processing method is adopted to layer the matrix grid, and the physical property parameters of different rock types and the corresponding apparent permeability model are introduced, so that the description of the fracture morphology and seepage characteristics is more accurate.

[0077] (2) The fracture surface is reduced in dimension according to the parameters of the artificial fractures, and the crossflow coefficient is used to characterize the pressure exchange between grids. The flow exchange between different grids is further calculated, which optimizes the exchange process of pressure and flow between grids. Combined with the layer tracking technology, the output of each lithology can be accurately obtained, thus solving the technical problem that the traditional model cannot accurately track the layer production capacity.

[0078] (3) The combination of the multiple continuous medium model and the embedded discrete fracture model overcomes the shortcomings of each model. The combination of the two can better reflect the widely distributed microcracks, large-scale artificial fractures and heterogeneous stratigraphic characteristics of transitional shale gas reservoirs.

[0079] (4) The method of the present invention can not only improve the accuracy of the production capacity evaluation of marine-continental transitional phase shale gas reservoirs, but also has low calculation cost and high model adaptability. It is suitable for the development and evaluation of unconventional oil and gas reservoirs and has wide application value.

[0080] Other advantages, objectives and features of the present invention will be embodied in part through the following description, and in part will be understood by those skilled in the art through study and practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 The present invention is a flow chart of the method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM.

[0082] Figure 2 The following is a comparison of the effects of the evaluation method of the present invention and commercial software.

[0083] Figure 3 It is a comparison chart of the output of each layer in the embodiment. DETAILED DESCRIPTION

[0084] The preferred embodiments of the present invention are described below in conjunction with the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0085] The data used in this example are from the public literature (Numerical simulation study on combined mining of marine-continental transitional shale reservoirs, Chen Xuezhong, 2024), and the model parameter settings are shown in Table 1.

[0086] Table 1 Model parameter settings

[0087] Model parameters Coal seam Shale formation Dense sandstone layer Matrix porosity / % 6 3 4 <![CDATA[Matrix permeability / 10 -3 um 2 > 0.1 0.001 0.01 <![CDATA[Langmuir volume / (m 3 / t)]]> 13.98 2.35 —— Langmuir pressure / MPa 2.5 3.6 —— Matrix water saturation 0 0 40 Hydraulic fracture opening 0.01 0.01 0.01 Hydraulic fracture permeability 1 1 1 Hydraulic fracture porosity 0.01 0.01 0.01

[0088] according to Figure 1The process shown in the figure is used to derive the differential equation of gas-water two-phase seepage based on the coupled p-EDFM and multiple continuous medium models. The mathematical model is obtained using the method described in the present invention, and the relevant numerical simulator is written using the python computer language. A gas reservoir model of 1800m×560m×40m (length×width×height) is established, with a grid size of 40m×40m×10m. Eight hydraulic fractures are set evenly distributed on both sides of the horizontal wellbore, and the simulated production time is 800 days. The numerical simulator is verified using the commercial software CMG. Figure 2 As shown, it can be seen that the model established by the method of the present invention has a good fitting effect with the commercial software, which shows that the model obtained by the present invention can accurately characterize the seepage conditions of each layer, and the calculation results of the numerical model obtained are reasonable and reliable. Figure 3 This is a comparison chart of the output of each layer, from which we can see the contribution rate of each lithology output, which is helpful for studying the complex seepage mechanism of marine-continental transitional shale gas reservoirs.

[0089] The above description is only a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Although the present invention has been disclosed as a preferred embodiment as above, 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 modification, equivalent change and modification 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 still falls within the scope of the technical solution of the present invention.

Claims

1. A method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM, characterized in that: The steps include: S1: Obtain basic parameters of reservoirs, fractures, and wells; S2: Use orthogonal grids to divide the matrix and use embedded discrete fracture models to couple natural fractures, artificial fractures and matrix grids; S3: The fracture grid where the wellbore and the artificial fracture intersect is set and processed as the wellbore grid; S4: Calculate the apparent permeability of the matrix, using the corresponding calculation method according to the different matrix types: (1) For shale layers, the calculation formula for apparent permeability is as follows: in, In the formula, ρ g is the gas density in the porous medium, kg / m 3 ;μ g is the gas viscosity in the porous medium, Pa·s; k0 is the permeability of the porous medium, μm 2 ; Kn is the Knudsen coefficient, dimensionless; λ is the molecular mean free path, m; D is the diameter of the roar channel, m; M g is the molar mass of gas molecules, g / mol; P is the pressure in the porous medium, MPa; Z is the compression factor, dimensionless; R is the gas constant, whose value is 8.31×10 -3 , J / (mol·K); T is the temperature of the porous medium, K; α is the gas rarefaction effect coefficient, dimensionless; b is the slip factor, dimensionless; C g is the gas compressibility coefficient, MPa -1 ;D sur is the gas surface diffusion coefficient, m 2 / s;P sc is the atmospheric pressure under standard conditions, MPa; T sc is the temperature under standard conditions, K; ρ r is the rock density, g / cm 3 ;P L is the Langmuir pressure, MPa; V L is the Langmuir volume, m 3 ; (2) For sandstone reservoirs, the calculation formula for apparent permeability is as follows: (3) For coal rock reservoirs, the calculation formula for apparent permeability is as follows: S5: Derive the seepage equations for matrix, natural fractures and artificial fractures based on the continuity equation and the equation of motion; S6: The seepage equation is discretized using the finite volume method, and the fully implicit Newton-Raphson iteration method is used to solve the discretization result and establish the final mathematical model; the numerical equation obtained by solving is as follows: In the formula, is the gas phase flow rate at step n+1, m 3 / s; is the liquid flow rate at step n+1, m 3 / s; Respectively represent the coefficients of gas phase pressure change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of liquid phase pressure change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of gas saturation change of the i-th and j-th grids at the current time step; Respectively represent the coefficients of liquid saturation change of the i-th and j-th grids at the current time step; g represents the gas phase, w represents the liquid phase; δP i Indicates the pressure change; δS i Indicates the change in saturation; d g is the coefficient of gas phase flow exchange of the i-th and j-th grids at the current time step; d w is the coefficient of liquid phase flow exchange of the i-th and j-th grids at the current time step; S7: Use the mathematical model established in step S6 to evaluate the productivity of marine-continental transitional shale gas, record the flow between various systems at each time step, and the systems include three systems: natural fractures, artificial fractures and matrix; further calculate the output of each layer.

2. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: In step S6, the result of discretizing the seepage equation using the finite volume method is as follows: in, ψ i,s =P i,s -r i,s gh i,s ,ψ j,s =P j,s -r j,s gh j,s Where, subscripts i and j represent adjacent grids i and j; subscript s represents the system type, m represents the matrix, f represents the natural fracture, and F represents the artificial fracture; subscript l represents the fluid type, w represents the water phase, and g represents the gas phase; T ij,s is the conductivity of grid i and grid j, kg / Pa·s; A i,j is the common area between grid i and grid j, m 2 ;k l,s is the fluid permeability of the system, m 2 ;μ l,s is the fluid viscosity of the system, Pa·s; B l,s is the fluid saturation of the system, dimensionless; is the vector from the center of grid i to the common surface, dimensionless; is the vector from the center of grid j to the common surface, dimensionless; is the normal vector of the common surface, dimensionless; ψ i,s is the potential of grid i, MPa; ψ j,s is the potential of grid j, MPa; n is the current time step, dimensionless; is the mass density in grid i of different systems at n+1 time steps, dimensionless; is the flow rate in grid i of different systems at n+1 time steps, m 3 / s; ΔV i is the volume of grid i; Δt is the time step, s; P i,s is the pressure of grid i in different systems, MPa; ρ i,s is the fluid density in the grid i of different systems, kg / m 3 ;h i,s is the height of grid i in different systems, m.

3. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: The parameters obtained in step S1 include geometric parameters of the matrix, fractures, and wells, basic physical parameters of the matrix and fractures, and gas-water relative permeability curves and capillary force curves of the matrix and fractures, as well as gas-water two-phase PVT data, wellbore parameters, and skin factors.

4. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: In step S2, the matrix and the crack are coupled using connection pairs, including the following three types of connection pairs: the first type, when the crack grid intersects with the matrix grid, the connection pair between the crack grid and the matrix grid; the second type, when different cracks intersect, the connection pair between different cracks; the third type, the connection pair between adjacent crack grids in the same crack.

5. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: In step S3, when the fracture intersects with the wellbore grid, the fracture grid intersecting with the wellbore trajectory is processed as a wellbore grid.

6. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: In step S5, the seepage equations include the gas phase seepage equation in the matrix, the gas phase seepage equation and the liquid phase seepage equation in the natural fracture, and the gas phase seepage equation and the liquid phase seepage equation in the artificial fracture.

7. The method for evaluating the productivity of shale gas in the marine-continental transition phase based on p-EDFM according to claim 1, characterized in that: In step S7, the calculation formula of each system flow is as follows: ∑Q=∑λGΔP ij In the formula, Q is the flow rate between different systems, m 3 / s; λ is the fluid velocity, m 2 ·s / kg; G is the connecting pair; ΔP ij is the pressure difference between the two grids, MPa.

Images