Projection-based embedded discrete fracture model using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method

Abstract

This invention presents a projection embedded discrete fracture model integrating a TPFA and MFD hybrid approach, creating a pEDFM framework for various anisotropic two-phase flow situations. It specifies the distribution of extra pressure freedoms on matrix grids for MFD implementation, maintains f-f connections in TPFA through a standard pEDFM workflow, and introduces a low-conductivity fracture treatment for MFD. It also outlines the derivation of numerical flux calculation formulas for effective m-m and m-f connections. The mixed TPFA-MFD design applies to numerical flux estimation across both K-orthogonal and non-K-orthogonal grids, enhancing computational efficiency and facilitating the spatial discretization of continuity equations for matrix and fracture grids under anisotropic permeability conditions. A global equation system is formulated based on the continuity of effective connections, with time discretization via the implicit backward Euler method and pressure and water saturation distributions determined by a Newton-Raphson based nonlinear solver.

Claims

1. A method for generating flow simulation of fracture formation, for an anisotropic media with full permeability tensor, by a projection-based embedded discrete fracture model (pEDFM) using a hybrid of two-point flux approximation (TPFA) method and mimetic finite difference (MFD) method, i.e., (TPFA-MFD) method, the method comprising: discretizing of a generic fracture model of the anisotropic media, by a data processing system, for: identifying a plurality of grids comprising matrix grids and fracture grids based on a location of each grid in the generic fracture model; identifying three types of connections in a fracture network based on the location of each grid in the generic fracture model, wherein the three types of connections comprises: m-m connection corresponding to connection between two matrix grids; m-f connection corresponding to connection between a matrix grid and a fracture grid; and f-f connection corresponding to a connection between two fracture grids; and identifying effective m-m connections and/or effective m-f connections, wherein each effective m-m connection and/or each effective m-f connection comprises at least one low-conductivity fracture f.sub.l; constructing a treatment method for each low-conductivity fracture by the data processing system, for: identifying an inter-grid connection set Cont.sub.1 comprising a set of connections obtained from a generic pEDFM to treat the low-conductivity fractures with a non-projection transmissibility multiplier method is used that provides a flow area between adjacent grids, wherein the low-conductivity fractures are ignored; identifying an efficient set of connections Cont.sub.1.sup.eff based on the inter-grid connection set Cont.sub.1 as follows:
Cont.sub.1.sup.eff=Cont.sub.1Cont.sub.1.sup.0(6) wherein Cont.sub.1.sup.0 is a set of connections with a flow area 0 in the inter-grid connection set Cont.sub.1; constructing a numerical flux calculation method for effective m-f and m-m connections, by the data processing system, for: determining flux continuity conditions of effective m-f and m-m connections using the hybrid TPFA-MFD method, wherein the hybrid TPFA-MFD method comprises: estimating a numerical flux on each effective m-f connection and/or effective m-m connection related to K-orthogonal grids using TPFA method; and estimating a numerical flux on each effective m-f connection and/or effective m-m connection related to non-K-orthogonal grids using MFD method; and performing a spatial discretization of a continuity equation of each matrix grid and each fracture grid in the effective m-m and m-f connections; and constructing global equations, by the data processing system, by combining the flux continuity conditions of the effective m-f and m-m connections and a time discretization scheme of implicit backward Euler scheme, and calculating distribution of pressure and water saturation by solving the global equations using a nonlinear solver based on Newton-Raphson method.

2. The method for generating the flow simulation of fracture formation, for the anisotropic media with full permeability tensor, by the projection-based embedded discrete fracture model (pEDFM) using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method according to claim 1, wherein the method for constructing the treatment method for each of the low-conductivity fractures by the data processing system comprises, determining a transmissibility of each fracture grid in pEDFM based on TPFA method, wherein the at least one low-conductivity fracture f.sub.l has a permeability k.sub.f and a fracture width w.sub.f, wherein determining transmissibility of fracture grids in the pEDFM includes determining the transmissibility in three types of connections that are affected by low-conductivity fractures, wherein the transmissibility based on low-conductivity fractures for the three types of connections are as follows: (1) for f-f connections, set to connect between i-th fracture grid f.sub.i and j-th fracture grid f.sub.j, determining a transmissibility of fracture grids in f.sub.i-f.sub.j connection that is affected by low-conductivity fracture f.sub.l comprises: determining a first transmissibility of f.sub.i-f.sub.j connection before being affected by the low-conductivity fracture grid f.sub.l as:
T.sub.f.sub.i.sub.f.sub.j=(T.sub.f.sub.i.sup.1+T.sub.f.sub.j.sup.1)(1) wherein, T.sub.f.sub.i is a transmissibility of fracture grid f.sub.i, T.sub.fj is a transmissibility of fracture grid f.sub.j; and determining a transmissibility of f.sub.i-f.sub.j connection after being affected by low-conductivity fracture f.sub.1 as: $\begin{array}{cc}{T}_{f_i{f}_j}={\left(T_{f_i{f}_j}^{-1}+T_{f_l}^{-1}\right)}^{-1},\mathrm{wherein}{T}_{f_l}=\frac{k_{f_l}A\left(f_i,f_j\right)}{\frac{w_{f_l}}{2}}& \left(2\right)\end{array}$ wherein, k.sub.fi is a permeability of fracture grid f.sub.1, w.sub.fi is a width of an opening of fracture grid f.sub.1, A.sub.(f.sub.i.sub., f.sub.j.sub.) is a flow area of (f.sub.i, f.sub.j); and updating the transmissibility T.sub.f.sub.i.sub.f.sub.j using the Eq. (1) if other low-conductivity fractures f.sub.n that affect the f-f connections exist, wherein, for the low-conductivity fracture grid f.sub.n that intersects (f.sub.i, f.sub.j) connection, the transmissibility is further updated based on Equation (2) as: $\begin{array}{cc}{T}_{f_i{f}_j}={\left(T_{f_i{f}_j}^{-1}+T_{f_n}^{-1}\right)}^{-1},& \left(2\right)\end{array}$ (2) for m-f connections, set to connection between i-th matrix grid m.sub.i and j-th fracture grid f.sub.j, determining a transmissibility of fracture grid in m.sub.i-f.sub.j connection that is affected by low-conductivity fracture f.sub.j comprises: identifying the transmissibility of j-th fracture grid f.sub.j in m.sub.i-f.sub.j connections not affected by low-conductivity fracture grid f.sub.l based on the transmissibility T.sub.fj obtaining the transmissibility T.sub.fj of j-th fracture grid f.sub.j after applying the effect of low-conductivity fracture f.sub.l is-applied to the j-th fracture grid f.sub.j, which is equivalent to the series connection of low-conductivity fractures f.sub.1 and f.sub.j, wherein, the transmissibility T.sub.fj of fracture grid f.sub.j is updated as:
T.sub.f.sub.j=(T.sub.f.sub.j.sup.1+T.sub.f.sub.l.sup.1).sup.1(3); and updating the transmissibility Td using Eq. (3) if other low-conductivity fractures f.sub.n that affect the m.sub.i-f.sub.j connections exist, wherein for the other low-conductivity fracture grid f.sub.n that intersect with m.sub.i-f.sub.j connection, the transmissibility T.sub.fj is further updated based on Equation (3) as follows; and (3) for m-m connections, set to connection between i-th matrix grid m.sub.i and j-th matrix grid m.sub.j, determining a transmissibility of fracture grid in the m.sub.i-m.sub.j connections that is affected by low-conductivity fracture f.sub.l comprises: splitting the m.sub.i-m.sub.j connections into m.sub.i-f.sub.l connections and m.sub.j-f.sub.l connections, wherein, flow area of m.sub.i-f.sub.l connections and m.sub.j-f.sub.l connections is same as a flow area A.sub.(mi, mj) of (m.sub.i, m.sub.j), wherein a transmissibility of the low conductivity fracture fi is as follows: $\begin{array}{cc}{T}_{f_l}^{}=\frac{k_{f_l}{A}_{\left(m_i,m_j\right)}}{\frac{w_{f_l}}{2}}& \left(4\right)\end{array}$ and updating the transmissibility T.sub.fj of the low-conductivity fracture f.sub.l using Eq. (3) if other low-conductivity fractures f.sub.k that affect (m.sub.i, m.sub.j) exist, wherein the transmissibility of the low-conductivity fracture f.sub.l is updated by using a transmissibility of the other low-conductivity fracture f.sub.k in Eq. (3) as follows: $\begin{array}{cc}{T}_{f_l}={\left(T_{f_l}^{-1}+T_{f_k}^{-1}\right)}^{-1},\mathrm{wherein}{T}_{f_k}=\frac{k_{f_k}A\left(m_i{m}_j\right)}{\frac{w_{f_k}}{2}}.& \left(5\right)\end{array}$

