A MINC-EDFM solution method considering isotopic fractionation

Abstract

This invention provides a MINC-EDFM solution method considering isotopic fractionation, relating to the field of fractured reservoir simulation. Specifically, it includes the following steps: establishing a three-medium dual-isotope coupled model and an initial discrete grid; constructing a joint error index driven by isotopic errors; adaptively correcting the transition point of MINC, the fracture neighborhood grid of EDFM, and the local time step based on the joint error index; performing isotopic conservation correction at the end of the macro-time step and outputting the final result. The technical solution of this invention overcomes the limitations of existing technologies where grid refinement, MINC sublayer adjustment, conductivity coefficient updating, and time step control cannot simultaneously serve pressure propagation and δ¹⁸O₅. 13 The accuracy of biisotopic exchange between the C fractionation front and the matrix and fracture is ensured after local refinement and multi-rate advancement. 12 CH4 13 The issue of mass conservation of CH4 and total methane.

Description

Technical Field

[0001] This invention relates to the field of fractured reservoir simulation, and more specifically to a MINC-EDFM solution method that considers isotopic fractionation. Background Technology

[0002] Post-fracture unconventional gas reservoirs typically exhibit multi-scale flow channels encompassing the matrix, natural fractures, and hydraulic fractures. Existing methods for simulating fractured reservoirs mainly include dual-medium models, discrete fracture models, and embedded discrete fracture models. Dual-medium models offer high computational efficiency but struggle to explicitly describe hydraulic fracture geometry; discrete fracture models offer high accuracy but require body-fitted meshes, which become difficult to mesh when fractures are complex. EDFM embeds fractures into the background mesh, avoiding complex body-fitting meshes and reducing the meshing complexity of complex fracture models by calculating the exchange between the matrix and fractures through non-adjacent connections. Yan et al. pointed out that EDFM does not require bedrock meshing based on fracture morphology; the fracture portion is processed according to its intersection with the bedrock mesh, thus improving the simulation efficiency of complex fractures.

[0003] The closest existing technologies can be summarized as either the "fixed MINC-EDFM numerical simulation method for fractured gas reservoirs" or the "pressure-driven adaptive EDFM method." These existing technologies have the following problems: 1. The pressure-driven adaptive EDFM method primarily targets the pressure or flow field for adaptive correction, and cannot guarantee the accuracy of methane carbon isotope fractionation simulations. For isotope fractionation problems, the spatial variation of δ¹³C is not always synchronized with the pressure gradient. In the later stages of production, the pressure field may have become relatively flat, but adsorbed gas desorption, diffusion replenishment, and exchange between the matrix and fractures can still cause the δ¹³C front to continue advancing. If the mesh refinement is still determined solely by the pressure gradient, it is easy to encounter problems where the pressure calculation results are reasonable, but the δ¹³C curve is numerically dissipated, and the fractionation stage is distorted.

[0004] 2. Existing fixed MINC methods typically pre-determine the number of nested layers, which cannot be dynamically adjusted based on the intensity of isotopic exchange between the matrix and fractures. When isotopic exchange in the fracture neighborhood is intense, a fixed number of layers may be insufficient to describe the unsteady-state recharge differences between ¹²CH4 and ¹³CH4; conversely, retaining too many sublayers when the exchange process slows down leads to ineffective calculations. Existing pressure-driven adaptive methods also cannot determine this. 12 CH4 and 13 Is the local mass balance of CH4 disrupted by discrete errors?

[0005] 3. Existing EDFM local refinement or conductivity corrections primarily focus on the pressure conduction capacity between the matrix and fractures, lacking simultaneous dual-isotope corrections for convective and diffusion exchange. For methane carbon isotope fractionation, ¹²CH₄ and ¹³CH₄ differ in diffusion coefficients, adsorption or desorption responses, and migration capabilities. Correcting only the total gas flux would result in a lack of reliable component conservation basis for δ¹³C calculations.

[0006] Therefore, there is a need for a method that enables mesh refinement, MINC sublayer adjustment, conduction coefficient update, and time step control to simultaneously serve pressure propagation and delta propagation. 13 The accuracy of biisotopic exchange between the C fractionation front and the matrix and fracture is ensured after local refinement and multi-rate advancement. 12 CH4 13 The MINC-EDFM solution method for conserving the mass of CH4 and total methane, considering isotopic fractionation. Summary of the Invention

[0007] The main objective of this invention is to provide a MINC-EDFM solution method that considers isotopic fractionation, thereby addressing the limitations of existing technologies that cannot simultaneously enable mesh refinement, MINC sublayer adjustment, conductivity coefficient updates, and time step control to serve pressure propagation and δ¹⁸O. 13 The accuracy of biisotopic exchange between the C fractionation front and the matrix and fracture is ensured after local refinement and multi-rate advancement. 12 CH4 13 The issue of mass conservation of CH4 and total methane.

[0008] To achieve the above objectives, this invention provides a method for solving the MINC-EDFM problem considering isotopic fractionation, specifically including the following steps: S1. Establish a three-medium dual-isotope coupling model and an initial discrete grid.

[0009] S2, constructing a joint error index driven by isotope error.

[0010] S3, based on the joint error index, adaptively corrects the transition point of the multi-interaction continuous medium model MINC, the crack neighborhood grid of the embedded discrete crack model EDFM, and the local time step.

[0011] S4 performs isotope conservation correction at the end of the macro time step and outputs the final result.

[0012] Furthermore, step S1 specifically includes the following steps: S1.1, the reservoir is divided into matrix media. Microcracked media and hydraulic fracture media .

[0013] S1.2, two isotopic components in methane Establish the component conservation equations separately, and characterize them uniformly as follows: ; in, as medium Middle Isotope Components Storage capacity per unit volume; Isotopic components In the medium Total flux in; For inter-media exchange items; ∇⋅ is the well source term; ∇⋅ is the divergence operator.

[0014] S1.3, Storage Items in the Matrix Medium Written as: .

[0015] Reserves in microfractured and hydraulically fractured media , Write them as follows: ; ; in, , , They are matrix media Microcracked media and hydraulic fracture media porosity; , , They are matrix media Microcracked media and hydraulic fracture media Gas phase saturation; This refers to the gas phase density. isotopic components in the gas phase mole fraction; Bulk density; This refers to the reserves of adsorbed isotopes.

