Fractured medium equivalent permeability full-tensor calculation method based on EDFM

Abstract

The invention discloses a fractured medium equivalent permeability full-tensor calculation method based on EDFM, and the method comprises the following steps: S1, generating a random fracture, and constructing a fractured medium physical model; s2, dividing a matrix grid and a crack grid based on EDFM, and endowing basic parameters to the grids; s3, constructing a fractured medium seepage model based on EDFM; s4, solving pressure distribution in the whole computational domain by adopting a three-class boundary condition loading mode; s5, calculating pressure gradient and flow velocity distribution, and calculating average pressure gradient and average flow velocity in the whole computational domain; and S6, calculating an equivalent permeability full tensor matrix, and carrying out symmetric processing on the equivalent permeability full tensor matrix to obtain a calculation result. According to the method, the problems of low calculation efficiency and high calculation cost caused by complex fracture grid form or too many fractures in the fractured medium can be solved, the calculated equivalent permeability full tensor can well represent heterogeneity and anisotropy of the fractured medium, and technical support is provided for development of the fractured medium.

Description

Technical Field

[0001] The invention relates to the technical field of fractured medium development, and in particular to an EDFM-based fractured medium equivalent permeability full tensor calculation method. Background Art

[0002] Fractured media are difficult to simulate because they contain bedrock pores and rich fracture structures, have many types of media, large scale differences, and complex flow laws. The commonly used mathematical model is the equivalent continuum medium model. The equivalent continuum medium model solves the equivalent permeability tensor of the matrix grid block containing fractures, and uses the equivalent continuum medium body represented by the equivalent permeability tensor to replace the original fracture medium body for simulation.

[0003] Numerical simulation of oil and gas reservoirs is a common method for evaluating reservoir dynamics. It can be used to solve seepage or engineering problems in complex situations. However, its simulation scale is large, and the grid scale is generally tens to hundreds of meters. A grid may contain many natural fractures. It is obviously difficult to display a large number of natural fractures in the entire oil and gas reservoir, and the calculation cost is too high. Therefore, for units with developed fractures, certain methods are used to calculate the comprehensive equivalent permeability of this unit including the matrix and fracture effects, which can reduce the calculation cost of numerical simulation of oil and gas reservoirs.

[0004] When calculating the comprehensive equivalent permeability, it is necessary to characterize the cracks within the unit. Common crack characterization methods include discrete crack model, boundary element method, dual medium method, etc. Among them, the discrete crack model method is the most popular. When using the DFM method to calculate the full tensor of equivalent permeability, the model directly calculates the conduction between the matrix and the cracks, so fine grid division is required. The DFM method divides the matrix and cracks by using irregular grids. Especially when dealing with cracks, it is necessary to treat the cracks as independent grid units, resulting in high computational cost when calculating the full tensor of equivalent permeability. Summary of the invention

[0005] In view of the above problems, the present invention aims to provide a full tensor calculation method for equivalent permeability of fractured media based on EDFM.

[0006] The technical solution of the present invention is as follows:

[0007] A method for calculating equivalent permeability full tensor of fractured media based on EDFM, comprising the following steps:

[0008] S1: Generate random cracks and build a physical model of fractured media;

[0009] S2: dividing the matrix grid and the fracture grid on the fracture medium physical model based on EDFM, and assigning basic grid parameters;

[0010] S3: Construct fracture medium seepage model based on EDFM;

[0011] S4: adopting three types of boundary condition loading modes, and carrying out steady-state seepage numerical simulations according to the fracture medium seepage model to solve the pressure distribution in the entire calculation domain;

[0012] S5: According to the pressure distribution in the entire calculation domain, the pressure gradient and the flow velocity distribution are calculated respectively, and then the average pressure gradient and the average flow velocity in the entire calculation domain are calculated respectively by volume weighted average;

[0013] S6: Based on the generalized Darcy's law, an equivalent permeability full tensor matrix is ​​calculated and obtained, and the equivalent permeability full tensor matrix is ​​symmetrically processed to obtain a calculation result of the equivalent permeability full tensor of the fractured medium based on EDFM.