3. The method for generating flow simulation of fracture formation, for the anisotropic media with full permeability tensor, by the projection-based embedded discrete fracture model (pEDFM) using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method according to claim 1, wherein constructing the numerical flux calculation method for the effective m-f and m-m connections, by the data processing system comprises: identifying a number set of a matrix grid which is adjacent to a matrix grid m.sub.i to be neigh.sub.m.sub.i, wherein a set of effective connections on each m.sub.i-side face of matrix grid m.sub.i is: Cont.sub.2.sup.eff(m.sub.i)=.sub.jneigh.sub.m.sub.i.sup. Cont.sub.2.sup.eff(l.sub.ij.sup.m.sup.i) according to Eq. (7): $\begin{array}{cc}{\stackrel{}{p}}_{{}_{ij}}^{m_i}=\frac{{.Math.}_{{{\mathrm{Cont}}_2^{\mathrm{eff}}\left({}_{ij}^{m_i}\right)}^{p_{}{A}_{}}}}{{.Math.}_{{{\mathrm{Cont}}_2^{\mathrm{eff}}\left({}_{ij}^{m_i}\right)}^{A_{}}}}& \left(7\right)\end{array}$ wherein: l.sub.ij.sup.m.sup.i is the m.sub.i-side of l.sub.ij, wherein l.sub.ij; is the interface between matrix grid m.sub.i and matrix grid m.sub.j, Cont.sub.2.sup.eff(m.sub.i) denotes the set of connections in Cont.sub.2.sup.eff that have additional pressure degree of freedom located on all m.sub.i-side faces of matrix grid m.sub.i in addition to Cont.sub.1.sup.eff, Cont.sub.2.sup.eff(l.sub.ij.sup.m.sup.i) denotes the set of connections in Cont.sub.2.sup.eff that have additional pressure degree of freedom is located on l.sub.ij.sup.m.sup.i in addition to Cont.sub.1.sup.eff, p.sub.l.sub.ij.sup.m.sup.i is an average pressure of a side near matrix grid m.sub.i on intersecting face of matrix grids m.sub.i and m.sub.j, p.sub. is a pressure degree of freedom added by th connection in Cont.sub.2.sup.eff(l.sub.ij.sup.m.sup.i), and A.sub. is a flow area of the th connection in Cont.sub.2.sup.eff(l.sub.ij.sup.m.sup.i); obtaining a conversion formula between the average pressure of m.sub.i-side of m.sub.i-faces and the value of the pressure degree of freedom added by each connection in Cont.sub.2.sup.eff(m.sub.i) as follows:
pm.sub.i=Ap.sub.m.sub.i(8) wherein: p.sub.m.sub.i is a column vector composed of p.sub.ij.sup.m.sup.i(jneigh.sub.m.sub.i), p.sub.ij.sup.m.sup.i is the average pressure of l.sub.ij.sup.m.sup.i, p.sub.ij.sup.m.sup.i, is a column vector composed of p.sub.j(jCont.sub.2.sup.eff(m.sub.i)) wherein, if jCont.sub.2.sup.eff, (m.sub.i) is a connection in Cont.sub.2.sup.eff(m.sub.i), then p.sub.j denotes the added pressure degree of freedom for this connection, and each j-line of A represents an area weight when using Eq. (7) to calculate the average pressure on l.sub.ij.sup.m.sup.i (i.e., the m.sub.i-side of the j-th surface of m.sub.i).

4. The method for generating flow simulation of fracture formation, for the anisotropic media with full permeability tensor, by the projection-based embedded discrete fracture model (pEDFM) using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method according to claim 3, wherein, constructing the numerical flux calculation method for the effective m-f and m-m connections, by the data processing system comprises: using a distribution matrix A.sup.T to relate a numerical flux on each surface of m.sub.i denoted by flux.sub.i to a numerical flux flux.sub.i on each effective connection in Cont.sub.2.sup.eff(m.sub.i) as follows:
flux.sub.i=A.sup.Tflux.sub.i(9), wherein: if the matrix grid m.sub.i is K orthogonal, then a transmissibility matrix T.sup.m.sup.i in the i-th matrix grid m.sub.i is calculated by using Eq. (10) based on TPFA as follows: $\begin{array}{cc}{T}_{,y}^{m_i}=\{\begin{array}{c}{T}_i\mathrm{if}=\\ 0\mathrm{if}\end{array},\mathrm{wherein}{T}_i=.Math.{}_i.Math.\frac{K_i{r}_i}{{.Math.r_i.Math.}^2}.Math.n_i& \left(10\right)\end{array}$ wherein: T.sub.,.sup.m.sup.i represents row and column of matrix T.sup.m.sup.i, .sub.i is portion of a control region of grid m.sub.i at -th row and -th column of the transmissibility matrix T.sub.,.sup.m.sup.i, |.sub.i| is volume of the portion of the control region of grid m.sub.i at -th row and -th column of the transmissibility matrix T.sub.,.sup.m.sup.i, r.sub.i is a vector from grid i-center to grid .sub.i-center, and .sub.i is a unit outward normal vector of .sub.i; and if the matrix grid m.sub.i is non-K orthogonal, then the transmissibility matrix T.sup.m.sup.i in m.sub.i is calculated by using Eq. (11) based on MFD as follows: $\begin{array}{cc}{T}_{i1}=\frac{1}{.Math.{}_i.Math.}{N}_i{K}_i{N}_i^T{T}_{i2}=\frac{6}{d}tr\left(K_i\right)A_i\left(I_i-Q_i{Q}_i^T\right)A_i& \left(11\right)\end{array}$ wherein: .sub.i is a control region of matrix grid m.sub.i, |.sub.i| is a volume of matrix grid m.sub.i, K.sub.i is a permeability tensor of matrix grid m.sub.i, X.sub.i=(x.sub.i1x.sub.i, . . . , x.sub.ix.sub.i, . . . , x.sub.in.sub.ix.sub.i).sup.T, x.sub.i is a vector of the center of grid i, x.sub.i is a vector of .sub.icenter, N.sub.i=(|.sub.i1|.sub.i1, . . . , |.sub.i|.sub.i, . . . , |.sub.in.sub.i|.sub.in.sub.i).sup.T wherein N.sub.i is the number of all the interfaces of the i-th grid, d is a grid dimension, tr denotes to calculate the trace of a matrix, A.sub.i=diag(|.sub.i1|, . . . , .sub.in.sub.i|)), and Q.sub.i=orth(A.sub.iX.sub.i) wherein Q.sub.i is a standard orthogonal basis for a column space of A.sub.iX.sub.i, wherein the number of columns in Q.sub.i is equal to the rank of A.sub.iX.sub.i; using an average pressure to participate in the calculation of the numerical flux flux.sub.i on each surface of the matrix grid, as:
flux.sub.i=T.sup.m.sup.ip.sub.m.sub.i=T.sup.m.sup.iAp.sub.m.sub.i=T.sup.m.sup.iA(p.sub.m.sub.iIp.sub.m.sub.i)(12) wherein: flux.sub.i is a column vector composed of numerical flux on each surface of the matrix grid, each j-line of A represents the area weight when the average pressure on the m.sub.i-side of the j-th surface of m.sub.i is calculated by Eq. (8), A.sup.T is a transverse of A and is a distribution matrix to relate the numerical flux on each surface of matrix grid m.sub.i to the flux on each connection in Cont.sub.2.sup.eff(m.sub.i), and I is a column vector with all elements being 1, and its length is equal to the number of effective connections in the matrix grid, obtaining the numerical flux and actual transmissibility matrix on the matrix grid m.sub.i by combining Eq. (9) and Eq. (12), as follows:
flux.sub.i=A.sup.TT.sup.m.sup.iA(p.sub.m.sub.iIp.sub.m.sub.i)(13), wherein flux.sub.i is the numerical flux on the matrix grid m.sub.i, and for the matrix grid m.sub.i, the actual transmissibility matrix is:
{tilde over (T)}.sup.m.sup.i=A.sup.TT.sup.m.sup.iA(14).

5. The method for generating flow simulation of fracture formation, for the anisotropic media with full permeability tensor, by the projection-based embedded discrete fracture model (pEDFM) using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method according to claim 4, wherein constructing global equations by the data processing system comprises: calculating the flow on each connection related to m.sub.i by using the actual transmissibility matrix given in Eq. (14), and obtaining a discrete scheme of the continuity equation in matrix grid m.sub.i as; $\begin{array}{cc}{.Math.}_{{\mathrm{Cont}}_1^{\mathrm{eff}}\left(m_i\right)}\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}{.Math.}_{Con{t}_1^{\mathrm{eff}}\left(m_i\right)}{\tilde{T}}_{,}^{m_i}\left(p_{m_i^{t+t}}-p_{}^{t+t}\right)+Q_a=\frac{.Math.{}_i.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}_{m_i}^{t+t}-{\left(\frac{S_a}{B_a}\right)}_{m_i}^t\right]& \left(15\right)\end{array}$ wherein: k.sub.r,ij is a relative permeability of the -th phase between the matrix grid m.sub.i and the matrix grid m.sub.j according to a single-point upwind scheme, .sub.,ij and B.sub.,ij are viscosity and volume coefficient of the ith phase calculated by an arithmetic average scheme between the matrix grid m.sub.i and the matrix grid m.sub.j, respectively, subscript refers to a serial number of an intersection surface of matrix grid m.sub.i and matrix grid m.sub.j in all surfaces of matrix grid m.sub.i, F.sub.i is an outward normal flux of matrix grid m.sub.i on the -th plane, T.sub.,.sup.m.sup.i is the -th row and -th column of the actual transmissibility matrix {tilde over (T)}.sup.m.sup.i, is the surface center pressure of the -th plane of the matrix grid m.sub.i, is a body center pressure of matrix grid m.sub.i, t is a time stepping, .sub.i, S.sub.,i and B.sub.,i are a porosity of matrix grid m.sub.i, a saturation of -th phase and the volume coefficient of -th phase, respectively, and superscripts t+t and t represent the time; and obtaining a discrete scheme of continuity equation for the fracture grid f.sub.j, wherein when the transmissibility is in a simple scheme based on TPFA, the effective connection Cont.sub.2.sup.eff(f.sub.j) related to f.sub.j includes f-f connection Cont.sub.f.sub.j.sup.f-f and m-f connection Cont.sub.f.sub.j.sup.m-f, and if the set of fracture grids adjacent to f.sub.j reflected from Cont.sub.f.sub.j.sup.f-f is denoted as neigh.sub.f.sub.j, then the discrete scheme of the continuity equation in fracture grid f.sub.j is: $\begin{array}{cc}\underset{Con{t}_{f_j}^{m-f}}{.Math.}\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}\underset{Con{t}_{f_j}^{m-f}}{.Math.}{\tilde{T}}_{,}^{m_i}\left(p_{m_i^{t+t}}-p_{}^{t+t}\right)+\underset{{\mathrm{neigh}}_{f_j}}{.Math.}\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}\left(p_{f_j}^{t+t}-p_{}^{t+t}\right)+Q_a=\frac{.Math.{}_{f_j}.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}_{f_j}^{t+t}-{\left(\frac{S_a}{B_a}\right)}_{f_j}^t\right].& \left(16\right)\end{array}$