[0016] Furthermore, step S1 also includes the following steps: S1.4, the matrix region is discretized using conventional volumetric meshes, the hydraulic fractures are discretized using EDFM embedded discretization, and the microfractures are characterized using dual-medium continuous characterization; for the exchange zone between the matrix and the fractures, at least one MINC transient nested sublayer is pre-placed near the fractures to describe the early unsteady recharge process.

[0017] Furthermore, step S2 specifically includes the following steps: S2.1, Constructing the pressure gradient index For any control body ,definition: ; in, To control the length of the volume feature; To control the average pressure of the body; To prevent scalar pressure from an excessively small denominator, For control body Internal pressure gradient magnitude.

[0018] S2.2, Definition Gradient index : ; in, The methane carbon isotope value of the current control body; It is an isotopic characteristic scale; It is a small positive number.

[0019] S2.3, for control volumes that are not adjacent to fracture elements, defines the isotopic exchange residual index between the matrix and the fracture. for: ; in, For time step Time control body Internal isotope components reserves, For the current time step, For control body A boundary surface, For control body The set of all boundary surfaces, To control the component flux at each interface of the volume; It serves as a source and sink for exchange between the matrix and the cracks.

[0020] Define the isotopic exchange residual index between the matrix and the fracture. : ; in, To prevent tiny positive numbers with a denominator of zero, , Control body middle 12 CH4 and 13 Exchange residuals between the matrix and cracks of CH4.

[0021] S2.4, after making the three types of indicators dimensionless, define the joint error index. : ; in, , , , These are the weights for pressure, isotopes, and exchange residuals, respectively. They are dimensionless versions , , .

[0022] Furthermore, step S3 specifically includes the following steps: S3.1, For matrix mesh blocks near cracks, the number of transient nested layers is determined by the joint error index. and transition point .

[0023] S3.2, for satisfying The control body performs local refinement; for several consecutive macro time steps, it satisfies... The control body performs local coarsening, where , These are the refinement threshold and the coarsening threshold, respectively.

[0024] S3.3, For different regions, the local time step is determined based on the local spatial scale, convection velocity, pressure diffusion rate, and isotope diffusion rate. : ; in, These are the time step control coefficients for convection, pressure diffusion, and isotope diffusion, respectively. Local velocity; To control the length of the volume feature; To prevent a small positive number with a zero denominator due to a zero velocity; For control body Pressure diffusion rate, For control body The effective diffusion coefficient of isotopes.

[0025] S3.4 divides the computational domain into coarse time step region, medium time step region, and fine time step region, and performs sub-loop iteration within a macro time step in the fine time step region; for the fine time step region, based on the local time step size... Progressing step by step, when the residuals of the pressure equation within a certain sub-time step are... 12 CH4 13When all CH4 mass-conserving residuals meet the preset nonlinear iterative convergence condition, the calculation result of that sub-time step is accepted and the process proceeds to the next sub-time step; if the convergence condition is not met, the local time step size is reduced and the sub-time step is recalculated; when the cumulative advancement time of the local sub-cycle reaches the macro time step... At this point, the sub-cycle in this region stops, and the dual isotope conservation correction begins at the end of the macro time step.

[0026] Furthermore, step S3.1 specifically includes the following steps: S3.1.1, when the control body that is not adjacent to the crack element satisfies: or At that time, the MINC sublayer is repartitioned, and the transition point is... Moving inwards into the matrix; among which, To exchange residual refinement thresholds, This is an index of isotopic exchange residuals between the matrix and the cracks.

[0027] When the adjacent control bodies of the crack are continuous Each time step satisfies: and At that time, the outer MINC sublayers are merged, and the transition point is pushed back towards the crack; among them, To exchange residual coarsening thresholds.

[0028] S3.1.2, Define the pressure diffusion length and isotope diffusion length : , ; in, For pressure diffusion rate, For the medium permeability, These are porosity, overall compressibility, and fluid viscosity, respectively. is the effective diffusion coefficient of the isotope.

[0029] The thickness of the innermost nested layer is taken as : ; in, This is the proportionality coefficient. This is a function that takes the minimum value.

[0030] S3.1.3, each sub-layer expands outward in a geometric progression: ; in, The layer thickness expansion factor is... For the first MINC sublayer thickness, This represents the total number of transient MINC sublayers retained in the current matrix block.

[0031] For transition points When the outer perimeter of a certain layer simultaneously meets the following two conditions, the outer perimeter of that layer is merged into an equivalent equilibrium region, and transient sublayers are no longer retained. The conditions are: and ; in, and These are small positive numbers used to prevent the denominator from being zero. and These are the allowable error thresholds for the pressure balance criterion and the isotope balance criterion, respectively. The MINC sublayer number corresponding to the candidate transition point; and The first The mean pressure and mean carbon isotope value of the MINC sublayer. and These are the volume-weighted average pressure and volume-weighted average carbon isotope value within the candidate equivalent equilibrium region, respectively. This is the function for finding the maximum value.

[0032] S3.1.4, MINC interlayer component exchange conductivity Calculated using equivalent resistance: ; in, Equivalent permeability of components; For component fluidity; For the exchange cross-sectional area.

[0033] Furthermore, after performing local refinement or local coarsening, the non-adjacent connectivity conductivity between the matrix and the fracture is recalculated; isotopic composition Exchange flux Written as: ; in: , ; in, Isotopic components The matrix-to-microcrack convection exchange conductivity, Isotopic components The diffusion exchange conductivity between the matrix and the crack, This represents the exchange area between the matrix and the crack. Equivalent normal distance; Normal penetration rate; These represent the component concentrations of the matrix medium and the microcrack medium, respectively. The effective diffusion coefficient; To control the average pressure of the body, The average pressure of the microcrack element connected to the control volume. Isotopic components The equivalent mobility; For matrix porosity; This represents the gas phase saturation of the matrix.

[0034] Furthermore, the time progression termination loop condition in step S3.4 is: .

[0035] That is, the cumulative advancement time of the local sub-loop reaches the macro time step. When this happens, the sub-loop in that region stops; Nonlinear iterative convergence condition: ; ; in, For the residuals of the pressure equation, , They are respectively 12 CH4 and 13 The mass conservation residual of CH4 and The first After the first nonlinear iteration, the computational region 12 CH4 and 13 Storage vector of CH4.

[0036] Furthermore, step S4 specifically includes the following steps: S4.1, Define the total residual of components within the macro time step: , ; in, , Each within a macro time step 12 CH4 13 The overall mass conservation residual of CH4 , The calculation region at the start of the macro time step 12 CH4 13 Total reserves of CH4 , Calculate the region at the end of the macro time step respectively 12 CH4 13 Total reserves of CH4 These are the macro time steps within which12 CH4 13 The net source and sink of CH4 are defined as positive when entering the calculation region and negative when leaving the calculation region.