[0014] Preferably, in step S1, when generating random fractures, they are randomly generated based on geological statistical analysis of fracture attribute parameters and representative statistical characteristics of fracture attribute parameters.

[0015] Preferably, in step S3, the fractured medium seepage model includes a flow control equation and a crossflow term calculation equation;

[0016] The flow control equation is:

[0017]

[0018] Where: φ m ,φ f are matrix and fracture porosity, dimensionless; ρ is density, kg / m 3 ; t is time, s; represents the divergence operator, dimensionless; u m 、u f are the seepage velocities in the matrix and fractures, m / s; α f is the aperture of the crack, m; Ψ f，m , m，f are the cross-flow terms between fracture units and matrix units, and the cross-flow terms between matrix units and fracture units, respectively. (m 3 / s) / m 3 ; f1，f2 is the cross-flow term between fracture units f1 and f2, (m 3 / s) / m 2 ;

[0019] The cross-flow term calculation equation is:

[0020]

[0021] Where: L m，fis the length of the fracture element f within the matrix element m, m; <d> m，f is the average distance between matrix unit m and fracture unit f, m; K m is the matrix permeability, m 2 ; μ is the viscosity of the fluid in the grid, Pa·s; P m , P f are matrix and fracture pressure, Pa; α f1 , α f2 are the crossflow coefficients of fracture unit f1 and fracture unit f2, dimensionless; P f1 , P f2 are the pressures of fracture unit f1 and fracture unit f2, Pa, respectively;

[0022] The crossflow coefficients of the crack unit f1 and the crack unit f2 are calculated by the following formula:

[0023]

[0024] In the formula: FA f is the aperture of the fracture unit f, m; λ t，f is the total fluidity of the fracture unit f, dimensionless; dx f is the spacing of crack units, m.

[0025] Preferably, in step S4, the three types of boundary condition loading modes refer to setting boundary conditions on two surfaces corresponding to the three directions of x, y and z respectively;

[0026] In the x direction, the two corresponding surfaces are yz and y′z′; in the y direction, the two corresponding surfaces are xz and x′z′; in the z direction, the two corresponding surfaces are xy and x′y′;

[0027] The boundary condition is that along one direction, a higher pressure one is applied on one of the corresponding surfaces, a lower pressure two is applied on the other corresponding surface, the matrix grids of the remaining four surfaces serve as closed boundaries, and the pressure two is applied on the remaining four fracture unit surfaces.

[0028] Preferably, in step S5, the flow velocity distribution is calculated by the following formula:

[0029]

[0030] Where: When the subscript i is equal to 1, 2, or 3, x i Represent the three directions of x, y, and z respectively; are the seepage velocities of the matrix and fracture units in the x, y, and z directions of the fracture medium in the Cartesian coordinate system, m / s; K f is the fracture permeability, m 2 ; and are the pressure gradients of the matrix along the x, y, and z directions of the grid, Pa / m; and are the pressure gradients of the cracks along the x, y, and z directions of the grid, Pa / m;

[0031] The average pressure gradient and average flow velocity in the entire calculation domain are calculated by the following formulas:

[0032]

[0033] Where: is the average pressure gradient in the entire calculation domain of the fractured medium, Pa / m; is the average flow velocity of the fractured medium in the entire calculation domain, m / s; V t is the total volume of all grids, m 3 ;∑ e is the summation symbol, indicating the summation of all elements in the computational domain; V e is the volume of each grid unit in the model, m 3 .

[0034] Preferably, in step S6, the equivalent permeability full tensor matrix is ​​calculated by the following formula:

[0035]

[0036] In the formula: superscripts Ⅰ, Ⅱ, Ⅲ represent the simulation results under three types of boundary conditions; is the equivalent permeability tensor, m 2 , when the subscript i or j is 1, 2, or 3, it represents the x, y, or z directions respectively;

[0037] The equivalent permeability full tensor matrix is ​​symmetric processed by the following formula:

[0038]

[0039] Where: i≠j.

[0040] The beneficial effects of the present invention are:

[0041] The present invention solves the full tensor of equivalent permeability of fractured media based on the EDFM (Embedded Discrete Fracture Model) method. By dividing the structured grid, the explicit processing of a large number of natural fractures is avoided, thereby reducing the computational cost of numerical simulation of oil and gas reservoirs. The problem of low computational efficiency and high computational cost caused by complex fracture grid morphology or too many fractures in fractured media can be solved. The calculated full tensor of equivalent permeability can well characterize the heterogeneity and anisotropy of fractured media, providing technical support for the development of fractured media. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0043] Figure 1 It is a schematic flow chart of the full tensor calculation method of equivalent permeability of fractured media based on EDFM of the present invention;

[0044] Figure 2 This is a schematic diagram of loading boundary conditions of the present invention;

[0045] Figure 3 It is a schematic diagram of meshing results of a fracture medium physical model in a specific implementation;

[0046] Figure 4 It is a schematic diagram of the pressure distribution results of fractured media under type I boundary conditions in a specific implementation;

[0047] Figure 5 It is a schematic diagram of the pressure distribution results of fractured media under type II boundary conditions in a specific implementation;

[0048] Figure 6 It is a schematic diagram of the pressure distribution results of fractured media under the Class III boundary condition in a specific implementation;

[0049] Figure 7 It is a schematic diagram of the full tensor ellipsoid of equivalent permeability of fractured media in a specific implementation. DETAILED DESCRIPTION

[0050] The present invention is further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, in the absence of conflict, the embodiments in this application and the technical features in the embodiments can be combined with each other. It should be noted that, unless otherwise specified, all technical and scientific terms used in this application have the same meanings as those generally understood by those of ordinary skill in the art to which this application belongs. The words "including" or "comprising" and the like used in the disclosure of the present invention mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects.

[0051] like Figure 1 As shown, the present invention provides a method for calculating the equivalent permeability full tensor of a fractured medium based on EDFM, comprising the following steps:

[0052] S1: Generate random cracks and construct a physical model of fractured media.

[0053] In a specific embodiment, when generating random fractures, they are randomly generated according to geological statistical analysis of fracture attribute parameters and representative statistical characteristics of fracture attribute parameters.

[0054] S2: Dividing the matrix grid and the fracture grid on the fracture medium physical model based on EDFM, and assigning basic grid parameters.

[0055] S3: Construct a fractured medium seepage model based on EDFM.

[0056] In a specific embodiment, the fractured medium seepage model includes a flow control equation and a crossflow term calculation equation;

[0057] The flow control equation is:

[0058]

[0059] Where: φ m ,φ f are matrix and fracture porosity, dimensionless; ρ is density, kg / m 3 ; t is time, s; represents the divergence operator, dimensionless; u m 、u f are the seepage velocities in the matrix and fractures, m / s; α f is the aperture of the crack, m; Ψ f，m , m，f are the cross-flow terms between fracture units and matrix units, and the cross-flow terms between matrix units and fracture units, respectively. (m 3 / s) / m 3 ; f1，f2 is the cross-flow term between fracture units f1 and f2, (m 3 / s) / m 2 ;

[0060] The cross-flow term calculation equation is:

[0061]

[0062] Where: L m，f is the length of the fracture element f within the matrix element m, m; <d> m，f is the average distance between matrix unit m and fracture unit f, m; K m is the matrix permeability, m 2 ; μ is the viscosity of the fluid in the grid, Pa·s; P m , P f are matrix and fracture pressure, Pa; α f1 , α f2 are the crossflow coefficients of fracture unit f1 and fracture unit f2, dimensionless; P f1 , P f2 are the pressures of fracture unit f1 and fracture unit f2, Pa, respectively;

[0063] The crossflow coefficients of the crack unit f1 and the crack unit f2 are calculated by the following formula:

[0064]

[0065] In the formula: FA f is the aperture of the fracture unit f, m; λ t，f is the total fluidity of the fracture unit f, dimensionless; dx f is the spacing of crack units, m.