6. The method for generating flow simulation of fracture formation, for the anisotropic media with full permeability tensor, by the projection-based embedded discrete fracture model using hybrid of two-point flux approximation and mimetic finite difference (TPFA-MFD) method according to claim 5, wherein when actual transmissibility of each matrix grid {tilde over (T)}.sup.m.sup.i is obtained, the transmissibility of all matrix grids in effective m-m connections and m-f connections is known, while a transmissibility of fracture grids in m-f connections still adopts a simple scheme based on TPFA wherein: for m-m connections, suppose m.sub.i and m.sub.k are matrix grids corresponding to i-th surface and k-th surface respectively, then, the corresponding numerical flux is:
flux.sub.m.sub.i.sub..fwdarw.m.sub.k={tilde over (T)}.sup.m.sup.i.sub.m.sub.k(p.sub.m.sub.iIp.sub.m.sub.i)flux.sub.m.sub.k.sub..fwdarw.m.sub.i={tilde over (T)}.sup.m.sup.i.sub.m.sub.k(p.sub.m.sub.iIp.sub.m.sub.i)(17) the corresponding flux continuity condition is:
{tilde over (T)}.sup.m.sup.i.sub.m.sub.k(p.sub.m.sub.iIp.sub.m.sub.i)+{tilde over (T)}.sup.m.sup.i.sub.m.sub.k(p.sub.m.sub.iI-p.sub.m.sub.i)=0(18) for m-f connections, suppose m.sub.i and f.sub.j are the matrix grid corresponding to i-th surface and the fracture grid corresponding to i-th surface respectively, then, the corresponding numerical flux is:
flux.sub.m.sub.i.sub..fwdarw.f.sub.j={tilde over (T)}.sup.m.sup.i.sub.f.sub.j(p.sub.m.sub.iIp.sub.m.sub.i)flux.sub.f.sub.j.sub..fwdarw.m.sub.i=T.sub.f.sub.j(p.sub.f.sub.jp.sub.(m.sub.i.sub.,f.sub.j.sub.))(19), and the corresponding flux continuity condition is:
{tilde over (T)}.sup.m.sup.i.sub.f.sub.j(p.sub.m.sub.iIp.sub.m.sub.i)+T.sub.f.sub.j(p.sub.f.sub.jp.sub.(m.sub.i.sub.,f.sub.j.sub.))=0(20); and for f-f connections, the transmissibility is obtained using the generic pEDFM without defining additional pressure degrees of freedom is adopted.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a model diagram;

(2) FIG. 2 is a sketch of some cases;

(3) FIG. 3 is a Cartesian background grid corresponding to the reservoir model and the reservoir model of example 1;

(4) FIG. 4 shows a water saturation distribution at 200 days calculated by different methods in case 1 of the example 1;

(5) FIG. 5 shows a water saturation distribution at 200 days calculated by different methods in case 2 of the example 1;

(6) FIG. 6 is a pressure distribution of 200 days calculated by different methods in case 2 of the example 1;

(7) FIG. 7 is a water cut curve of production wells in case 1 and case 2 of the example 1;

(8) FIG. 8 is a distribution of water saturation at 200 days calculated by different methods in case 2 of the example 2;

(9) FIG. 9 shows a pressure distribution at 200 days calculated by different methods of the example 2;

(10) FIG. 10 is a water cut curve of production wells calculated by different methods of the example 2;

(11) FIG. 11 is a reservoir model and grid of the new pEDFM and DFM of the example 5;

(12) FIG. 12 shows a water saturation and pressure distribution after 400 days and 600 days obtained by different calculation methods of the example 5; and

(13) FIG. 13 is a comparison of water saturation and pressure distribution after 200 days and 600 days obtained by different calculation methods in case 2 of the example 5.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(14) In order to illustrate the technical effect of the invention, the existing technology and the improvement of the invention are explained first.

(15) 1.1 Reservoir Two-Phase Flow Control Equation

(16) Continuity equation of each phase:

(17) 0$\begin{array}{cc}-.Math.v_a+q_{a,\mathrm{sc}}=\frac{}{t}\left(\frac{S_a}{B_a}\right),a=o\left(\mathrm{oil}\right)\mathrm{or}w\left(\mathrm{water}\right)& \left(1\right)\end{array}$

(18) Wherein v.sub. is the flow velocity, q.sub.,sc is the source and sink term of phase a under the ground condition, S.sub. and B.sub. are the saturation and volume coefficient of phase a, respectively, t is time, is porosity.

(19) The flow velocity satisfies the Darcy's law:

(20) $\begin{array}{cc}{v}_a=-{}_{a}K{p}_a{}_a=\frac{k_{ra}}{B_a{}_a}& \left(2\right)\end{array}$

(21) Wherein K is the permeability tensor, .sub. is the mobility of phase , p.sub., k.sub.r, .sub. and B.sub. are the pressure, relative permeability, viscosity and volume coefficient of phase , respectively.

(22) Taking Eq. (2) into Eq. (1), get:

(23) $\begin{array}{cc}.Math.\left(\frac{k_a}{{}_a{B}_a}K{p}_a\right)+q_{a,\mathrm{sc}}=\frac{}{t}\left(\frac{S_a}{B_a}\right)& \left(3\right)\end{array}$
1.2 Discretization of Control Equations

(24) Integrating both sides of Eq. (3) in grid i control volume .sub.i and time period [t, t+Vt], get:

(25) $\begin{array}{cc}{}_t^{t+t}{}_{{}_i}.Math.\left(\frac{k_a}{{}_a{B}_a}K{p}_a\right)d\mathrm{dt}+Q_a={}_{{}_i}{\left(\frac{S_a}{B_a}\right)}^{t+t}-{\left(\frac{S_a}{B_a}\right)}^{t}d& \left(4\right)\end{array}$
wherein .sub..sub.iq.sub.,scdA=Q.sub.,sc

(26) The rectangular formula is used to estimate the time integral on the left side and the space integral on the right side of Eq. (4), get:

(27) $\begin{array}{cc}{}_{{}_i}.Math.\left(\frac{k_a}{{}_a{B}_a}K{p}_a\right)d+Q_a=\frac{.Math.{}_i.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}_i^{t+t}-{\left(\frac{S_a}{B_a}\right)}_i^t\right]& \left(5\right)\end{array}$

(28) By using the divergence theorem, the first term on the left side of (5) can be rewritten as follows:

(29) $\begin{array}{cc}{}_{{}_i}.Math.\left(\frac{k_a}{{}_a{B}_a}Kp\right)d=\underset{=1}{\overset{n_i}{.Math.}}{}_{{}_i}\frac{k_{\mathrm{ra}}}{{}_a{B}_a}Kp.Math.\mathrm{ndA}=-\underset{=1}{\overset{n_i}{.Math.}}{F}_i& \left(6\right)\end{array}$ wherein, subscript i represents the th edge of grid i (if it is three-dimensional, it is the th surface of grid i), and F.sub.i is the outward normal flux on .sub.i.

(30) For two adjacent grids, it may be written as grid i and grid j, and the intersection of grid i and grid j is written as e.sub.ij, , s.t., .sub.i=.sub.j=e.sub.ij. In Eq. (5), the single-point upstream weight scheme in (6) is generally used in k.sub.r, and .sub. and B.sub. generally adopt the arithmetic average scheme in Eq. (7).