[0037] S4.2, in the high-error local region Internally, solve the minimum correction problem: ; in, As weight.

[0038] Satisfy constraints: , ; in, , This is the local element quality correction amount.

[0039] S4.3, after correction, recalculate the isotope ratios and δ. 13 C, output pressure field, gas production rate 12 CH4 and 13 CH4 component distribution, exchange flux between matrix and fracture and δ 13 C is distributed in space and time.

[0040] The present invention has the following beneficial effects: 1. This invention constructs an adaptive criterion for MINC–EDFM based on isotope fractionation error. Instead of solely relying on pressure or saturation gradients for mesh correction, this invention constructs an adaptive error index by combining the pressure gradient, δ¹³C gradient, and matrix-fracture dual isotope exchange residuals. The pressure gradient characterizes the intensity of flow disturbance, and δ¹³C gradient... 13 The C-gradient is used to identify the isotopic fractionation front, and the matrix-fracture isotopic exchange residual is used to determine... 12 CH4 and 13 The local conservation error of CH4 during inter-medium exchange. The core of this invention lies in incorporating isotopic fractionation information itself into the correction basis of numerical grids and discrete structures, rather than treating isotopes merely as the output result after calculation.

[0041] 2. This invention proposes an adaptive adjustment method for MINC sublayers and transition points based on isotope exchange residuals. This invention addresses the unsteady matrix-fracture exchange process and proposes a method based on... 12 CH4 and 13The local exchange residuals of CH4 dynamically adjust the number of MINC sublayers and the location of transition points. When the isotopic exchange residuals are large in the near-fracture region of the stratum, the number of MINC sublayers is increased to refine the unsteady replenishment process of the matrix to the fracture; when both pressure and δ¹³C tend to level off, the outer sublayers are merged to form an equivalent equilibrium zone. The core of this invention is that the nested MINC structure is no longer a fixed number of layers, but a dynamic discrete structure driven by dual isotopic exchange errors.

[0042] 3. This invention proposes a method for local refinement, multi-rate advancement, and conservation correction of EDFM based on dual isotope conservation. During the local refinement of the fracture neighborhood in EDFM, this invention simultaneously updates the matrix-fracture convection exchange coefficient and diffusion exchange coefficient, and applies different time steps to different regions based on a joint error index. After the macro-time step, the method is further refined... 12 CH4 and 13 CH4 is subjected to mass conservation correction, and then δ is calculated from the ratio of the two. 13 C. The core of this invention lies in the fact that local mesh refinement, time step adjustment, and conservation correction all revolve around the conservation of dual isotope mass, thus avoiding spurious δ values ​​caused by numerical discretization errors. 13 C fractionation response. Attached Figure Description

[0043] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1 A flowchart of a method for solving MINC-EDFM considering isotopic fractionation according to the present invention is shown.

[0044] Figure 2 The diagram illustrates the variation of produced gas δ¹³C1 with production time and a comparison of the prediction results of the method of this invention with those of existing methods. Detailed Implementation

[0045] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. 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.

[0046] like Figure 1 The MINC-EDFM solution method shown includes the following steps, taking into account isotopic fractionation: S1. Establish a three-medium dual-isotope coupling model and an initial discrete grid.

[0047] S2, constructing a joint error index driven by isotope error.

[0048] S3, based on the joint error index, adaptively corrects the transition point of the multi-interaction continuous medium model MINC, the crack neighborhood grid of the embedded discrete crack model EDFM, and the local time step.

[0049] S4 performs isotope conservation correction at the end of the macro time step and outputs the final result.

[0050] This invention provides a MINC-EDFM solution method considering isotopic fractionation. This method addresses numerical simulation problems of unconventional gas reservoirs that include three types of media: matrix, microfractures, and hydraulic fractures. 12 CH4 and 13 Based on the independent conservation solution of CH4, a joint error index driven by isotope error is constructed using the pressure gradient, δ¹³C gradient, and matrix-fracture isotope exchange residual. This index adaptively adjusts the number of nested MINC layers and transition point positions on the matrix side, the local grid scale of the EDFM in the fracture neighborhood, and the local time step. Isotope conservation correction is performed at the end of the macro time step, thereby improving the solution efficiency while ensuring the synchronization accuracy of the pressure field, isotope front, and matrix replenishment process.

[0051] Specifically, step S1 includes the following steps: S1.1, the reservoir is divided into matrix media. Microcracked media and hydraulic fracture media Among them, matrix media are used to describe the occurrence and migration of adsorbed gas, free gas and water, while microcrack media and hydraulic crack media are used to describe the rapid flow and transport in the crack system.

[0052] S1.2, two isotopic components in methane Establish the component conservation equations separately, and characterize them uniformly as follows: ; in, as medium Middle Isotope Components Storage capacity per unit volume; Isotopic components In the medium Total flux in; For inter-media exchange items; ∇⋅ is the well source term; ∇⋅ is the divergence operator.

[0053] S1.3, Storage Items in the Matrix Medium Written as: ; Reserves in microfractured and hydraulically fractured media , Write them as follows: ; ; in, , , They are matrix media Microcracked media and hydraulic fracture media porosity; , , They are matrix media Microcracked media and hydraulic fracture media Gas phase saturation; This refers to the gas phase density. isotopic components in the gas phase mole fraction; Bulk density; This refers to the reserves of adsorbed isotopes.

[0054] Specifically, step S1 also includes the following steps: S1.4, Establish the initial discrete mesh. The matrix region is discretized using a conventional volume mesh, the hydraulic fractures are discretized using EDFM embedded discretization, and the microfractures are characterized using a dual-medium continuous model. For the exchange zone between the matrix and the fractures, at least one MINC transient nested sublayer is pre-placed near the fractures to describe the early unsteady recharge process.

[0055] This part forms the physical basis for the subsequent adaptive solution. In particular... 12 CH4 and 13 CH4 must be solved independently by conservation, and cannot be solved directly by... Establish the equation. Because... It is not a conserved quantity itself and cannot be directly used for mass balance; it can only be obtained by conversion from the calculation results of 12CH4 and 13CH4. The conversion relationship is as follows: , ; in, The standard ratio, and They are respectively 13 CH4 and 12 molar amount of CH4 For the sample or control volume13 CH4 and 12 The molar ratio of CH4; The ratio of 13C to 12C isotopic abundance in the standard sample; The carbon isotope offset relative to the standard ratio is expressed in per mille (‰).