[0066] It should be noted that the above cross-flow term calculation equation only lists the cross-flow term Ψ between the matrix unit and the fracture unit. m，f The calculation formula of Ψ f，m is the cross-flow term between the fracture unit and the matrix unit, and its value is related to Ψ m，f On the contrary, in m，f Add a minus sign before the calculation formula to get Ψ f，m The calculation formula for .

[0067] S4: Three types of boundary condition loading modes are adopted, and steady-state seepage numerical simulations are carried out according to the fracture medium seepage model to solve the pressure distribution in the entire calculation domain.

[0068] In a specific embodiment, Figure 2 As shown, the three types of boundary condition loading modes refer to setting boundary conditions on two surfaces corresponding to the three directions of x, v and z respectively;

[0069] In the x direction, the two corresponding surfaces are yz and y′z′; in the y direction, the two corresponding surfaces are xz and x′z′; in the z direction, the two corresponding surfaces are xy and x′y′;

[0070] The boundary condition is that along one direction, a higher pressure one is applied on one of the corresponding surfaces, a lower pressure two is applied on the other corresponding surface, the matrix grids of the remaining four surfaces serve as closed boundaries, and the pressure two is applied on the remaining four fracture unit surfaces.

[0071] In a specific embodiment, the pressure one and the pressure two can be determined according to two points of the target formation, that is, according to the pressure gradient of the target formation.

[0072] S5: According to the pressure distribution in the entire calculation domain, the pressure gradient and flow velocity distribution are calculated respectively, and then the average pressure gradient and average flow velocity in the entire calculation domain are calculated respectively by volume weighted average.

[0073] In a specific embodiment, the flow velocity distribution is calculated by the following formula:

[0074]

[0075] Where: When the subscript i is equal to 1, 2, or 3, x i Represent the three directions of x, y, and z respectively; are the seepage velocities of the matrix and fracture units in the x, y, and z directions of the fracture medium in the Cartesian coordinate system, m / s; K f is the fracture permeability, m 2 ; and are the pressure gradients of the matrix along the x, y, and z directions of the grid, Pa / m; and are the pressure gradients of the cracks along the x, y, and z directions of the grid, Pa / m;

[0076] The average pressure gradient and average flow velocity in the entire calculation domain are calculated by the following formulas:

[0077]

[0078] Where: is the average pressure gradient in the entire calculation domain of the fractured medium, Pa / m; is the average flow velocity of the fractured medium in the entire calculation domain, m / s; V t is the total volume of all grids, m 3 ;∑ e is the summation symbol, indicating the summation of all elements in the computational domain; V e is the volume of each grid unit in the model, m 3 .

[0079] S6: Based on the generalized Darcy's law, an equivalent permeability full tensor matrix is ​​calculated and obtained, and the equivalent permeability full tensor matrix is ​​symmetrically processed to obtain a calculation result of the equivalent permeability full tensor of the fractured medium based on EDFM.

[0080] In a specific embodiment, the equivalent permeability full tensor matrix is ​​calculated by the following formula:

[0081]

[0082] In the formula: superscripts Ⅰ, Ⅱ, Ⅲ represent the simulation results under three types of boundary conditions; is the equivalent permeability tensor, m 2 , when the subscript i or j is 1, 2, or 3, it represents the x, y, or z directions respectively;

[0083] The equivalent permeability full tensor matrix is ​​symmetric processed by the following formula:

[0084]

[0085] Where: i≠j.

[0086] In a specific embodiment, taking a certain reservoir as an example, the equivalent permeability full tensor of the fractured medium of the reservoir is calculated by using the EDFM-based fractured medium equivalent permeability full tensor calculation method of the present invention, which specifically includes the following steps:

[0087] (1) Generate random cracks and construct a physical model of fractured media.

[0088] In this embodiment, random fractures are randomly generated based on geological statistical analysis of fracture attribute parameters and representative statistical characteristics of fracture attribute parameters.

[0089] (2) Dividing the matrix grid and the fracture grid on the fracture medium physical model based on EDFM, and assigning basic grid parameters.

[0090] In this embodiment, the grid division result is as follows: Figure 3 As shown, the basic grid parameters include reservoir matrix permeability 0.1mD, matrix porosity 0.25, and fracture permeability 1D.