(31) $\begin{array}{cc}{k}_{\mathrm{ra},\mathrm{ij}}=\{\begin{array}{cc}{k}_{\mathrm{ra}}\left(S_{a,j}\right)& \mathrm{if}{F}_i<0\\ k_{\mathrm{ra}}\left(S_{a,i}\right)& \mathrm{if}{F}_{i}0\end{array}& \left(7\right)\end{array}$

(32) $\begin{array}{cc}{}_{a,\mathrm{ij}}=\frac{{}_{a,i}+{}_{a,j}}{2}{B}_{a,\mathrm{ij}}=\frac{B_{a,i}+B_{a,j}}{2}& \left(8\right)\end{array}$
2.1 Two-Point Flux Approximation (TPFA)

(33) In the finite volume method with two-point flux approximation, the numerical flux of grid i on e.sub.ij is calculated as:

(34) $\begin{array}{cc}{F}_i=\frac{k_{\mathrm{ra},\mathrm{ij}}}{{}_{a,ij}{B}_{a,ij}}{T}_i\left(p_i-p_i\right)& \left(9\right)\end{array}$ wherein

(35) $T_i=.Math.{}_i.Math.\frac{K_i{r}_i}{{.Math.r_i.Math.}^2}.Math.n_i,$
r.sub.i is a vector from grid i-center to grid, .sub.i-center, n.sub.i is the unit outward normal vector of .sub.i, and p.sub.i is the pressure value at the center of .sub.i.

(36) According to Eq. (8), the expression of the outward normal flux of .sub.i expressed by the transmissibility matrix T.sup.i can be obtained:

(37) 0$\begin{array}{cc}{F}_i=\frac{k_{\mathrm{ra},\mathrm{ij}}}{{}_{a,ij}{B}_{a,ij}}\underset{=1}{\overset{n_i}{.Math.}}{T}_{,}^i\left(p_i-p_i\right)T_{,}^i=\{\begin{array}{cc}{T}_i& \mathrm{if}=\\ 0& \mathrm{if}\end{array}& \left(10\right)\end{array}$ wherein, T denotes the row and columns of matrix T.sup.i.
2.2 Simulating Finite Difference

(38) Generally, TPFA can only achieve high-precision flux approximation in K-orthogonal grids, and K-orthogonal meets that:
K.sub.ir.sub.in.sub.i,=0(11)

(39) MFD can be applied to both K-orthogonal and non-K-orthogonal grids, the transmissibility matrix T.sup.i given by MFD is generally more denser than that given by TPFA, which is different from the diagonal transmissibility matrix in (10).

(40) In MFD, the calculation expression of T.sup.i is:

(41) $\begin{array}{cc}{T}^i=T_{i1}+T_{i2}{T}_{i1}=\frac{1}{.Math.{}_i.Math.}{N}_i{K}_i{N}_i^T{T}_{i2}=\frac{6}{d}tr\left(K_i\right)A_i\left(I_i-Q_i{Q}_i^T\right)A_i& \left(12\right)\end{array}$

(42) Wherein X.sub.i=(x.sub.i1x.sub.i, L, x.sub.ix.sub.i, L, x.sub.in.sub.ix.sub.i).sup.T x.sub.i is the vector of the center of grid i, x.sub.i is the position vector of .sub.i center, N.sub.i=(|.sub.i1|n.sub.i1, L, |.sub.i|n.sub.i, L, |.sub.in.sub.i|n.sub.in.sub.i).sup.T, d is the grid dimension, A.sub.i=diag(|.sub.i1|, L, |.sub.i|, L, |.sub.in.sub.i|)) Q.sub.i=orth(A.sub.iX.sub.i)

(43) 2.3 Flux Continuity Conditions

(44) It can be seen from 2.1 and 2.2 that grid i and grid j have an outward normal flux F.sub.i and F.sub.j on the intersection e.sub.ij, respectively, the flux continuity condition should be satisfied, that is, the algebraic sum of F.sub.i and F.sub.j should be equal to 0.

(45) $\begin{array}{cc}{F}_i+F_j=\underset{=1}{\overset{n_i}{.Math.}}{T}_{,}^i\left(p_i^{t+t}-p_i^{t+t}\right)+\underset{=1}{\overset{n_j}{.Math.}}{T}_{,}^j\left(p_j^{t+t}-p_j^{t+t}\right)=0& \left(13\right)\end{array}$

(46) In fact, when using TPFA, p.sub.i can be eliminated according to Eqs. (9) and (13), the numerical flux based on TPFA can be calculated as Eq. (14), that is, the transmissibility formula based on the harmonic average scheme in the common TPFA, in essence, there is no need to add additional pressure degrees of freedom between grid i and grid j.

(47) $\begin{array}{cc}{F}_{ij}=\frac{k_{\mathrm{ra},\mathrm{ij}}}{{}_{a,ij}{B}_{a,ij}}{T}_{\mathrm{ij}}\left(p_i-p_j\right)T_{\mathrm{ij}}={\left(T_i^{-1}+T_j^{-1}\right)}^{-1}& \left(14\right)\end{array}$

(48) It can be considered that TPFA is a form of MFD in the case of K-orthogonality without adding additional pressure degrees of freedom to the connection between grids, however, in the case of non-K orthogonal, the definition of additional pressure degree of freedom is difficult to avoid. With the additional pressure degree of freedom, there will be a flux continuity condition, so that the global equations can still be solved in a closed form.

(49) 2.4 Hybrid of TPFA and MFD

(50) TPFA can be used to estimate the numerical flux on K-orthogonal grids, for example, it has a regular grid of isotropic permeability, while MFD-based flux approximation is used only in the remaining grids that do not satisfy K-orthogonality. According to Eq. (5), Eq. (6), Eq. (10) and Eq. (12), the discrete scheme of Eq. (3) can be obtained as follows:

(51) $\begin{array}{cc}\frac{k_{\mathrm{ra},\mathrm{ij}}}{{}_{a,ij}{B}_{a,ij}}\underset{=1}{\overset{n_i}{.Math.}}\underset{=1}{\overset{n_i}{.Math.}}\left[T_{,}^i\left(p_i^{t+t}-p_i^{t+t}\right)\right]+Q_a=\frac{.Math.{}_i.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}^{t+t}-{\left(\frac{S_a}{B_a}\right)}^t\right]& \left(15\right)\end{array}$

(52) The Eq. (13) and Eq. (15) constitute the global equations, and the Newton-Raphson (NR) iteration method is used to solve the equations to obtain the pressure, saturation distribution and production dynamic data.

(53) (3) EDFM and Generic pEDFM Workflow

(54) 3.1 EDFM

(55) EDFM mainly includes three types of connections, namely, the connections between two adjacent matrix grids, the connections between the matrix grid and the inter-grid fracture it contains, and the connection between inter-grid fractures. In order to facilitate the description, the connections between two grids are denoted by (.Math., .Math.). As shown in FIG. 1(a), the first type of connection includes (m.sub.i, m.sub.j) and (m.sub.i, m.sub.k), the second type of connection includes (m.sub.i, f.sub.i), (m.sub.j, f.sub.j) and (m.sub.k, f.sub.k), and the third type of connection includes (f.sub.i, f.sub.k).

(56) For the first type of connection, taking (m.sub.i, m.sub.j) as an example, when TPFA is used, similar to Eq. (14), it can be calculated that the transmissibility of (m.sub.i, m.sub.j) is half of the harmonic mean of the transmissibility of m.sub.i and the transmissibility of m.sub.j:

(57) $\begin{array}{cc}{T}_{m_i,m_j}={\left(T_{m_i}^{-1}+T_{m_j}^{-1}\right)}^{-1}{T}_{m_i}^{}=\frac{k_{m_i}A\left(1_{\mathrm{ij}}\right)}{d_{m_i}}{T}_{m_j}^{}=\frac{k_{m_j}A\left(1_{\mathrm{ij}}\right)}{d_{m_j}}& \left(16\right)\end{array}$ wherein A(.Math.) denotes the area operator and 1.sub.ij is the interface of m.sub.i and m.sub.j.

(58) Similarly, for the second type of connection,

(59) $\begin{array}{cc}{T}_{f_i,m_i}={\left(T_{f_i}^{-1}+T_{m_i}^{-1}\right)}^{-1}{T}_{f_i}=\frac{k_{f_i}{A}_{f_i}}{w_{f_i}/2}{T}_{m_i}=\frac{2k_{m_i}{A}_{f_i}}{.Math.d_{f_i{m}_i}.Math.}.Math.d_{f_i{m}_i}.Math.=\frac{{}_{{}_{m_i}}.Math.x.Math.n.Math.\mathrm{dV}}{V_{m_i}}& \left(17\right)\end{array}$

(60) For the third type of connection, it is generally necessary to subdivide it into two cases. The first case is the connection between inter-grid fractures on the same fracture surface. For example (f.sub.i, f.sub.k), the transmissibility is calculated as:

(61) $\begin{array}{cc}{T}_{f_i{f}_k}={\left(T_{f_i}^{-1}+T_{f_k}^{-1}\right)}^{-1},T_{f_i}=\frac{k_{f_i}{w}_{f_i}{h}_{f_i}}{d_{f_i}},T_{f_k}=\frac{k_{f_k}{w}_{f_k}{h}_{f_k}}{d_{f_k}}& \left(18\right)\end{array}$

(62) The second case is the connection between multiple fracture grids that may exist in the same matrix grid, which is generally calculated by star-delta transformation.

(63) 3.2 Generic pEDFM Process

(64) pEDFM is a model between EDFM and DFM, as shown in FIG. 1(b), taking f; as an example, pEDFM projects f.sub.i along the x and y directions to 1.sub.ij and 1.sub.ik, respectively, at this time, pEDFM transforms the original model in FIG. 1(a) into the model in FIG. 1(c). At this point, compared with EDFM, f.sub.i needs to be connected to m.sub.j and m.sub.k, the corresponding flow area is the projection area of f.sub.i meanwhile, since the connection of (m.sub.j, f.sub.i) and (m.sub.k, f.sub.i) occupies the area of 1.sub.ij and 1.sub.ik, the flow area of the original (m.sub.j, f.sub.i) and (m.sub.k, f.sub.i) is weakened. Meanwhile, because f.sub.i and f.sub.j are projected onto 1.sub.ij, there is also a connection between f.sub.i and f.sub.j.