[0056] Specifically, step S2 includes the following steps: S2.1, Constructing the pressure gradient index For any control body ,definition: ; in, To control the length of the volume feature; To control the average pressure of the body; To prevent scalar pressure from an excessively small denominator, For control body The internal pressure gradient magnitude represents the intensity of pressure change per unit distance near the control volume. The control volume is a general concept in the finite volume method, referring to a discrete volume element after the computational domain has been divided. It can belong to the matrix, microfractures, or hydraulic fractures.

[0057] This index characterizes the strength of changes in the local pressure field. Numerically, if... A large value indicates that there is a significant pressure drop, near-fracture transient flow, or pressure disturbance at the fracture tip in the region. Conventional coarse meshes are prone to causing distortion in pressure wave propagation.

[0058] S2.2, Definition Gradient index : ; in, The methane carbon isotope value of the current control body; It is an isotopic characteristic scale; It is a small positive number.

[0059] This index is specifically used to identify isotopic leading regions. This is because in isotopic fractionation problems, when the pressure field has already leveled off, δ0... 13 C-forward may continue to advance; if the grid is adjusted only based on the pressure gradient, the isotope fractionation interface will be "flattened", resulting in a distortion of the fractionation intensity and stage boundary judgment.

[0060] S2.3, for control volumes that are not adjacent to fracture elements, defines the isotopic exchange residual index between the matrix and the fracture. for: ; in, For time step Time control body Internal isotope components reserves, For the current time step, For control body A boundary surface, For control body The set of all boundary surfaces, To control the component flux at each interface of the volume; It serves as a source and sink for exchange between the matrix and the cracks.

[0061] Define the isotopic exchange residual index between the matrix and the fracture. : ; in, To prevent tiny positive numbers with a denominator of zero, , Control body middle 12 CH4 and 13 Exchange residuals between the matrix and cracks of CH4.

[0062] The physical meaning of this index is: it directly measures the transport of the matrix into the fracture under the current discrete method. 12 CH4 and 13 Is the conservation error of CH4 too large? If the value is too large, it indicates that it is not simply a matter of mesh size, but rather that the local exchange process itself has been distorted by discretization error. In this case, the discretization of the MINC sublayer or the crack neighborhood must be corrected first.

[0063] S2.4, after making the three types of indicators dimensionless, define the joint error index. : ; in, , , , These are the weights for pressure, isotopes, and exchange residuals, respectively. They are dimensionless versions , , Prioritize improving the area adjacent to the crack. In the far field region, the speed can be appropriately increased. ; Improve in the fractionation front region .

[0064] Traditional adaptive methods mainly refine based on pressure, saturation, or velocity gradients, while this invention refines based on δ. 13The introduction of C gradient and isotope exchange residuals into the error control index enables local refinement to serve not only the flow field, but also the analysis of the fractionation front and the matrix replenishment mechanism.

[0065] Specifically, step S3 includes the following steps: S3.1, Adaptive Update of MINC Nesting Layers and Transition Points: For matrix mesh blocks near cracks, a fixed number of MINC layers is not used; instead, the number of transient nesting layers is determined by the joint error index. and transition point .

[0066] S3.2, Embedded Discrete Crack Model EDFM Crack Neighborhood Local Refinement and Non-Adjacent Connection Transmission Coefficient Update: For satisfying The control body performs local refinement; for several consecutive macro time steps, it satisfies... The control body performs local coarsening, where , These are the refinement threshold and the coarsening threshold, respectively.

[0067] S3.3, Multi-rate Time Advancement: A uniform time step is not used for different regions; instead, the local time step is determined based on local spatial scale, convection velocity, pressure diffusion rate, and isotope diffusion rate. : ; in, These are the time step control coefficients for convection, pressure diffusion, and isotope diffusion, respectively. Local velocity; To control the length of the volume feature; To prevent a small positive number with a zero denominator due to a zero velocity; For control body Pressure diffusion rate, For control body The effective diffusion coefficient of isotopes.

[0068] S3.4 divides the computational domain into coarse time step region, medium time step region, and fine time step region, and performs sub-loop iteration within a macro time step in the fine time step region; for the fine time step region, based on the local time step size... Progressing step by step, when the residuals of the pressure equation within a certain sub-time step are... 12 CH4 13 When all CH4 mass-conserving residuals meet the preset nonlinear iterative convergence condition, the calculation result of that sub-time step is accepted and the process proceeds to the next sub-time step; if the convergence condition is not met, the local time step size is reduced and the sub-time step is recalculated; when the cumulative advancement time of the local sub-cycle reaches the macro time step... At this point, the sub-cycle in this region stops, and the dual isotope conservation correction begins at the end of the macro time step.

[0069] Smaller time steps are used in the crack neighborhood and high isotope gradient region, while larger time steps are maintained in the far field region.

[0070] The physical significance of this step lies in the fact that it is not simply about "refining the mesh near the crack," but rather about simultaneously determining, through three types of error indices: whether a MINC transient sublayer needs to be added to the matrix side; whether the transition point should be moved outward or inward; whether the background mesh in the crack neighborhood needs to be refined; and whether the local time step needs to be reduced. Therefore, this invention corrects the three core processes of pressure propagation, fractionation front, and matrix replenishment, rather than simply correcting the pressure field.

[0071] Specifically, step S3.1 includes the following steps: S3.1.1, when the control body that is not adjacent to the crack element satisfies: or This indicates pressure disturbances, isotope fractionation leading edges, or other factors within the control body. If the value is too large, it triggers a re-partitioning of the MINC sublayer and sets the transition point accordingly. Moving inwards into the matrix; among which, To exchange residual refinement thresholds, This is an index of isotopic exchange residuals between the matrix and the cracks.

[0072] When the adjacent control bodies of the crack are continuous Each time step satisfies: and This indicates that the changes in pressure and isotopes within the control body tend to level off or become more gradual. The smaller value triggers the merging of the outer MINC sublayers and regresses the transition point towards the crack direction; among which, To exchange residual coarsening thresholds.

[0073] In a preferred embodiment, the following can be taken: =1.0, =0.35, =5.0×10 −3 =1.0×10 −3 , =3.

[0074] S3.1.2, Define the pressure diffusion length and isotope diffusion length : , ; in, For pressure diffusion rate, This refers to the permeability of the medium, which can be the permeability of the matrix or the equivalent permeability. These are porosity, overall compressibility, and fluid viscosity, respectively. is the effective diffusion coefficient of the isotope.