[0091] (3) Based on EDFM, a fractured medium seepage model as shown in equations (1)-(5) is constructed.

[0092] (4) Three types of boundary condition loading modes are adopted, and steady-state seepage numerical simulations are carried out according to the fracture medium seepage model to solve the pressure distribution in the entire calculation domain.

[0093] In this embodiment, the pressure 1 and the pressure 2 are determined according to the pressure gradient of the target reservoir, the pressure 1 is 5MPa, and the pressure 2 is 1MPa. The first type of boundary condition loading is along the z direction, the second type of boundary condition loading is along the x direction, and the third type of boundary condition loading is along the y direction. The fracture medium pressure distribution simulated by the fracture medium physical model after loading three types of boundary conditions is as follows: Figure 4-Figure 6 As shown. Figure 4-Figure 6 It can be seen that the pressure distribution is uneven due to the presence of a large number of natural fractures in the fractured medium. When the pressure distribution is uneven, the seepage will be affected by the fractures, resulting in uneven permeability distribution. The pressure gradient of the entire calculation domain can be obtained through the calculated pressure distribution, thereby calculating the equivalent permeability full tensor.

[0094] (5) According to the pressure distribution in the entire calculation domain, the pressure gradient and flow velocity distribution are calculated respectively, and then the average pressure gradient and average flow velocity in the entire calculation domain are calculated respectively by volume weighted average.

[0095] In this embodiment, the flow velocity distribution is calculated by equations (6)-(7), and the average pressure gradient and average flow velocity in the entire calculation domain are calculated by equations (8)-(9).

[0096] (6) Based on the generalized Darcy's law, the equivalent permeability full tensor matrix is ​​calculated by equation (10), and the equivalent permeability full tensor matrix is ​​symmetrically processed by equations (11)-(12) to obtain the calculation result of the equivalent permeability full tensor of the fractured medium based on EDFM.

[0097] In this embodiment, the calculation result of the equivalent permeability full tensor of the fractured medium based on EDFM is:

[0098]

[0099] In fractured media, the direction and velocity of fluid flow can be described by the permeability tensor. In this embodiment, the equivalent permeability tensor of fractured media calculated based on EDFM is a third-order tensor. By solving the eigenvalue equation of the third-order square matrix, three permeabilities are obtained, namely, the maximum principal permeability k max =9.2×10 -15 , minimum principal permeability k min =2.1×10 -15 and the average principal permeability k med =5.53×10 -15 , according to the three main permeabilities calculated, Figure 7 The full tensor ellipsoid diagram of equivalent permeability of fractured media is shown, through which the calculation results can be displayed intuitively.

[0100] In summary, the present invention can quickly calculate the equivalent permeability full tensor of fractured media based on EDFM, which can be used to characterize the heterogeneity and anisotropy of fractured media. Compared with the prior art, the present invention has significant progress.

[0101] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modify the technical contents disclosed above into equivalent embodiments 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 according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.< / d> < / d>

Claims

1. A full tensor calculation method for equivalent permeability of fractured media based on EDFM, characterized in that: The following steps are involved: S1: Generate random cracks and build a physical model of fractured media; S2: dividing the matrix grid and the fracture grid on the fracture medium physical model based on EDFM, and assigning basic grid parameters; S3: Construct fracture medium seepage model based on EDFM; S4: adopting three types of boundary condition loading modes, and carrying out steady-state seepage numerical simulations according to the fracture medium seepage model to solve the pressure distribution in the entire calculation domain; S5: According to the pressure distribution in the entire calculation domain, the pressure gradient and the flow velocity distribution are calculated respectively, and then the average pressure gradient and the average flow velocity in the entire calculation domain are calculated respectively by volume weighted average; S6: Based on the generalized Darcy's law, an equivalent permeability full tensor matrix is ​​calculated and obtained, and the equivalent permeability full tensor matrix is ​​symmetrically processed to obtain a calculation result of the equivalent permeability full tensor of the fractured medium based on EDFM.