(65) The core idea of the generic pEDFM workflow is to use the micro-translation projection method to deal with high-conductivity fracture units, and to use the non-projection transmissibility multiplier method to model low-conductivity fractures.

(66) When the matrix permeability is isotropic, the permeability of the m-f connections added to the pEDFM in FIG. 1 is:

(67) $\begin{array}{cc}{T}_{f_i{m}_i}={\left(T_{f_i}^{-1}+T_{m_j}^{-1}\right)}^{-1},T_{f_i}=\frac{k_{f_j}A\left(f_{j,1_{\mathrm{ij}}}^p\right)}{w_{f_j}/2},T_{m_j}=\frac{k_{m_j}A\left(f_{j,1_{\mathrm{ij}}}^p\right)}{d_{m_j,1_{\mathrm{ij}}}}& \left(19\right)\end{array}$

(68) $\begin{array}{cc}{T}_{f_j{m}_i}={\left(T_{f_j}^{-1}+T_{m_i}^{-1}\right)}^{-1},T_{f_j}=\frac{k_{f_j}A\left(f_{j,1_{\mathrm{ij}}}^p\right)}{w_{f_j}/2},T_{m_j}=\frac{k_{m_j}A\left(f_{j,1_{\mathrm{ij}}}^p\right)}{d_{m_j,1_{\mathrm{ij}}}}& \left(20\right)\end{array}$ wherein f.sub.i,1.sub.ij.sup.p and f.sub.j,1.sub.ij.sup.p are the projections of f.sub.i and f.sub.j on 1.sub.ij.

(69) The transmissibility of the original m-m connections is reduced to:

(70) 0$T_{m_i,m_j}={\left(T_{m_i}^{-1}+T_{m_j}^{-1}\right)}^{-1},T_{m_i}^{}=\frac{k_{m_i}\left(A\left(1_{\mathrm{ij}}\right)-A\left(f_{i,1_{\mathrm{ij}}}^p\right)\right)}{d_{m_i,1_{\mathrm{ij}}}},T_{m_j}^{}=\frac{k_{m_j}\left(A\left(1_{\mathrm{ij}}\right)-A\left(f_{j,1_{\mathrm{ij}}}^p\right)\right)}{d_{m_j,1_{\mathrm{ij}}}}$

(71) $T_{m_i,m_k}={\left(T_{m_i}^{-1}+T_{m_k}^{-1}\right)}^{-1},T_{m_i}^{}=\frac{k_{m_i}\left(A\left(1_{\mathrm{ij}}\right)-A\left(f_{i,1_{\mathrm{ik}}}^p\right)\right)}{d_{m_i,1_{\mathrm{ik}}}},T_{m_k}^{}=\frac{k_{m_k}\left(A\left(1_{\mathrm{ij}}\right)-A\left(f_{j,1_{\mathrm{ik}}}^p\right)\right)}{d_{m_k,1_{\mathrm{ik}}}}$

(72) By halving the transmissibility of the matrix grid, the transmissibility of the original m-f connections is weakened. Taking (m.sub.i, f.sub.i) as an example, the transmissibility is calculated as:

(73) $\begin{array}{cc}{T}_{f_i,m_i}={\left(T_{f_i}^{-1}+T_{m_i}^{-1}\right)}^{-1},T_{f_i}=\frac{k_{f_i}{A}_{f_i}}{w_{f_i}/2},T_{m_i}=\frac{k_{m_i}{A}_{f_i}}{.Math.d_{f_i{m}_i}.Math.},.Math.d_{f_i{m}_i}.Math.=\frac{{}_{{}_{m_i}}.Math.x.Math.n.Math.\mathrm{dV}}{V_{m_i}}& \left(23\right)\end{array}$

(74) Added transmissibility (f.sub.i, f.sub.i):

(75) $\begin{array}{cc}{T}_{f_i{f}_j}={\left(T_{f_i}^{-1}+T_{f_j}^{-1}\right)}^{-1},T_{f_i}=\frac{k_{f_i}{A}_{f_i}^p}{w_{f_i}/2},T_{f_j}=\frac{k_{f_j}{A}_{f_j}^p}{w_{f_j}/2}& \left(24\right)\end{array}$

(76) Then, the generic pEDFM gives a non-projection method to model low-conductivity fractures, as shown below: the updated transmissibility of the connection is:

(77) $\begin{array}{cc}{T}_{\mathrm{ij}}^{}={\left(T_{\mathrm{ij}}^{-1}+T_f^{-1}\right)}^{-1},T_f=\frac{k_f{A}_{\mathrm{ij}}}{w_f}& \left(25\right)\end{array}$

(78) Wherein T.sub.ij and T.sub.i are the updated and original transmissibility of the connection; connected updated and original transmissive semi-transmissibility of low-conductivity fractures; k.sub.f and w.sub.f are the permeability and pore size of low-conductivity fractures, respectively; A.sub.ij is the flow cross-sectional area corresponding to the connection; if k.sub.f=0, the updated transmissibility in Eq. (25) will become zero, so the flow barrier can block the flow. It should be noted that if multiple low-conductivity fractures intersect with the same connection, the transportability of the connection will be updated multiple times by using Eq. (25).

(79) 4.1 the Basic Idea of the Present Invention

(80) The connection between the matrix and the fracture in EDFM is very simple, and only the connection between the matrix grid and the fracture grid contained in the matrix grid is established, it is essentially a local dual-medium model, which makes the m-f connections in EDFM not affect the original m-m connections between adjacent matrix inter-grid. Therefore, when using MFD to deal with the full tensor permeability of the matrix grid in EDFM, it is only necessary to use the transfer matrix based on dense MFD in the structured background matrix grid to replace the diagonal transmissibility matrix obtained by TPFA in the original matrix grid. Of course, it can be further considered that the permeability of the matrix grid is the influence of the full tensor case on the m-f transmissibility in EDFM. In general, due to the simplicity of the grid connection structure in EDFM, it is a simple task to construct an EDFM based on MFD. However, pEDFM, which performs better than EDFM, has a much more complex connection relationship than EDFM. The most important point is that the m-f connections in pEDFM will weaken the original m-m connections, or even directly cover the original m-m connections and make them disappear, which makes the above strategies that can work in EDFM not directly work in pEDFM. The specific explanation is as follows:

(81) As shown in FIG. 2(a), m, contains a fracture unit f.sub.i, m.sub.i and m.sub.j are adjacent, in EDFM, the pressure is continuous on the intersection I.sub.ij of m.sub.i and m.sub.j, this facilitates the application of the aforementioned MFD theory to construct flux continuity conditions at I.sub.ij. However, because the essence of pEDFM is to project fractures onto the interface of matrix grids, the original embedded discrete scheme is transformed into an approximate DFM to deal with. In pEDFM, as shown in FIG. 2(b), the f.sub.i in the matrix grid m.sub.i needs to be projected on I.sub.ij, but not fully covered I.sub.ij, the area covered by the fracture projection is denoted as I.sub.ij.sup.1, the remaining region on the I.sub.ij plane is I.sub.ij.sup.2, the pressure is continuous only on the m.sub.i-side and m.sub.j-side of I.sub.ij.sup.2, the pressure on the m.sub.i-side and m.sub.i-side of I.sub.ij.sup.1 may be discontinuous, for example, if the permeability of the fracture is only slightly higher than that of the matrix grid, at this time, fluid is displaced to m.sub.j, from m.sub.i, and there is a pressure difference that cannot be ignored on both sides of I.sub.ij.sup.1. Therefore, the pressure values on different connections may be distributed on the intersection of the matrix grids that undertake the fracture projection in pEDFM, and it is necessary to distinguish the possible discontinuous pressures on the left and right sides on the intersection of the matrix grids.

(82) The above analysis gives the following important inspirations:

(83) i) The concept of effective connectivity needs to be introduced. As shown in FIG. 2(c), the projection of fracture f.sub.j on I.sub.ij will completely cover I.sub.ij, at this time, there will be no need to define additional pressure degrees of freedom between m.sub.i and m.sub.j. If the additional definition is given, the global equations will not have the flow continuity equation related to the pressure degree of freedom, so that the global equations can not be solved in closed form. Therefore, it is necessary to introduce the definition of effective connection, that is, the connection whose actual flow area is not 0. Only such a connection can define an additional pressure degree of freedom.

(84) ii) The pressure degree of freedom is only added to the effective connection related to the matrix grid. Since TPFA can still be used to calculate the numerical flux in the fracture. Therefore, the transmissibility information of the f-f connections in the generic pEDFM workflow can be kept unchanged without the need to define additional pressure degrees of freedom on the original f-f connections to carry out MFD. In other words, by obtaining effective m-m connections and m-f connections, an additional pressure degree of freedom is defined only on these effective connections related to the matrix grid.