[0075] The thickness of the innermost nested layer is taken as : ; in, This is the proportionality coefficient. The function is to minimize the value; this formula means that the first MINC sublayer must be able to simultaneously distinguish the penetration scale of pressure perturbation and isotopic perturbation within one time step.

[0076] S3.1.3, each sub-layer expands outward in a geometric progression: ; in, The layer thickness expansion factor is... For the first MINC sublayer thickness, This represents the total number of transient MINC sublayers retained in the current matrix block.

[0077] For transition points The following criteria shall be used to determine this: and ; in, and These are small positive numbers used to prevent the denominator from being zero. and These are the allowable error thresholds for the pressure balance criterion and the isotope balance criterion, respectively. The MINC sublayer number corresponding to the candidate transition point; and The first The mean pressure and mean carbon isotope value of the MINC sublayer. and These are the volume-weighted average pressure and volume-weighted average carbon isotope value within the candidate equivalent equilibrium region, respectively. This is the function for finding the maximum value.

[0078] That is, when the region outside a certain layer simultaneously satisfies "approximate isobaric pressure" and "approximate isotopic equivalence", the region outside that layer is merged into an equivalent equilibrium block, and the multi-layer transient nesting is no longer maintained.

[0079] S3.1.4, MINC interlayer component exchange conductivity Calculated using equivalent resistance: ; in, Equivalent permeability of components; For component fluidity; For the exchange cross-sectional area.

[0080] It should be noted that, Used to describe the internal exchange conduction capability between adjacent MINC sublayers; and They are used to describe the non-adjacent connection exchange capability between the matrix and fracture elements, and their objects of action are different.

[0081] With this treatment, the total drag remains unchanged when the sublayer decreases, and the near-fracture transient resolution improves when the sublayer increases.

[0082] Specifically, after performing local refinement or local coarsening, the non-adjacent connectivity (NNC) conductivity between the matrix and fractures is recalculated; isotopic composition Exchange flux Written as: ; in: , ; in, Isotopic components The matrix-to-microcrack convective exchange conductivity coefficient is used to characterize the ability of this component to migrate from the matrix to the microcracks under pressure difference. Isotopic components The diffusion exchange conductivity coefficient between the matrix and the cracks is used to characterize the ability of the component to diffuse and migrate between the matrix and microcracks under concentration gradient drive. This represents the exchange area between the matrix and the crack. Equivalent normal distance; Normal penetration rate; These represent the component concentrations of the matrix medium and the microcrack medium, respectively. The effective diffusion coefficient; To control the average pressure of the body, The average pressure of the microcrack element connected to the control volume. Isotopic components The equivalent mobility can be expressed as the product of the gas phase mobility and the mole fraction of the component; For matrix porosity; This represents the gas phase saturation of the matrix.

[0083] After local refinement or local coarsening, and Update according to the new geometric relationship, thereby correcting both convection and diffusion exchange simultaneously.

[0084] Specifically, the time progression termination loop condition in step S3.4 is: ; That is, the cumulative advancement time of the local sub-loop reaches the macro time step. When this happens, the sub-loop in that region stops.

[0085] Nonlinear iterative convergence condition: ; ; in, For the residuals of the pressure equation, , They are respectively 12 CH4 and 13 The mass conservation residual of CH4 and The first After the first nonlinear iteration, the computational region 12 CH4 and 13 Storage vector of CH4.

[0086] When the cumulative time of the local sub-loop reaches the macro time step, and both the pressure residual and the dual isotope mass residual meet the convergence threshold, the sub-loop within that macro time step ends; otherwise, the local time step is reduced and the calculation is repeated.

[0087] Specifically, step S4 includes the following steps: S4.1, Define the total residual of components within the macro time step: , ; in, , Each within a macro time step 12 CH4 13 The overall mass conservation residual of CH4 , The calculation region at the start of the macro time step 12 CH4 13 Total reserves of CH4 , Calculate the region at the end of the macro time step respectively 12 CH4 13 Total reserves of CH4 These are the macro time steps within which 12 CH4 13 The net source and sink of CH4 are defined as positive when entering the calculation region and negative when leaving the calculation region.

[0088] S4.2, in the high-error local region Internally, solve the minimum correction problem: ; in, As a weight, it is preferentially allocated to areas with high error rates or high reserves.

[0089] Satisfy constraints: , ; in, , This is the local element quality correction amount.

[0090] S4.3, after correction, recalculate the isotope ratios and δ. 13 C, output pressure field, gas production rate 12 CH4 and 13 CH4 component distribution, exchange flux between matrix and fracture and δ 13 C is distributed in space and time.

[0091] To verify the beneficial effects of the present invention, the following experiments were conducted: The study focuses on fracturing-enhanced coal-gas reservoirs. These reservoirs contain matrix porosity, microfractures, and artificial hydraulic fractures. Methane exists in both free and adsorbed states in the matrix, and primarily migrates in the free state within the microfractures and hydraulic fractures. During production, 12 CH4 and 13 CH4 exhibits varying migration capabilities during adsorption, desorption, diffusion, and exchange between the matrix and cracks, resulting in different production gas δ 13 C changes over time.

[0092] Conventional MINC–EDFM methods typically adjust the mesh based on the pressure field or saturation field, making it difficult to accurately describe δ. 13 The C-front and the isotopic exchange process between the matrix and the fracture. This embodiment proposes a process involving pressure gradient, δ 13 An adaptive correction method driven by both the C gradient and the isotopic exchange residual between the matrix and the crack provides a clear physical basis for mesh refinement, MINC layer number adjustment, and time step control.

[0093] The technical solution of this embodiment includes the following steps: S1: Establishing the initial model of dual isotope MINC–EDFM The reservoir is divided into matrix media. Microcracked media and hydraulic fracture media The unsteady-state exchange between the matrix and microcracks is described using the multiple interaction continuum model (MINC), while the embedded fracture exchange between microcracks and hydraulic fractures is described using the embedded discrete fracture model (EDFM). The model dimensions are 600m × 300m × 30m, the basic mesh is 60 × 30 × 3, and the basic mesh size is 10m × 10m × 10m.

[0094] A horizontal production well is located at the center of the model, with a horizontal section length of 400m. Five hydraulic fractures are arranged along the horizontal well. The half-length of each hydraulic fracture is 80m, the height is 30m, the fracture spacing is 80m, and the initial fracture aperture is 3mm.