2. The method for calculating equivalent permeability of fractured media based on EDFM according to claim 1, characterized in that: In step S1, when generating random fractures, they are randomly generated according to geological statistical analysis of fracture attribute parameters and representative statistical characteristics of fracture attribute parameters.

3. The method for calculating equivalent permeability of fractured media based on EDFM according to claim 1, characterized in that: In step S3, the fractured medium seepage model includes a flow control equation and a crossflow term calculation equation; The flow control equation is: Where: φ m ,φ f are matrix and fracture porosity, respectively, dimensionless; ρ is density, kg / m 3 ; t is time, s; represents the divergence operator, dimensionless; u m 、u f are the seepage velocities in the matrix and fractures, m / s; a f is the aperture of the crack, m; Ψ f,m , m,f are the cross-flow terms between fracture units and matrix units, and the cross-flow terms between matrix units and fracture units, respectively. (m 3 / s) / m 3 ; f1,f2 is the cross-flow term between fracture units f1 and f2, (m 3 / s) / m 2 ; The cross-flow term calculation equation is: Where: L m,f is the length of the fracture element f within the matrix element m, m; <d> m,f is the average distance between matrix unit m and fracture unit f, m; K m is the matrix permeability, m 2 ; μ is the viscosity of the fluid in the grid, Pa·s; P m , P f are matrix and fracture pressure, Pa; α f1 , α f2 are the crossflow coefficients of the fracture unit f1 and the fracture unit f2, dimensionless;< / d> P f1 , P f2 are the pressures of fracture unit f1 and fracture unit f2, Pa, respectively; The crossflow coefficients of the crack unit f1 and the crack unit f2 are calculated by the following formula: In the formula: FA f is the aperture of the fracture unit f, m; λ t,f is the total fluidity of the fracture unit f, dimensionless; dx f is the spacing of crack units, m.

4. The method for calculating equivalent permeability of fractured media based on EDFM according to claim 1, characterized in that: In step S4, the three types of boundary condition loading modes refer to setting boundary conditions on two surfaces corresponding to the three directions of x, y and z respectively; In the x direction, the two corresponding surfaces are yz and y'z'; in the y direction, the two corresponding surfaces are xz and x'z'; in the z direction, the two corresponding surfaces are xy and x'y'; The boundary condition is that along one direction, a higher pressure one is applied on one of the corresponding surfaces, a lower pressure two is applied on the other corresponding surface, the matrix grids of the remaining four surfaces serve as closed boundaries, and the pressure two is applied on the remaining four fracture unit surfaces.

5. The method for calculating equivalent permeability full tensor of fractured media based on EDFM according to claim 3, characterized in that: In step S5, the flow velocity distribution is calculated by the following formula: Where: When the subscript i is equal to 1, 2, or 3, x i Represent the three directions of x, y, and z respectively; are the seepage velocities of the matrix and fracture units of the fractured medium in the x, y, and z directions in the Cartesian coordinate system, m / s; K f is the fracture permeability, m 2 ; and are the pressure gradients of the matrix along the x, y, and z directions of the grid, Pa / m; and are the pressure gradients of the cracks along the x, y, and z directions of the grid, Pa / m; The average pressure gradient and average flow velocity in the entire calculation domain are calculated by the following formulas: Where: is the average pressure gradient in the entire calculation domain of the fractured medium, Pa / m; is the average flow velocity of the fractured medium in the entire calculation domain, m / s; V t is the total volume of all grids, m 3 ; ∑ e is the summation symbol, indicating the summation of all elements in the computational domain; V e is the volume of each grid unit in the model, m 3 .

6. The method for calculating equivalent permeability of fractured media based on EDFM according to claim 5, characterized in that: In step S6, the equivalent permeability full tensor matrix is ​​calculated by the following formula: In the formula: superscripts Ⅰ, Ⅱ, Ⅲ represent the simulation results under three types of boundary conditions; is the equivalent permeability tensor, m 2 , when the subscript i or j is 1, 2, or 3, it represents the x, y, or z directions respectively; The equivalent permeability full tensor matrix is ​​symmetric processed by the following formula: Where: i≠j.

Images