(85) iii) The distribution of the pressure degree of freedom on which side of the matrix grid surface should be solved. As mentioned above, the essence of pEDFM is to project the fracture onto the interface of the matrix grid, and transform the original embedded discrete scheme into an approximate DFM to deal with. Therefore, the added pressure degrees of freedom can be considered to be distributed on the matrix grid surface. Meanwhile, because the pressure on both sides of the fracture projection area on the matrix grid surface may no longer be continuous, the pressure degrees of freedom added by different connections need to be distinguished to which side of the matrix grid surface is located. As shown in FIG. 2(b), the projection area of fracture f.sub.j on I.sub.ij is not 0, therefore, there are two effective connections of (m.sub.i, f.sub.i) and (m.sub.j, f.sub.i). Then the pressure degree of freedom added for (m.sub.i, f.sub.i) is actually located on the side of I.sub.ij near m, side, the pressure degree of freedom added for (m.sub.j, f.sub.i) is actually located on the side of I.sub.ij near m.sub.j side, and the values of these two pressure degrees of freedom are generally different. Therefore, for an effective m-f connection, if the fracture f is projected on the interface 1 of the matrix grid, the additional degree of freedom added for the connection is judged to be located on the m side of 1. For an effective (m.sub.i, m.sub.j) connection, the pressure degree of freedom is located on both sides of I.sub.ij. In general, it is based on the matrix grid in the connection to solve the pressure degree of freedom added to the connection on which side of the corresponding matrix grid surface is located.

(86) (iv) The concept of average pressure on the matrix grid surface needs to be introduced. Since the pressure on the grid surface is needed to obtain the MFD flow operator in the matrix grid, pEDFM is different from EDFM as mentioned above. It may have multiple connections on the matrix grid surface, not only the m-m connections in EDFM, but also the original m-m connections may be completely replaced by other m-f connections. At this time, if the value of the pressure degree of freedom added to the different connections is subdivided on each grid surface to participate in the construction of the transmissibility matrix based on MFD in the matrix grid, the complexity of the algorithm will be significantly increased and the robustness and practicability of the algorithm will be reduced. Therefore, when constructing the transmissibility matrix based on MFD in the matrix grid, the average pressure on the side of the matrix grid surface close to the matrix grid surface should be used (it should be noted that, as mentioned above, the matrix grid surface may have different average pressures on two sides of the matrix grid surface). The average pressure is a weighted average of the additional pressure degrees of freedom added by each effective connection on the side of the matrix grid surface near the matrix grid with the corresponding flow area as the weight. As shown in FIG. 2(b), the average pressure p.sub.1.sub.ij.sup.m.sup.i of the side near m.sub.i of 1.sub.ij and the average pressure p.sub.1.sub.ij.sup.m.sup.i of the side near m, of 1 are calculated as follows:

(87) $\begin{array}{cc}{\stackrel{}{p}}_{1_{\mathrm{ij}}}^{m_i}=\frac{p_{\left(m_i,f_i\right)}{A}_{\left(m_i,f_i\right)}+p_{\left(m_i,m_j\right)}{A}_{\left(m_i,m_j\right)}}{A_{\left(m_i,f_i\right)}+A_{\left(m_i,m_j\right)}},{\stackrel{}{p}}_{1_{\mathrm{ij}}}^{m_j}=\frac{p_{\left(m_i,f_i\right)}{A}_{\left(m_i,f_i\right)}+p_{\left(m_i,m_j\right)}{A}_{\left(m_i,m_j\right)}}{A_{\left(m_i,f_i\right)}+A_{\left(m_i,m_j\right)}}& \left(26\right)\end{array}$ wherein, p.sub.(m.sub.i.sub., f.sub.i.sub.), p.sub.(m.sub.j.sub., f.sub.i.sub.) and p.sub.(m.sub.i.sub., m.sub.j.sub.), and are the pressure degrees of freedom added for (m.sub.i, f.sub.i), (m.sub.j, f.sub.i) and (m.sub.i, m.sub.j), respectively, A.sub.(m.sub.i.sub., f.sub.i.sub.), A.sub.(m.sub.j.sub., f.sub.i.sub.) and A.sub.(m.sub.i.sub., m.sub.j.sub.) are the corresponding flow areas of (m.sub.i, f.sub.i), (m.sub.j, f.sub.i) and (m.sub.i, m.sub.j), respectively. It can be seen that the pressure on the left and right sides of 1.sub.ij.sup.1. is discontinuous, the average pressure on the left and right sides of 1.sub.ij is no longer continuous. Of course, if there is no fracture projection on the matrix grid surface, the average pressure on both sides of 1.sub.ij can be calculated to be continuous, which degenerates into the case in EDFM.
4.2 Basic Concepts
Concept 1: effective connection

(88) Since the generic pEDFM workflow adopts the non-projection processing method of transmissibility multiplier for low-conductivity fractures, this section first ignores low-conductivity fractures, and at this time, the connection set between grids Cont.sub.1 can be obtained, the transmissibility of the corresponding connection is calculated according to Eq. (19) to Eq. (25), and the corresponding flow area is also reflected in these transmissibility calculation formulas. In Cont.sub.1, the set of connections with flow area of 0 is Cont.sub.1.sup.0, it can be seen that the transmissibility of the connection in Cont.sub.1.sup.0 must be 0, and it has no effect on the results of simulation calculation. Taking FIG. 2 as an example, C(m.sub.i, m.sub.j)Cont.sub.1.sup.0, therefore, the effective set in definition Cont.sub.1 is Cont.sub.1.sup.eff.
Cont.sub.1.sup.eff=Cont.sub.1Cont.sub.1.sup.0(27)

(89) Concept 2: Cont.sub.1.sup.eff (1.sub.i.sub.j.sup.m.sup.i) and the corresponding pressure degree of freedom

(90) It can be seen from the basic theory of MFD in 2.2, since the permeability of fracture grip is generally isotropic, when pEDFM is constructed based on MFD, there is no need to change the original f-f connections and the corresponding transmissibility in Cont.sub.1.sup.eff. However, the permeability of the matrix grid may be in the full tensor form, so it is necessary to construct an additional pressure degree of freedom for each m-m and f-m connection in Cont.sub.1.sup.eff. As mentioned above, pEDFM essentially projects the fracture grid in the matrix grid to the intersection surface of the matrix grid to form an approximate DFM to deal with. It can be considered that each m-m and f-m effective connection increases. A pressure degree of freedom is located on the matrix grid surface associated with the connection. For example, (m.sub.i, f.sub.i) in FIG. 2(b), The additional pressure degree of freedom p.sub.(m.sub.i.sub., f.sub.i.sub.) is considered to be on the side of 1.sub.ij near m.sub.i. Let 1.sub.ij.sup.m.sup.i be the side of 1.sub.ij near m.sub.i, 1.sub.ij.sup.m.sup.j is the side of 1.sub.ij near m.sub.j. Then, it is possible to define the effective connection set associated with 1.sub.ij.sup.m.sup.i as a set of connections in which the added pressure degree of freedom in Cont.sub.1.sup.eff is located in 1.sub.ij.sup.m.sup.i, noting Cont.sub.1.sup.eff (1.sub.ij.sup.m.sup.i), the flow area of the corresponding connection and the increased pressure degree of freedom are A.sub. and p.sub., respectively, Cont.sub.1.sup.eff (1.sub.ij.sup.m.sup.i). It should be pointed out that the increased pressure degree of freedom of the original m-f connection may fall on two substrate grids, such as f.sub.i in FIG. 2(c), it is due to the projection in the x direction and they direction, respectively, therefore, (f.sub.i, m.sub.i)Cont.sub.1.sup.eff (1.sub.ij.sup.m.sup.i) and (f.sub.i, m.sub.i)Cont.sub.1.sup.eff (1.sub.ik.sup.m.sup.i), the corresponding flow area should not be fracture area, but A(f.sub.i,1.sub.ij.sup.p) and A(f.sub.i,1.sub.ik.sup.p), respectively.

(91) Concept 3: Average Pressure of Matrix Grid Surface

(92) From the MFD theory of 2.2, it can be seen that in the case of full permeability tensor, the outward normal flux on one side of the grid is also related to the pressure on other sides of the grid. Therefore, the average pressure p.sub.ij.sup.i on one side of the matrix grid is defined to participate in the calculation of the outward normal flux on the side of the matrix grid based on MFD and the calculation is as follows:

(93) $\begin{array}{cc}{\stackrel{}{p}}_{1_{\mathrm{ij}}}^{m_i}=\frac{\underset{{\mathrm{Cont}}_1^{\mathrm{eff}}\left(1_{\mathrm{ij}}^{\mathrm{mj}}\right)}{.Math.}{p}_{}{A}_{}}{\underset{{\mathrm{Cont}}_1^{\mathrm{eff}}\left(1_{\mathrm{ij}}^{\mathrm{mj}}\right)}{.Math.}{A}_{}}& \left(28\right)\end{array}$
4.3 Treatment of Low-Conductivity Fractures

(94) As described in 3.2, the generic pEDFM uses the transmissibility multiplier method to deal with low-conductivity fractures, so that pEDFM can generally model various low-conductivity fractures. The transmissibility multipliers are based on the harmonic average scheme in TPFA, and the transmissibility of each connection affected by the low-conductivity fractures is updated with the transmissibility of the low-conductivity fractures. Let f.sub.l be a low-conductivity fracture, and its permeability and fracture width are k.sub.f and w.sub.f, respectively. Considering that the fracture grid in the pEDFM of the invention still adopts the transmissibility based on TPFA, the connection affected by low-conductivity fractures is divided into the following three treatment methods according to different connection types: i) for f-f connections, it may be set to f.sub.i and f.sub.j, and the f.sub.i-f.sub.j transmissibility before being affected by f.sub.l is:
T.sub.f.sub.i.sub.f.sub.j=(T.sub.f.sub.i.sup.1+T.sub.f.sub.j.sup.1).sup.1(29) then the transmissibility affected by f.sub.l is:

(95) $\begin{array}{cc}{T}_{f_i{f}_j}=\left(T_{f_i{f}_j}^{-1}+T_{f_l}^{-1}\right),T_{f_l}=\frac{k_{f_l}{A}_{\left(f_i,f_j\right)}}{w_{f_l}/2}& \left(30\right)\end{array}$ wherein, A.sub.(f.sub.i.sub., f.sub.j.sub.) is the flow area of (f.sub.i, f.sub.j). if there are other low-conductivity fractures that affect the connection, then the Eq. (29) continues to be updated. ii) For m-f connections, it may be set to m.sub.i and f.sub.j, suppose that a transmissibility off in the m.sub.i-f.sub.j connections is not affected by f.sub.l, the original transmissibility of f.sub.j is based on the TPFA's T.sub.f.sub.j, the influence of f.sub.l is applied to f.sub.j, which is equivalent to connecting the low-conductivity fracture f.sub.l with f.sub.j in series, that is, the transmissibility off is updated as:
T.sub.f.sub.j=(T.sub.f.sub.j.sup.1+T.sub.f.sub.l.sup.1).sup.1(31) if there are other low-conductivity fractures that affect the m.sub.i-f.sub.j connections, then use Eq. (31) to continue to update T.sub.f.sub.j. iii) for m-m connections, it may be set to m.sub.i and f.sub.j, since there is no fracture grid available for the harmonic averaging scheme, it is necessary to split the m.sub.i-m.sub.j connections into two m-f connections, they are m.sub.i-f.sub.l and m.sub.j-f.sub.l respectively, and the flow area of these two connections is (m.sub.i, m.sub.j) flow area A.sub.(m.sub.i.sub., m.sub.j.sub.), wherein the transmissibility off is:

(96) $\begin{array}{cc}{T}_{f_l}=\frac{k_{f_l}{A}_{\left(m_i,m_j\right)}}{w_{f_l}/2}& \left(32\right)\end{array}$

(97) If there are other low-conductivity fractures (such as f.sub.k), then the transmissibility off is updated with the transmissibility of f.sub.k by Eq. (31), that is:

(98) $\begin{array}{cc}{T}_{f_l}={\left(T_{f_l}^{-1}+T_{f_k}^{-1}\right)}^{-1},T_{f_k}=\frac{k_{f_k}{A}_{\left(m_i,m_j\right)}}{w_{f_k}/2}& \left(33\right)\end{array}$

(99) It can be seen that the treatment of low-conductivity fractures in this section will add new connections in Cont.sub.1.sup.eff, and the updated Cont.sub.1.sup.eff is recorded as Cont.sub.2.sup.eff.

(100) 4.4 Numerical Flux Calculation of Effective Connection

(101) Let the number set of the matrix grid which is adjacent to the matrix grid m.sub.i be neigh.sub.m.sub.i, then the set of effective connections on each side of m.sub.i near m.sub.i is:

(102) 0${\mathrm{Cont}}_2^{\mathrm{eff}}\left(m_i\right)=\underset{j{\mathrm{neigh}}_{m_i}}{U}{\mathrm{Cont}}_2^{\mathrm{eff}}\left(1_{\mathrm{ij}}^{m_i}\right),$
according to Eq. (28), the conversion formula between the average pressure of each side of m.sub.i near m.sub.i and the value of pressure degree of freedom added by each connection in Cont.sub.1.sup.eff (m.sub.i) can be obtained as follows:
p.sub.m.sub.i=Ap.sub.m.sub.i(34) wherein p.sub.m.sub.i is a column vector composed of p.sub.ij.sup.m.sup.i (jneigh.sub.m.sub.i), p.sub.m.sub.i is a column vector composed of p.sub.j(jCont.sub.2.sup.eff (m.sub.i)), each j row of A represents the area weight when the average pressure of the jth surface of m.sub.i near m.sub.i is calculated by Eq. (28), meanwhile, A.sup.T is the distribution matrix of the numerical flux on each surface of m.sub.i to the flow on each connection in Cont.sub.1.sup.eff(m.sub.i):
flux.sub.i=A.sup.T{square root over (flux.sub.i)}(35)

(103) Firstly, it is judged whether the matrix grid m.sub.i is K orthogonal, if yes, then the transmissibility matrix T.sup.m.sup.i in m.sub.i is calculated by using the Eq. (9) based on TPFA, if not, then the transmissibility matrix T.sup.m.sup.i in m.sub.i is calculated by using the MFD-based Eq. (12). As mentioned above, the average pressure is used to participate in the calculation of the numerical flow on each surface of the matrix grid, so it can be obtained that:
flux.sub.i=T.sup.m.sup.ip.sub.m.sub.i=T.sup.m.sup.iAp.sub.m.sub.i=T.sup.m.sup.iA(p.sub.m.sub.iIp.sub.m.sub.i)(36) wherein, I is a column vector, wherein one of the elements is 1 and the length is the number of effective connections of the matrix grid.

(104) Furthermore, combining Eq. (35) and Eq. (36), it is obtained that:
flux.sub.i=A.sup.TT.sup.m.sup.iA(p.sub.m.sub.iIp.sub.m.sub.i)(37)

(105) Therefore, for the matrix grid, the actual transmissibility matrix is:
=A.sup.TT.sup.m.sup.iA(38)
4.4 Global Equation

(106) By using the transmissibility matrix given in Eq. (38), the flux on each connection related to m.sub.i can be calculated, by taking into Eq. (15), the discrete scheme of the continuity equation in m; can be obtained:

(107) $\begin{array}{cc}\underset{{\mathrm{Cont}}_1^{\mathrm{eff}}\left(m_i\right)}{.Math.}\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}\underset{{\mathrm{Cont}}_1^{\mathrm{eff}}\left(m_i\right)}{.Math.}\left(p_{}^{t+t}-p_{m_i}^{t+t}\right)+Q_a=\frac{.Math.{}_i.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}_{m_i}^{t+t}-{\left(\frac{S_a}{B_a}\right)}_{m_i}^t\right]& \left(39\right)\end{array}$

(108) For a fracture grid f.sub.j, its transmissibility still adopts a simple TPFA scheme, the efficient connections Cont.sub.2.sup.eff (f.sub.j) related to f.sub.j include f-f connections (Cont.sub.fj.sup.f-f) and m-f connections (Cont.sub.fj.sup.m-f), the set of fracture grids adjacent to f.sub.j reflected from Cont.sub.fj.sup.f-f is denoted as neigh.sub.f.sub.j, then the discrete scheme of the continuity equation in f.sub.j can be written as Eq. (40). It should be emphasized that the reason why the form of the mass transfer part corresponding to the f-f connections in Eq. (40) is different from the form of the mass transfer part corresponding to the m-f connections is that: in this work, in order to reduce the computational cost and algorithm complexity, no additional pressure degree of freedom is added to the f-f connections. Therefore, the mass transfer corresponding to the f-f connections is directly expressed in a numerical flux expression similar to that in Eq. (14) using grid average pressure and TPFA.

(109) $\begin{array}{cc}\underset{Con{t}_{f_j}^{m-f}}{.Math.},\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}\underset{Con{t}_{f_j}^{m-f}}{.Math.}\left(p_{}^{t+t}-p_{f_j}^{t+t}\right)+\underset{neig{h}_{f_j}}{.Math.}\frac{k_{\mathrm{ra},i}}{{}_{a,i}{B}_{a,i}}{T}_{f_j}\left(p_{}^{t+t}-p_{f_j}^{t+t}\right)+Q_a=\frac{.Math.{}_{f_j}.Math.}{t}\left[{\left(\frac{S_a}{B_a}\right)}_{f_j}^{t+t}-{\left(\frac{S_a}{B_a}\right)}_{f_j}^t\right]& \left(40\right)\end{array}$

(110) Meanwhile, when of each matrix grid is obtained, the transmissibility of all effective m-m connections and m-f connections matrix grids is known, and the transmissibility of the fracture grid in the m-f connections is still based on the TPFA scheme. Therefore, the continuity conditions of each m-m connection and m-f connection in Cont.sub.2.sup.eff can be given.

(111) For m-m connections, suppose that is m.sub.i and m.sub.k, then:
flux.sub.m.sub.i.sub..fwdarw.m.sub.k=(p.sub.m.sub.iIp.sub.m.sub.i),flux.sub.m.sub.k.sub..fwdarw.m.sub.i=(p.sub.m.sub.iIp.sub.m.sub.i)(41)

(112) The corresponding flux continuity condition is: (p.sub.m.sub.iIp.sub.m.sub.i)+(p.sub.m.sub.iIp.sub.m.sub.i)=0 (42)

(113) For m-f connections, suppose that is m.sub.i and f.sub.j, then:

(114) $\begin{array}{cc}{\mathrm{flux}}_{m_i.fwdarw.f_j}=\left(p_{m_i}I-p_{m_i}\right),{\mathrm{flux}}_{f_j.fwdarw.m_i}=T_{f_j}\left(p_{f_j}-p_{\left(m_i,f_j\right)}\right)& \left(43\right)\end{array}$

(115) The corresponding flux continuity condition is:
(p.sub.m.sub.iIp.sub.m.sub.i)+T.sub.f.sub.j(p.sub.f.sub.jp.sub.(m.sub.i.sub.,f.sub.j.sub.))=0(44)

(116) For f-f connections, because the present invention still uses the transmissibility formula in Eq. (18) and Eq. (24) without defining additional pressure degrees of freedom in the generic pEDFM, there is no additional flow continuity condition for f-f connections.

(117) In general, let the reservoir calculation domain contain n.sub.m matrix grids and n.sub.f fracture grids, the number of m-m and m-f connections contained in Cont.sub.2.sup.eff is n.sub.c. The pressure degrees of freedom of the new pEDFM include n.sub.m matrix grid center pressure, n.sub.f fracture grid average pressure, and n.sub.c additional pressure degrees of freedom, a total of n.sub.m+n.sub.f+n.sub.c, the degree of freedom of water saturation includes n.sub.m matrix grid average saturation and n.sub.f fracture grid average saturation, a total of n.sub.m+n.sub.f, therefore, the global degree of freedom is 2(n.sub.m+n.sub.f)+n.sub.c. The global equations include a continuity equation of 2n.sub.m matrix grids (including oil phase and water phase, and therefore 2n.sub.m, not n.sub.m), that is also the Eq. (39), the continuity equation of n.sub.f fracture grids (including oil phase and water phase), namely Eq. (40), and n.sub.c flux continuity conditions composed of Eq. (42) and Eq. (44), Therefore, the global equations are also 2(n.sub.m+n.sub.f)+n.sub.c, and it can be closed solution. The nonlinear solver based on Newton-Raphson (NR) iteration is used to solve the global equations to obtain the pressure and water saturation distribution at each time.

Example 1

(118) FIGS. 3(a) and (b) show two reservoir models in this case, in which the first model contains only one high-conductivity fracture consistent with the coordinate line, and the second model has an additional flow barrier. There is a water injection well and a production well in the lower left corner and the upper right corner of the reservoir respectively. FIGS. 3(c) and (d) show the Cartesian background grid of the two models. The reservoir permeability is the anisotropic permeability tensor of the principal axis in the x and y directions in Eq. (45). The relative permeability function is shown in Eq. (46), and the other physical parameters are shown in table 1.

(119) $\begin{array}{cc}K=\left[\begin{array}{cc}80& 0\\ 0& 10\end{array}\right]mD& \left(45\right)\end{array}$

(120) $\begin{array}{cc}{k}_{\mathrm{rw}}=0.8{\left(\frac{S_w-0.2}{1-0.2-0.2}\right)}^2,k_{ro}=0.552{\left(\frac{1-0.2-S_w}{1-0.2-0.2}\right)}^2& \left(46\right)\end{array}$
Table 1 Physical Parameters of Reservoir Model

(121) TABLE-US-00001 Physical parameters The value Physical parameters The value Original formation 15 MPa Water phase viscosity 0.6 mPa .Math. s pressure Initial water saturation 0.2 Reference pressure 15 MPa Thickness 5 m Oil phase compressibility 1 10.sup.3 MPa.sup.1 coefficient Matrix porosity 0.12 Water phase compressibility 5 10.sup.4 MPa.sup.1 coefficient Fracture porosity 0.4 Matrix compressibility 1 10.sup.4 MPa.sup.1 coefficient Fracture permeability 100000 mD Fracture compressibility 1 10.sup.4 MPa.sup.1 coefficient Fracture aperture 0.01 m Oil phase volume coefficient 1.0 under reference pressure Fracture permeability of 0 mD Water phase volume coefficient 1.0 blocking fracture under reference pressure Oil phase viscosity 2 mPa .Math. s

(122) At this time, the Cartesian grids are K-orthogonal, and the reference solution can be obtained by local grid refinement (LGR) and the finite volume method based on TPFA. FIG. 4 compares the water saturation distribution at 200 days calculated by LGR, new pEDFM, EDFM and DFM at case 1. FIG. 5 and FIG. 6 compare the water saturation and pressure distribution at 200 days calculated by LGR, new pEDFM, EDFM and DFM at case 2, respectively. FIG. 7 compares the water cut curves of production wells in case 1 and case 2. It can be seen that the new pEDFM can achieve almost the same saturation, pressure distribution and well response results as LGR and DFM. However, EDFM produces significant errors when simulating two-phase flow across high-conductivity fractures in case 1. The flow barrier in case 2 did not prevent the flow of injected water in EDFM, and the pressure distribution calculated by EDFM did not reflect the discontinuity of pressure on both sides of the flow barrier. The dynamic data of oil wells calculated by EDFM are also significantly different from the reference solution. The analysis of the results of this example shows that the new pEDFM can effectively solve the limitations of EDFM in the applicability of flow scenarios.

Example 2

(123) Keeping the reservoir model and physical parameters of case 2 in case 1 unchanged, only the permeability tensor in Eq. (46) is rotated by 30 degrees to obtain the full permeability tensor in Eq. (47). Theoretically, TPFA will not be able to obtain high-precision numerical flux approximation at this time. Taking the DFM solution based on MFD as the reference solution, FIG. 8 and FIG. 9 compare the water saturation and pressure distribution at 200 days calculated by pEDFM using TPFA, new pEDFM, and EDFM using MFD, respectively. FIG. 10 compares the water cut curves of production wells calculated by different methods. It can be seen that: (i) the new pEDFM can obtain almost the same calculation results as the reference solution, and the calculated water flooding front and pressure distribution accurately reflect the principal axis direction of the permeability tensor in Eq. (47) (i.e., the angle with the coordinate axis is 30 degrees). (ii) Although the pEDFM using TPFA can simulate the blocking effect of the flow barrier on the flow, the difference with the reference solution is still very obvious. From the calculated water drive front and pressure distribution, it can be seen that the use of TPFA to estimate the numerical flux leads to its failure to accurately grasp the principal axis direction of the permeability tensor in Eq. (47), and mistakenly believe that the principal axis direction of the permeability tensor is x direction and y direction. (iii) The EDFM using MFD fails to effectively characterize the role of the flow barrier, which is also reflected in case 2 of the first case, but the numerical flux approximation based on MFD enables EDFM to simulate a water flooding front with an inclination of about 30 degrees on the left side of the high-conductivity fracture.

(124) $\begin{array}{cc}K=\left[\begin{array}{cc}62.5& 30.3109\\ 30.3109& 27.5\end{array}\right]mD& \left(47\right)\end{array}$

(125) The above results demonstrate that the new pEDFM can achieve high-precision simulation results for the case with full permeability tensor, while the pEDFM using TPFA and the EDFM using MFD will have significant errors.

Example 3

(126) In this example, a more complex implementation example will be used to show the application effect of the new pEDFM, and to test that the new pEDFM can achieve the same calculation accuracy as DFM in the case of more complex slit network, and to show the advantages of the new pEDFM in grid generation compared with DFM. As shown in FIG. 11, the reservoir model of this implementation case contains 36 high-conductivity fractures marked with solid lines and 4 flow barriers marked with dotted lines. The permeability of anisotropic reservoir is the full permeability tensor in Eq. (47). The initial pressure and reference pressure of the reservoir are both 20 MPa. The other physical parameters are the same as those in table 1, and the relative permeability data are the same as that in Eq. (46). Wherein, the injection-production wells in case1 are located in the lower left corner and the upper right corner of the reservoir, and the injection-production wells in case2 are located in the upper left corner and the lower right corner of the reservoir. FIGS. 13(c) and (d) show the matching grid used by DFM and the non-matching grid used by the new pEDFM, respectively. It can be seen that when generating the matching grid for the complex fracture network, a large number of small grids will be generated in the narrow area between the fractures, so that the number of generated grids is large and the generation of grids is difficult. Practice shows that in the case of more complex fracture grids, even some mature triangular (tetrahedral) grid generation software cannot generate grids that match the geometric structure of the fracture network.

(127) FIG. 12 and FIG. 13 compare the water saturation and pressure distribution after 400 days and 600 days calculated by DFM and new pEDFM in case1 and case 2, respectively. It can be seen intuitively that the new pEDFM can achieve the same calculation accuracy as DFM in case1 and case 2. It is shown that the new pEDFM is more practical than the DFM for implementations with complex grids.

(128) Therefore, the present invention adopts a projection embedded DFM based on the hybrid method of TPFA and MFD, wherein the new pEDFM can deal with anisotropic reservoirs with full permeability tensor. For the first time, the implementation of TPFA-MFD (or MFD) is realized in the pEDFM framework, which significantly expands the original generic pEDFM using TPFA as a special case of the new pEDFM in the case of K-orthogonal grids, and significantly expands the application scope of the pEDFM framework.

(129) Finally, it should be noted that the above implementation examples are only used to illustrate the technical scheme of the invention rather than to restrict it, although the invention is described in detail with reference to the better implementation example. The ordinary technical personnel in this field should understand that they can still modify or replace the technical scheme of the invention, and these modifications or equivalent substitutions can not make the modified technical scheme out of the spirit and scope of the technical scheme of the invention.