[0095] right 12 CH4 and 13 Establish the mass conservation equations for CH4 respectively: ; in, =12, 13, respectively represent 12 CH4 and 13 CH 4。

[0096] The component reserves in the matrix medium are written as: ...

[0097] Reserves in microfractured and hydraulically fractured media , Write them as follows: ; .

[0098] The isotopic values ​​of the produced gas are from 12 CH4 and 13 Calculation of the molar amount of CH4 produced: , ;

[0099] Reservoir and fluid parameters are shown in Table 1: Table 1 Reservoir and fluid parameters

[0100] in, 13 The diffusion coefficient of CH4 is slightly less than 12 CH4 is used to characterize the physical properties of heavy isotopes with weak migration ability.

[0101] S2: Constructing a joint error index driven by isotopic errors After each macro time step, for each control volume Calculate the three error indices.

[0102] (1) Pressure gradient index ;

[0103] This indicator is used to identify areas of intense pressure propagation, such as fracture tips, near-wellbore areas, and pressure drop fronts.

[0104] (2) δ 13 C-gradient index ; in, =10‰, =10 −6 .

[0105] This index is used to identify isotopic fractionation fronts. A larger index indicates a higher δf between adjacent grids. 13 The C difference is significant, and coarse meshes weaken the intensity of isotope fractionation.

[0106] (3) Isotope exchange residual index between matrix and fracture For the control volumes adjacent to the crack, calculate separately 12 CH4 and 13 Local exchange residuals of CH4: .

[0107] Further definition: ; in, =10 −12 .

[0108] This index is used to determine whether the current discretization method accurately describes the isotopic supply from the matrix to the fracture. If A large value indicates a mismatch between the exchange flux between the matrix and the crack and the changes in the quality of local components, requiring adjustments to the number of MINC layers or the crack neighborhood mesh.

[0109] The three types of indicators are combined into a joint error index: .

[0110] In this embodiment, δ 13 The gradient weight C is set to 0.40 because this method mainly serves isotope fractionation simulation, and it is necessary to prioritize the accuracy of isotope front analysis.

[0111] S3: Adaptive correction of the MINC–EDFM discrete structure based on the joint error index.

[0112] (1) Adaptive adjustment of the number of MINC sub-layers and transition points

[0113] Set the number of sublayers N in the MINC sublayer to be near the crack. MINC In this embodiment: , That is, the minimum number of MINC sub-layers N MINC The maximum number of MINC sub-layers is 2, N. MINC It is 8.

[0114] When a certain matrix block satisfies: >1.0 or >5.0×10 −3 Then, a MINC sublayer is added inside the matrix block.

[0115] When a matrix block satisfies the following for three consecutive macro time steps: <0.35 and <1.0×10 −3 Then the outermost MINC sublayer is merged.

[0116] The thickness of the innermost layer of MINC is determined by both the pressure diffusion length and the isotopic diffusion length: ;

[0117] in: ; ; Among them, D 12 for 12 The effective diffusion coefficient of CH4 in the corresponding medium; D 13 for 13 The effective diffusion coefficient of CH4 in the corresponding medium.

[0118] In this embodiment, the thickness of the MINC sublayer expands outward in a geometric progression: .

[0119] transition point Determined by the combined effect of approximate isobaric pressure and approximate isotopic equivalence: ; and .

[0120] When the two conditions mentioned above are met simultaneously on the outer side of a certain layer, the outer side of that layer is merged into an equivalent equilibrium region, and transient sublayers are no longer retained.

[0121] Near-fracture zone pressure and δ 13The C-means change drastically, requiring the retention of multiple MINC layers to describe the unsteady supply; the pressure and isotope changes in regions far from the fracture tend to be gradual, and can be merged into an equivalent equilibrium region to reduce the computational load.

[0122] (2) Local refinement and exchange coefficient update of crack neighborhood in EDFM

[0123] Local refinement is performed on the background mesh that meets the following conditions: >1.0; Alternatively, the grid may be located within 20m of the tip of a hydraulic fracture.

[0124] Local refinement employs a bisection method, where the original mesh is divided into 2×2 sub-mesh sections in the planar direction. After refinement, the non-adjacent connectivity conduction coefficient between the matrix and the cracks is recalculated.

[0125] For component a, the exchange flux between the matrix and the fracture is written as: ;

[0126] The convective exchange conduction coefficient is: ;

[0127] The diffusion-exchange conductivity is: ; After local refinement, and Update according to the new geometric relationship, thereby correcting both convection and diffusion exchange simultaneously.

[0128] (3) Multi-rate time propagation Based on the local joint error index, the computational domain is divided into three categories: Table 2. Regional Type Table

[0129] Within a macro time step, the high error region performs 8 sub-loop calculations, the medium error region performs 4 sub-loop calculations, and the low error region performs only 1 calculation.

[0130] The area near the crack and the isotope front region changes rapidly, so small time steps must be used to ensure the accuracy of isotope exchange; the far-field region changes slowly, so large time steps can be used to reduce computational costs.

[0131] S4: Perform dual isotope conservation correction at the end of the macro time step.

[0132] Because local refinement, coarsening, and multi-rate advancement can lead to inconsistent computational step sizes in different regions, this embodiment performs a step-by-step adjustment at the end of each macro time step. 12 CH4 and 13 CH4 was subjected to conservation correction.

[0133] Define the total residual of the components within the macro time step: ; in, =12,13.

[0134] When the following conditions are met: .

[0135] Initiate local conservation correction.

[0136] The correction region is selected to satisfy η c >0.35 control body set The correction amount is determined by the following minimum perturbation problem: ;

[0137] The constraints are: , ; Among them, weight Take as: That is, the larger the reserve, the greater the amount of minute correction required. and The control body at the end of the macro time step Inside 12 CH4 and 13 CH4 reserves.

[0138] Recalculate after correction: , ; Thus guarantee 12 CH4 13 CH4 and δ 13 The C output is based on the principle of mass conservation.

[0139] S5: To verify the effectiveness of the present invention, this embodiment sets up three sets of comparative calculations: Option A: Fixed MINC–EDFM method: MINC is fixed at 3 layers, the background mesh is not locally refined, and the time step is uniformly 10d.

[0140] Option B: Pressure-driven adaptive method: Local refinement and time step adjustment are performed only based on the pressure gradient.

[0141] Scheme C: The method of the present invention: joint adaptive correction based on pressure gradient, δ¹³C gradient and matrix-crack isotope exchange residual.

[0142] A high-precision reference model was also set up: MINC was fixed at 10 layers, global and local densification was applied to the crack neighborhood, and a unified time step of 1 day was used. Using this high-precision model as an error reference, the calculation results in a set of exemplary numerical examples are shown in Table 3.

[0143] Table 3 Comparison of Calculation Results by Different Methods

[0144] From Table 3 and Figure 2 The calculation results show that although the fixed MINC–EDFM method is faster, δ 13 The large C error indicates that coarse grids and a fixed number of MINC layers are insufficient to accurately describe the isotopic fractionation front. Pressure-driven adaptive methods can improve pressure field and gas production errors, but they fail to identify δ... 13 C gradient and isotope exchange residuals, the resulting gas δ 13 The error C still reaches 0.82‰.

[0145] The method of this invention, when the average number of active grids is lower than that of the pressure-driven adaptive method, will produce gas δ 13 The root mean square error of C decreased to 0.29‰, and the maximum isotope exchange residual decreased to 4.8 × 10⁻⁶. −4 Meanwhile, the calculation speed is 5.8 times faster than that of the high-precision model. This demonstrates that the present invention can significantly reduce computational costs while maintaining the accuracy of isotope fractionation simulation.

[0146] Existing MINC–EDFM numerical simulation methods mostly use pressure or saturation fields as the basis for grid correction, which can improve the accuracy of flow calculations, but they have difficulty accurately identifying the methane carbon isotope fractionation front and determining the matrix-fracture relationship. 12 CH4 13 Does the CH4 exchange process have local conservation errors that could easily cause δ 13 C. Prediction bias. This invention utilizes pressure gradient, δ... 13 The C gradient and matrix-crack isotope exchange residuals together construct a joint error index, enabling mesh correction to serve not only pressure propagation, but also the pressure field, isotope field, and inter-medium exchange process.

[0147] Compared with the fixed MINC–EDFM method, the present invention can significantly reduce the produced gas δ 13 C-prediction error is reduced, improving the descriptive ability of matrix replenishment and isotope fractionation processes; compared with the simple pressure-driven adaptive method, this invention can avoid the problem of "accurate pressure calculation but distorted isotope prediction". At the end of the macro time step, respectively... 12 CH4 and 13 CH4 is used for conservation correction to ensure the final δ 13The C result originates from actual component mass variations, rather than numerical artifacts caused by local mesh variations or multi-rate propagation. Therefore, this invention has the combined effects of improving the accuracy of isotope fractionation simulations, reducing invalid encryption, increasing computational efficiency, and enhancing the physical reliability of the results.

[0148] First, the pressure gradient and δ in the joint error index 13 The C-gradient and isotope exchange residuals do not necessarily require fixed weights and can be adaptively adjusted according to the simulation stage. For example, in the early stages of production, when pressure propagation is intense, the pressure gradient weight can be increased; in the middle and later stages, when adsorbed gas desorption and diffusion are enhanced, the δ-gradient weight can be increased. 13 The C-gradient and isotope-exchange residual weights are also relevant. Secondly, the MINC sublayer adjustment method can be replaced by a continuous equivalent transfer function method, i.e., without explicitly adding nested layers, but by dynamically modifying the matrix-fracture conduction coefficient to describe the unsteady supply. Thirdly, EDFM local refinement can be achieved using quadtree meshes, unstructured mesh local reconstruction, or multi-mesh local correction methods, and is not limited to regular mesh bisection refinement. Multi-rate time advancement can also be replaced by local time steps, adaptive implicit time steps, or a hybrid method of local explicit-global implicit time steps.

[0149] This invention can also be extended to other numerical simulation fields involving "multi-field coupling, cross-medium exchange, component fractionation, or tracer response." For example, in carbon dioxide geological storage, δ 13 The C gradient is replaced with a CO2 concentration gradient or a carbon-oxygen isotope gradient, and the matrix-fracture isotope exchange residual is replaced with a CO2 dissolution-diffusion-mineralization reaction residual, which is used to adaptively characterize CO2 plume migration and reaction fronts. In shale oil or tight oil development, isotopic indices can be replaced with light and heavy component concentration gradients to describe the selective transport of multi-component fluids in nanopores.

[0150] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A MINC-EDFM solution method considering isotopic fractionation, characterized in that, Specifically, the steps include the following: S1, Establish the three-medium dual-isotope coupling model and the initial discrete grid; S2, constructing a joint error index driven by isotope error; S3, based on the joint error index, adaptively corrects the transition point of the multi-interaction continuous medium model MINC, the crack neighborhood grid and local time step of the embedded discrete crack model EDFM. S4 performs isotope conservation correction at the end of the macro time step and outputs the final result; Step S2 specifically includes the following steps: S2.1, Constructing the pressure gradient index For any control body ,definition: ； in, To control the length of the volume feature; To control the average pressure of the body; To prevent scalar pressure from an excessively small denominator, For control body Internal pressure gradient magnitude; S2.2, Definition Gradient index : ； in, The methane carbon isotope value of the current control body; It is an isotopic characteristic scale; It is a tiny positive number; S2.3, for control volumes that are not adjacent to fracture elements, defines the isotopic exchange residual index between the matrix and the fracture. for: ； in, For time steps Time control body Internal isotope components reserves, For the current time step, For control body A boundary surface, For control body The set of all boundary surfaces, To control the component flux at each interface of the volume; It serves as a source and sink for exchange between the matrix and the cracks; Define the isotopic exchange residual index between the matrix and the fracture. : ； in, To prevent tiny positive numbers with a denominator of zero, , Control body middle 12 CH4 and 13 Exchange residuals between the matrix and cracks of CH4; S2.4, after making the three types of indicators dimensionless, define the joint error index. : ； in, , , , These are the weights for pressure, isotopes, and exchange residuals, respectively. They are dimensionless versions , , .

2. The MINC-EDFM solution method considering isotopic fractionation according to claim 1, characterized in that, Step S1 specifically includes the following steps: S1.1, the reservoir is divided into matrix media. Microcracked media and hydraulic fracture media ; S1.2, two isotopic components in methane Establish the component conservation equations separately, and characterize them uniformly as follows: ； in, as medium Middle Isotope Components Storage capacity per unit volume; Isotopic components In the medium Total flux in; For inter-media exchange items; For well source and sink terms; ∇⋅: ​​for divergence operator; S1.3, Storage Items in the Matrix Medium Written as: ； Reserves in microfractured and hydraulically fractured media , They are written as: ； ； in, , , They are matrix media Microcracked media and hydraulic fracture media porosity; , , They are matrix media Microcracked media and hydraulic fracture media Gas phase saturation; This refers to the gas phase density. isotopic components in the gas phase mole fraction; Bulk density; This refers to the reserves of adsorbed isotopes.

3. The MINC-EDFM solution method considering isotopic fractionation according to claim 2, characterized in that, Step S1 also includes the following steps: S1.4, the matrix region is discretized using conventional volumetric meshes, the hydraulic fractures are discretized using EDFM embedded discretization, and the microfractures are characterized using dual-medium continuous characterization; for the exchange zone between the matrix and the fractures, at least one MINC transient nested sublayer is pre-placed near the fractures to describe the early unsteady recharge process.

4. The MINC-EDFM solution method considering isotopic fractionation according to claim 1, characterized in that, Step S3 specifically includes the following steps: S3.1, For matrix mesh blocks near cracks, the number of transient nested layers is determined by the joint error index. and transition point ; S3.2, for satisfying The control body performs local refinement; For a consecutive number of macro time steps, satisfy The control body performs local coarsening, where , These are the refinement threshold and the coarsening threshold, respectively. S3.3, For different regions, the local time step is determined based on the local spatial scale, convection velocity, pressure diffusion rate, and isotope diffusion rate. : ； in, These are the time step control coefficients for convection, pressure diffusion, and isotope diffusion, respectively. Local velocity; To control the length of the volume feature; To prevent a small positive number with a zero denominator due to a zero velocity; For control body Pressure diffusion rate, For control body The effective diffusion coefficient of isotopes; S3.4 divides the computational domain into coarse time step region, medium time step region, and fine time step region, and performs sub-loop iteration within a macro time step in the fine time step region; for the fine time step region, based on the local time step size... Progressing step by step, when the residuals of the pressure equation within a certain sub-time step are... 12 CH4 13 When all CH4 mass conservation residuals meet the preset nonlinear iteration convergence conditions, the calculation results of this sub-time step are accepted and the next sub-time step is entered; if the convergence conditions are not met, the local time step size is reduced and the sub-time step is recalculated. When the cumulative execution time of a local sub-loop reaches the macro time step At this point, the sub-cycle in this region stops, and the dual isotope conservation correction begins at the end of the macro time step.

5. The MINC-EDFM solution method considering isotopic fractionation according to claim 4, characterized in that, Step S3.1 specifically includes the following steps: S3.1.1, when the control body that is not adjacent to the crack element satisfies: or At that time, the MINC sublayer is repartitioned, and the transition point is... Moving inwards into the matrix; among which, To exchange residual refinement thresholds, This is an index of isotopic exchange residuals between the matrix and the fracture. When the adjacent control bodies of the crack are continuous Each time step satisfies: and At that time, the outer MINC sublayer is merged, and the transition point is pushed back towards the crack direction; among them, To exchange residual coarsening thresholds; S3.1.2, Define the pressure diffusion length and isotope diffusion length : ， ； in, For pressure diffusion rate, For the medium permeability, These are porosity, overall compressibility, and fluid viscosity, respectively. The effective diffusion coefficient of the isotope; The thickness of the innermost nested layer is taken as : ； in, This is the proportionality coefficient. The function is for finding the minimum value; S3.1.3, each sub-layer expands outward in a geometric progression: ； in, The layer thickness expansion factor is... For the first MINC sublayer thickness, This represents the total number of transient MINC sublayers retained in the current matrix block; For transition points When the outer perimeter of a certain layer simultaneously meets the following two conditions, the outer perimeter of that layer is merged into an equivalent equilibrium region, and transient sublayers are no longer retained. The conditions are: and ; in, and These are small positive numbers used to prevent the denominator from being zero. and These are the allowable error thresholds for the pressure balance criterion and the isotope balance criterion, respectively. The MINC sublayer number corresponding to the candidate transition point; and The first The mean pressure and mean carbon isotope value of the MINC sublayer. and These are the volume-weighted average pressure and volume-weighted average carbon isotope value within the candidate equivalent equilibrium region, respectively. This is a function to find the maximum value. S3.1.4, MINC interlayer component exchange conductivity Calculated using equivalent resistance: ； in, Equivalent permeability of components; For component fluidity; For the exchange of cross-sectional areas.

6. The MINC-EDFM solution method considering isotopic fractionation according to claim 1, characterized in that, After performing local refinement or local coarsening, the non-adjacent connectivity conductivity between the matrix and the fracture is recalculated; isotopic composition. Exchange flux Written as: ； in: , ； in, Isotopic components The matrix-to-microcrack convection exchange conductivity, Isotopic components The diffusion exchange conductivity between the matrix and the crack, This represents the exchange area between the matrix and the crack. Equivalent normal distance; Normal penetration rate; These represent the component concentrations of the matrix medium and the microcrack medium, respectively. The effective diffusion coefficient; To control the average pressure of the body, The average pressure of the microcrack element connected to the control volume. Isotopic components The equivalent mobility; For matrix porosity; This represents the gas phase saturation of the matrix.

7. The MINC-EDFM solution method considering isotopic fractionation according to claim 4, characterized in that, The time progression termination loop condition in step S3.4 is: ； That is, the cumulative advancement time of the local sub-loop reaches the macro time step. When this happens, the sub-loop in that region stops; Nonlinear iterative convergence condition: ; ; in, For the residuals of the pressure equation, , They are respectively 12 CH4 and 13 The mass conservation residual of CH4 and The first After the first nonlinear iteration, the computational region 12 CH4 and 13 Storage vector of CH4.

8. The MINC-EDFM solution method considering isotopic fractionation according to claim 1, characterized in that, Step S4 specifically includes the following steps: S4.1, Define the total residual of components within the macro time step: ， ； in, , Each within a macro time step 12 CH4 13 The overall mass conservation residual of CH4 , The calculation region at the start of the macro time step 12 CH4 13 Total reserves of CH4 , Calculate the region at the end of the macro time step respectively 12 CH4 13 Total reserves of CH4 These are the macro time steps within which 12 CH4 13 The net source and sink of CH4 are defined as positive when entering the calculation region and negative when leaving the calculation region. S4.2, in the high-error local region Internally, solve the minimum correction problem: ； in, As weight; Satisfy constraints: ， ； in, , This is the local element quality correction amount; S4.3, after correction, recalculate the isotope ratios and δ. 13 C, output pressure field, gas production rate 12 CH4 and 13 CH4 component distribution, exchange flux between matrix and fracture and δ 13 C is distributed in space and time.

Images