A Numerical Modelling Method of Fractured Reservoirs with Embedded Meshes and Topological Fracture Projection Configurations

Xiang Raoand Yina Liu

1Cooperative Innovation Center of Unconventional Oil and Gas (Ministry of Education &Hubei Province),Yangtze University,Wuhan,China

2School of Petroleum Engineering,Yangtze University,Wuhan,China

ABSTRACT Projection-based embedded discrete fracture model(pEDFM)is an effective numerical model to handle the flow in fractured reservoirs,with high efficiency and strong generalization of flow models.However,this paper points out that pEDFM fails to handle flow barriers in most cases,and identifies the physical projection configuration of fractures is a key step in pEDFM.This paper presents and proves the equivalence theorem,which explains the geometric nature of physical projection configurations of fractures,that is,the projection configuration of a fracture being physical is equivalent to it being topologically homeomorphic to the fracture,by analyzing the essence of pEDFM.Physical projection configurations of fractures may be rigorously established based on this theorem,allowing pEDFM to obtain physical numerical results for many flow models,particularly those with flow barriers.Furthermore,a natural idea emerges of employing flow barriers to flexibly‘cut’the formation to quickly handle the flow problems in the formation with complex geological conditions,and several numerical examples are implemented to test this idea and application of the improved pEDFM.

KEYWORDS Numerical modeling method;fractured reservoirs;fracture-projection configuration;flow barriers;topologically homeomorphism

1 Introduction

Accurate and efficient numerical simulation of flow in fractured media is a challenging and important topic,which is of great interest in hydrocarbon recovery [1–4],aquifer management,and carbon dioxide geological sequestration,etc.The discrete fracture model (DFM),which explicitly handles large-scale fractures,has been a common technique to numerically model the flow in fractured media over the past two decades.Unstructured matrix meshes were generated to conform to the geometry of complex fracture network such that fractures lied at the interfaces between matrix cells.DFM thus accurately accounted for the effect of complex fracture geometries on flow,yielding physical and accurate simulation results.For single-phase flow,multiphase flow,and compositional flow models,A large variety of DFMs with different based numerical calculation methods had been presented in literature for single-phase flow,multiphase flow,and compositional flow models,including finite difference method (FDM) [5,6],Galerkin finite element method (FEM) [6,7],finite volume method (FVM) [7],control volume finite element method (cvFEM) [8,9],and mixed finite element method (mFEM) [10].However,a high-quality unstructured mesh to match complex fracture networks was challenging to generate [11],and often comprised of huge numbers of small-size cells near fractures.These resulted in high computing costs,limiting the field applications of DFM.

The above problems led to the emergence of embedded discrete fracture model (EDFM) which made use of non-conforming mesh.Li et al.[12] proposed EDFM,in which fractures were treated as source or sink terms in matrix cells,allowing for the use of a structured background mesh.Then,inter-cell connections were obtained by figuring out the geometric relationships between embedded discrete fractures and structured matrix cells.That is,in contrast to DFM,EDFM avoided the complexity of creating a high-quality conforming mesh.EDFM has been used to numerically simulate a variety of flow problems in fractured media [13–21].To increase computational performance in heterogeneous formations,Hajibeygi et al.[22] presented a hierarchical embedded fracture model utilizing iterative multiscale FVM.Moinfar et al.[23] developed an efficient EDFM for three-dimensional (3D) compositional reservoir simulation.Mimetic FDM and its multiscale technique-based EDFM were reported by Yan et al.[24] and Zhang et al.[25],respectively,in which velocity and phase pressure fields were solved simultaneously.Zeng et al.[26]presented a coupled model of EDFM and extended finite element method (XFEM) to study dynamic behaviors of hydraulic fractures.Rao et al.[27] presented a multi-layer virtual-cell EDFM to improve the early-time simulation accuracy for the production forecast of fractured horizontal wells.Wu et al.[28] proposed a Green element method-based EDFM that may yield high-accuracy results for single-phase flow.

Although EDFM had been validated by the above works and used in various porous flow problems [29,30],Tene et al.[31] found that classical EDFM could not effectively handle cases when fracture permeability was lower than matrix permeability,especially when flow barriers were present.To resolve the limitation,they developed a projection-based EDFM (pEDFM) by projecting the fracture cells to matrix-cell interfaces,constructing additional fracture-matrix (fm) connections while weakening the original matrix-matrix (m-m) connections.Jiang et al.[32]pointed out that large errors could be induced for EDFM to handle multiphase flow across fractures,and presented an improved pEDFM to effectively resolve the limitation.Rao et al.[33]developed a modified pEDFM that included a practical approach for determining fracture cell projection configuration and new fracture-fracture (f-f) connections to make the model framework more complete.Liu et al.[34] used pEDFM to model carbon dioxide geological sequestration in depleted complex-shape shale reservoirs.Wang et al.[35] presented a 3D pEDFM for compositional simulation of fractured reservoirs.Ren et al.[36,37] coupled pEDFM and XFEM to study flow in fractured reservoirs by considering geomechanics.Rao et al.[38] constructed a pEDFM based two-phase flow heat and mass transfer model,and found that pEDFM could eliminate the computational errors of temperature profiles in EDFM.These works indicated that,pEDFM had a better generalization of flow models than EDFM,and a lower computational cost than DFM.That is to say,pEDFM possessed the advantages of both DFM and EDFM,and hence deserved more attention.

However,this work points out that,in most cases,pEDFM still cannot effectively handle embedded flow barriers,implying that Tene et al.[31] had not completely accomplished their initial goal.By analysing the essence of pEDFM,this paper presents and proves an equivalence theorem that shows that the projection configuration of a fracture being physical is equivalent to it being topologically homeomorphic to fracture.Based on the equivalence theorem,pEDFM can yield physical results for flow barriers.Then,we present a practical method to handle complex geological conditions of reservoir model with a 3D Cartesian grid.

2 A Brief Review of pEDFM

2.1 Governing Equations

Take the immiscible multiphase flow problem as an example to illustrate the governing equations.For phasea,the mass conservation equation is written as

wherekis the absolute permeability;kra,Ba,μa,pa,saare relative permeability,volume factor,viscosity,the converted pressure by considering gravity and saturation of phasea;is the flow term of phaseafor wells (source or sink term);

The discrete form of Eq.(1) is always obtained by block-center finite volume method.For theith cell,it is written in Eq.(2).

The flux across the interface of two connected cells is approximated by two-point linear flux approximation,that is Eq.(3).

whereTijandλa,ijare the transmissibility and phase mobility of the connection.Tijis the half of the harmonic mean of two half-transmissibilities,that is Eq.(4).

λa,ijis expressed as Eq.(5).

The upstream and arithmetic average schemes are used for the terms (relative permeability) subject to the saturation and the terms (viscosity and volume factor) subject to pressure,respectively,these are in Eq.(6).

As a result,the ultimate discrete scheme of Eq.(1) is expressed as Eq.(7).

2.2 pEDFM Framework

Tene et al.[31] presented pEDFM as a way to expand traditional EDFM to efficiently treat fractures with a wide variety of permeabilities (i.e.,from flow barriers to highly conductive fractures).Jiang et al.[32] pointed out large errors will be induced for EDFM to handle multiphase flow across fractures,and presented an improved pEDFM to resolve the limitation.These studies revealed that pEDFM had higher physical model generalization than traditional EDFM,indicating that this new model required more attention.

Firstly,the basic ideas of pEDFM are briefly introduced below.As shown in Fig.1,pEDFM,like EDFM,divides the reservoir calculation domain using a structured or orthogonal mesh(recorded as background matrix mesh) and embeds fractures into the background mesh to create fracture cells divided by intersecting lines (surfaces) between matrix cells.In reservoir simulation,the reservoir boundary condition is always assumed to be closed.The essential difference between pEDFM and EDFM lies in the different connections between cells.

Figure 1:Sketch of neighboring matrix cells,a contained fracture cell and its projections on interfaces

As shown in Fig.1,the grey squares represent three matrix cellsmi,mjandmk,respectively.The yellow segment denotes a fracture cellfinmi.Two red segments represent thex-direction projection offon interfaceΓij(denoted byfx) andy-direction projection offon interfaceΓik(denoted byfy),respectively.

2.2.1 Construction of Additional f-m Connections

For these three matrix cells and a fracture cell,there is only one fracture-matrix (f-m)connection in classical EDFM,i.e.,f-mi.However,additionalf-mconnectionsf-mjandf-mkare constructed in pEDFM,and the areas offxandfyare considered as the flow areas off-mjandf-mkconnections,respectively.Besides,the originalf-miconnection needs to be weakened.

In pEDFM,forf-miconnection,the ‘matrix half’and ‘fracture half’transmissibilities are given as Eq.(8).

wherekmiandkfare the permeabilities ofmiandf,respectively,wfis the aperture off.Afis the area off,〈d〉f miis the average normal distance betweenfandmi,which is generally approximated by assuming that phase pressure within a matrix cell distributes linearly along the direction normal to the contained fracture cell,that is Eq.(9).

For addedf-mconnections,the areas offxandfyare used to determine the transmissibilities of these connections,by takingf-mjconnection as an example,the ‘matrix half’and ‘fracture half’transmissibilities are given as Eq.(10).

wherekmjis the permeability ofmj,Afxdenotes the area offx,df mjis the distance between the centers ofmjandf.

ThenTf mjis calculated as Eq.(11).

However,assuming that there is fluid flowing frommjtof,the physical flowing route ismjmi-f,so Jiang et al.[32] discussed that the transmissibilityTf mjshould take the effect ofmi(i.e.,Tmi) into account,which is Eq.(12).

Similarly,the transmissibility of thef-mkconnection can be obtained in Eq.(13).

wherekmkis the permeability ofmk,Afxdenotes the area offx,df mkis the distance between the centers ofmkandf.

2.2.2 Weakening of m-m Connections

Because the flow areas ofm-mconnections are blocked by additionalf-mconnections,the transmissibility of correspondingm-mconnections should be weakened.Takemi-mjconnection as an example,the effective flow area ofmi-mjconnection is equal to the area ofΓijminus the area offx,so the half-transmissibilityTiandTjforTmimjare given as Eq.(14).

whereAijis the area ofΓij,diis the distance from the center ofmitoΓij,djis the distance from the center ofmjtoΓij.

2.3 Determination of Projection Configuration of Fracture Cells

In the original paper by Tene et al.[31],there was no obvious statement about the determination criteria of the projection configuration of fracture cells,but it indicated that in each direction we should select the interface that is closest to the fracture-cell center,and this implicit criterion is simply denoted as the “smaller distance criterion”.Jiang et al.[32] pointed out that,smaller distance criterion may give an unphysical fracture-cell projection configuration in some cases when the distances between the fracture-cell center and two neighboring interfaces are the same in at least one direction.To resolve this problem,they presented an additional criterion that the intersections between the fracture-cell projections and the fracture cell should not be on the same side of the fracture-cell center,which is simply denoted as the “different side criterion”.Based on these two criteria,Rao et al.[33] presented a practical method namely the ‘microtranslation method’to efficiently select the projected interfaces of fracture cells to determine the fracture-cell projection configurations.

3 Not-Fully-Achieved Motivation

The motivation of the original pEDFM proposed by Tene et al.[31] is to modify EDFM to effectively handle cases when fracture permeability is lower than matrix permeability,especially those cases with flow barriers.

However,it is found that pEDFM still fails when dealing with flow barriers in most cases,that is,it can not accurately characterize the block effect of flow barriers.As shown in Fig.2,a flow barrier (plotted in yellow) is divided into three cellsf1,f2,andf3by matrix cellsm1,m2andm3.By following smaller distance and different sides criteria,obtained fracture-cell projection configurations are plotted in Fig.3,in which red segments,blue segments and green segments denote projection configurations off1,f2andf3,respectively.A1,A2andA3denote areas ofx-direction projections off1,f2andf3,respectively.The total projection areas on interfacesΓ12andΓ23can be calculated as Eq.(15).

Figure 2:Schematics of an example in which pEDFM fails to handle the flow barrier

Figure 3:Schematics of the fracture-cell projection configurations of the given example

Then,by following Eq.(7),the ‘matrix half’transmissibilities for them1-m2connection are Eqs.(16) and (17).

Consequently,the transmissibility of them1-m2connection is calculated as Eq.(18).

Similarly,the transmissibility of them2-m3connection is obtained as Eq.(19).

In Fig.4,the blue arrows denote the connections whose transmissibility is not 0.Imagining that there is fluid flowing into matrix cellm1along the positive direction ofx-axis,the fluid can still pass through the flow barrier along the path “l-m1-m2-m3-r”,which is unphysical,meaning the flow barrier is not effectively handled.

Figure 4:Schematics of the route along which fluid flow across the flow barrier

Fig.5 depicts another example.A flow barrier is divided into four cellsf1,f2,f3andf4by matrix cellsm1,m2,m3andm4.According to the same analysis as the previous example,fluid can still flow across the flow barrier along the route “l-m1-m2-m3-m4-r”,which is also unphysical.Furthermore,we give an example in Fig.6 to validate the above analysis.As is shown,a flow barrier whose geometry is the same as that in Fig.5 is embedded in the 60 m×60 m reservoir model,and a water injection well and a producer are at position (10,30) and position(50,30),respectively.The initial pressure is 20 MPa,and the constant injection rate is 1m3/d,and the constant BHP of the producer is 10 MPa.Fig.6 shows the water saturation profiles with different mesh sizes at the 300thday from the original pEDFM,which illustrates the unphysical computational results that injected water flow across the flow barrier and the unphysical results are independent of the mesh size.Moreover,it can be predicted similarly that when the number of matrix cells in a row continues to increase in Fig.5,the result will still be unphysical,which shows that,in most cases pEDFM cannot effectively handle the flow barrier,implying that Tene et al.[31] have not fully achieved their original motivation.

Figure 5:Schematics of the route along which fluid flow across the flow barrier when the matrixcell number increases to 4

Figure 6:Water saturation profiles with different mesh sizes (a) 5 m (b) 4 m (c) 2 m (d) 1 m

4 What is Physical Projection Configuration?

The issue of selecting a physical projection configuration was first noticed by Jiang et al.[32].They pointed out that in some special cases,unphysical projection configurations may be implemented by only following the “smaller distance criterion”.For an example shown in Fig.7,a fracture cell is located on the diagonal of a rectangular matrix cell,and the fracture-cell center coincides with the matrix-cell center.By following the “smaller distance criterion”,edge 2 or edge 4 can be selected inx-direction,and edge 1 or edge 3 can be selected iny-direction,resulting in four projection configurations (i.e.,edges 2-1,edges 2-3,edges 4-1,edges 4-3).However,Jiang et al.[32] pointed out that the projection configurations of edges 1-4 and edges 2-3 do not conform to the actual effect of fracture geometry orientation on flow and thus are unphysical.To resolve this limitation,Jiang et al.proposed the different sides criterion to avoid those unphysical projection configurations.

Figure 7:The example given by Jiang et al.

However,the examples in Section 3 show that “smaller distance” and “different sides” criteria are insufficient to obtain physical projection configurations.Moreover,it can be seen that these two criteria are subjective and lack rigorous analysis.Therefore,a profound and rigorous analysis of “physical fracture-projection configuration” is required,as well as clarification of its geometric nature.

4.1 Introduction of the Equivalence Theorem

In this paper,it will be pointed out that determining physical projection configuration is the most crucial step of pEDFM.In addition,when determining the physical projection configuration,it should focus on an entire fracture,not a single fracture cell.It should be noted that a fracture cell that does not penetrate the entire matrix cell (as shown in Fig.8) cannot completely block the flow across the matrix cell,so the “shortest distance” and “different sides” criterion can be directly applied to determine the projection configuration of the fracture cell.Therefore,unless otherwise noted,the fractures mentioned below are those removing fracture cells at both ends that do not penetrate the whole matrix cell.

In this paper,a necessary and sufficient condition is given to judge whether the projection configuration of a fracture is physical,that is,the equivalent geometric description of “projection configuration of a fracture is physical”,and is simply denoted as “equivalence theorem” in the following,

Equivalence theorem: the projection configuration of a fracture (regardless of the fracture cells which are at both ends of the fracture and do not penetrate the entire matrix cell) being physical is equivalent to it being topologically homeomorphic to the fracture.That is,from the geometric view,the topological properties of the projection configuration are the same as those of the fracture.

The conception of homeomorphism in topology is used in the statement of the above theorem.The following is a quick description of the conception,and detailed explanations can be found in the textbook [39].

Figure 8:Schematics of those fracture cells which are not penetrating the matrix cell

In topology,two topological spaces {X,TX} and {Y,TY} are called homeomorphic when a mapf:X Ysatisfies the following conditions: (i)fis bijection;(ii)fis continuous;(iii) the inverse offis continuous.

Then f is called the homeomorphism between {X,TX} and {Y,TY},and {X,TX} and {Y,TY} have the same topological properties.

Colloquially,topological space is a geometric object,and homeomorphism is to extend and continuously bend the object to make it a new object.So a square and a circle are homeomorphic,but a sphere and a torus are not.

To prove the equivalence theorem,this paper firstly studies the essence of pEDFM,and then obtains a lemma to prove the theorem.

4.2 Analysis of the Essence of pEDFM

EDFM handles fracture cells as the source or sink terms in matrix cells,whereas DFM uses unstructured matrix cells to match fractures,resulting in fractures being located on the matrix cell interfaces.Due to the two distinct fracture treatments,these two models exhibit distinct characteristics: (i) Because EDFM does not require a high-quality unstructured mesh to match the fracture geometry,it can be applied only to the physical models in which fractures are the source or sink terms of matrix cells (e.g.,production or injection simulation of a multi-stage fractured horizontal well),EDFM has significant limitations in the generality of physical models;(ii) DFM has good generalization of physical models,but it is difficult to generate a high-quality conforming mesh.

The essence of pEDFM is to combine the advantages of EDFM and DFM,that is,in the case of structured matrix cells (the same as EDFM),the fractures are projected onto the interfaces of matrix cells,and the essence of projection is to remove the fracture to the interfaces of matrix cells (the same as DFM).That is,pEDFM changes the “original model” shown in Fig.9a to the“approximate model”shown in Fig.9b,wherein the fracture cellfyis they-direction projection of the fracture cellfon the interfaceΓik,andfxis thex-direction projection offon the interfaceΓij.It can be seen that,in thex-direction,fxoccupies a part of the area ofΓij,so it is necessary to weaken the connection betweenmiandmj,and establish the additional connection betweenmjandfx.Similarly,in they-direction,themi-mkconnection needs to be weakened,and an additionalmk-fyconnection needs to be established.Thus,pEDFM is essentially a hybrid of EDFM and DFM,combining their advantages while avoiding their shortcomings.

Figure 9:Analysis of the essence of pEDFM from the view of only one fracture cell (a) original model (b) approximate model

If considering a fracture,take Fig.10 as an example,the essence of pEDFM is to change the “original model” shown in Fig.10a to the “approximate model” shown in Fig.10b.

Figure 10:Analysis of the essence of pEDFM from the view of a fracture (a) original model(b) approximate model

Therefore,a lemma can be intuitively obtained but is difficult to prove rigorously,that is:

Lemma: For a fracture (regardless of the fracture cells that do not penetrate the entire matrix cell at both ends),the essence of pEDFM is to continuously deform the fracture to a “new fracture” located on matrix cell interfaces,so as to obtain an approximate model of the original model for running simulation.

4.3 Proof of Equivalence Theorem

Using the above lemma,the equivalence theorem given in Section 4.1 can be proved.

Proof: According to the above lemma,the “new fracture” (i.e.,the projection configuration of the fracture) in the approximate model is obtained by continuous deformation of the fracture in the original model.Thus,from the geometric view,the projection configuration of the fracture is topologically homeomorphic to the fracture,i.e.,the projection configuration of the fracture and the fracture have the same topological properties.

For the example given by Jiang et al.[32],by following the equivalence theorem mentioned above,to judge whether the projection is physical,it is necessary to focus on the fracture,not fracture cells.Therefore,the fracture cell is stretched into two matrix cells (see Fig.11) to obtain an “imaginary fracture”.If edges 2-3 or edges 1-4 are selected as the projection configuration of the fracture cell,four types of the projection configurations of the “imaginary fracture” will be obtained in Fig.12,in which,(a),(b),and (c) are not connected,but the “imaginary fracture”is connected,so (a),(b) and (c) are not topologically homeomorphic to the fracture.(d) is the path-connected space,and its fundamental group is non-trivial,but the fundamental group of the straight-line fracture is trivial,so (d) is also not homeomorphic to the fracture.Because (a),(b),(c) and (d) in Fig.12 do not satisfy the equivalence theorem,the selection of edges 2-3 or edges 1-4 as the projection configuration of the fracture cell is unphysical.On the contrary,when edges 1-2 or edges 3-4 are selected,four other types of the projection configurations of the “imaginary fracture”will be obtained in Fig.13.It can be seen that (a),(b),(c),and (d) are all homeomorphic to the “imaginary fracture”,which means that edges 1-2 or edges 3-4 should be selected as the projection configuration of the fracture cell.

Figure 11:Extending the fracture to two matrix cells to obtain an imaginary fracture

Figure 12:Possible fracture-projection configurations when edges 2-3 or edges 1-4 are selected as projected edges (a) the first type;(b) the second type;(c) the third type;(d) the fourth type

Figure 13:Possible fracture-projection configurations when edges 1-2 or edges 3-4 are selected as projected edges (a) the first type;(b) the second type;(c) the third type;(d) the fourth type

5 Application of Equivalence Theorem to Effectively Handle Flow Barriers

In this section,the equivalence theorem will be used to analyze the example from Section 3 in which the flow barrier cannot be effectively handled by pEDFM,to illustrate that the proposed equivalence theorem is required to determine physical projection configuration in pEDFM that makes it a truly practical model.

By following “smaller distance” and “different sides” criteria,the projection configuration of three fracture cells can be obtained in Fig.14,where the yellow line segment represents the projection configuration off1,the green line segment is the projection configuration off2,and the red line segment is the projection configuration off3.

Figure 14:The projection configuration of each fracture cell according to the criteria of “smaller distance” and “different sides”

The geometry of the projection configuration of the entire flow barrier is:

As illustrated in Fig.15,the projection configuration is not connected and contains an intersection of two line segments,indicating that the projection configuration is not topologically homeomorphic to the flow barrier,which is a connected straight-line segment.According to the equivalence theorem,the projection configuration is not physical,which results in the flow barrier failing to perform its flow-blocking effect.

Figure 15:Geometry of the fracture projection configuration

It should be noted that the closer the flow barrier is to its projected configuration,the more accurate the approximate model is.Due to the proximity of the first or second type of projection configuration to the fracture,they should be chosen.

From Fig.16a,it can be seen that thex-direction projection off3is onΓ12,instead ofΓ23on which thex-direction projection off3should be if using “smaller distance” and “different sides”criteria to only determine the projection configuration off3;From Fig.16b,it can be similarly obtained that,thex-direction projection off3is onΓ23,instead ofΓ12.These results fully confirm the above statement that “when determining the physical projection configuration,it should focus on a whole fracture,not a fracture cell”.May as well select the first type of projection configuration,then,the total projection areas on interfacesΓ12andΓ23can be calculated as Eq.(20).

Figure 16:Four types of fracture-projection configurations (a) the first type (b) the second type(c) the third type (d) the fourth type

whereA12is the area ofΓ12.

Then,to calculate the transmissibility of them1-m2connection,Eq.(21) is obtained:

Similarly,to calculate the transmissibility of them2-m3connection,Eq.(22) is obtained:

Numerous connections between cells are denoted in the figure by arrows,with a red arrow indicating that the connection is not transmissible and a blue arrow indicating that the connection is transmissible.As can be seen,there is no path for the fluid to pass through the flow barrier;in other words,the flow barrier effectively blocks the flow from left to right,demonstrating that the proposed equivalence theorem resolves the limitation that flow barriers cannot be effectively handled in most cases.

To illustrate the correctness of the equivalence theory more vividly,the following test case is presented,and the performance of pEDFM based on and without the equivalence theorem is compared.As illustrated in Fig.17,the reservoir contains an “impermeable fracture” (i.e.,flow barrier),and the matrix cells are 2 m×2 m×2 m in size.The injector is operating at a constant flow rate of 1 m3/d,the producer is operating at a constant BHP of 10 MPa,and other physical parameters are provided in Table 1.Fig.18 compares the oil saturation profiles from pEDFM using the equivalence theorem and without using it.As can be observed,when the flow barrier is not based on the equivalent theorem,the flow barrier does not effectively block the flow.However,when the flow barrier is based on the equivalence theorem,the flow barrier effectively blocks the flow.

Figure 17:Schematic of the flow barrier and reservoir

Table 1:Physical properties of the formation and fluid in Example 1

Figure 18:Oil saturation profiles in the validity case (a) oil saturation profile obtained from original pEDFM (b) oil saturation profile obtained from pEDFM based on equivalence theory

6 Efficient Applications of pEDFM Based on Equivalence Theorem in Reservoir Models with Complex Geological Conditions

6.1 A Natural Idea

Because pEDFM based on the equivalence theorem can effectively handle cases where fracture permeability is less than matrix permeability,particularly flow barriers,an intriguing idea occurs naturally: pEDFM can be highly efficient when applied in reservoir cases with complex geological conditions (complex boundary shape or inclined interlayers) by embedding flow barriers to flexibly cut the reservoir model.

In comparison to obliterating the permeability of the matrix cells that contain fractures,there are two significant advantages to this approach that can be examined directly: (i) One is flexibility and efficiency.As we all know,when the reservoir boundary shape is complex or the reservoir interlayers are sloped,providing zero permeability to the non-reservoir cells is a tedious process.However,if applying the above idea,the reservoir model can be efficiently and flexibly cut to meet the needed boundary shape or inclined interlayers.(ii) The other is accuracy;because the cells employed in reservoir simulations are always a little coarse (10 m or 20 m),if we just assign them zero permeability,the cells become completely invalid,with the invalid region potentially considerably greater than the actual region.However,using the proposed idea,it is possible to obtain a more precise invalid region by modifying the aperture of the embedded barrier (or lowpermeable fracture).

6.2 A Numerical Example for Validation of the Above Idea

This section gives a numerical example to verify the above idea.Fig.19 illustrates the reservoir model in this example,which has a matrix-cell size of 2 m×2 m×2 m and a mesh division of 30×30×2.Four vertical flow barriers are used to cut the cuboid reservoir to construct a polygonal reservoir boundary,and a 20000 mD fracture is embedded in the reservoir,with a corresponding fracture or flow barrier cells plotted in red.A water-injection vertical well penetrates matrix cells (15,5,1) and (15,5,2) with a constant injection rate 1 m3/day,and an oil-production vertical well penetrates matrix cells (15,25,1) and (15,25,2) with a constant BHP 10 MPa.Other physical parameters of reservoir and fluids are listed in Table 1.

Figure 19:Schematics of the reservoir model

Figs.20 and 21 illustrate the oil saturation profiles over 50,100,200,and 300 days using different mesh sizes (2 m×2 m×2 m and 1 m×1 m×2 m).The results demonstrate that the four flow barriers do indeed define the reservoir model’s polygonal boundary and that a denser mesh can improve the accuracy of the results.The cost of computing is compared in Table 2 for various mesh sizes in this example,and it is discovered that as the number of meshes increases,the number of Newton iterations and calculation time increase as well.

Table 2:Computational parameters of this example

Figure 20:(Continued)

Figure 20:Oil saturation profiles over different days when the mesh size is 2 m×2 m×2 m (a) 50 days (b) 100 days (c) 200 days (d) 300 days

Figure 21:Oil saturation profiles over different days when the mesh size is 1 m×1 m×2 m (a) 50 days (b) 100 days (c) 200 days (d) 300 days

6.3 A Numerical Example with a Thin Interlayer in the Longitudinal Direction

As shown in Fig.22,the mesh division in this example is 30×30×3.The size of matrix cells in all directions is 10 m,and an inclined thin interlayer is embedded in the second layer.Water injection wells and production wells are only penetrated in the third layer,and other physical properties of the reservoir and fluid are the same as those in the above example.The oil saturation over 2000 days is shown in Fig.23.As a result of the thin impermeable interlayer,the injected water in the third layer cannot flow into the first layer,resulting in no obvious change in the oil saturation profile of the first layer,and the injected water in the third layer will flow in the area unaffected by the thin interlayer in the second layer.Table 3 shows the computational costs of this case.

Table 3:Computational parameters of this example

Figure 22:Schematic of the reservoir model of this example (a) 3D perspective (b) left view (c)front view

Figure 23:Oil saturation profiles over 2000 days in this example (a) the first layer (b) the second layer (c) the third layer

6.4 An Application Case with a Fracture Network

Fig.24 shows the 3D two-layer reservoir model in this example,including 8 hydraulic fractures,two natural fractures,and one fault (i.e.,flow barrier).Five hydraulic fractures penetrate the upper and lower layers,whereas three fractures and natural fractures only penetrate the lower layer.There are four water injection wells with a constant water injection rate of 25 m3/day and a fractured horizontal well with a constant bottom flow pressure of 10 MPa.The real field situation generally requires history matching,which is another important topic [40],and this paper focuses on the improvement of numerical algorithm.Therefore,this case directly gives the physical parameters of reservoir model,which are identical to those in Example 1.The simulation is set to run for 800 days.Figs.25 and 26 show the oil saturation distributions of the upper and lower layers at various times,while Fig.27 depicts the oil production and water production rates of the production well.As can be observed,water channeling of injected water occurs along natural fractures and hydraulic fractures,increasing the water cut of the production well,while the fault prevents injected water from flowing directly to the production wells.In addition,since the natural fracture and the hydraulic fractures near injection well #1 only penetrate the lower layer,the water channeling phenomena are more visible in the lower layer.Table 4 summarizes the computational costs of this example [40].

Table 4:Computational parameters of this example

Figure 24:The reservoir model and mesh division in this example

Figure 25:Oil saturation profiles of the upper layer (a) 50 days (b) 200 days (c) 500 days (d) 800 days

Figure 26:Oil saturation profiles of the lower layer (a) 50 days (b) 200 days (c) 500 days (d) 800 days

Figure 27:Production rates of the production well in this example (a) oil production rate (b) water production rate

7 Conclusions and Future Work

Throughout the whole paper,several key conclusions can be obtained:

(i) The original pEDFM fails to handle flow barriers in most cases due to unphysical fracture projection configurations,and determining the physical projection configuration of fractures is the critical step in pEDFM.

(ii) When determining the physical projection configuration,the focus should be on the entire fracture,not on individual fracture cells.Then,the equivalence theorem is presented and proved to explain the geometric equivalent nature of physical fracture projection configurations.The equivalence theorem-based pEDFM can truly handle cases where fracture permeability is lower than matrix permeability,especially flow barriers.

(iii) pEDFM can efficiently handle reservoir models with complex geological conditions (complex boundary shapes or inclined interlayers) by embedding flow barriers to flexibly ‘cut’the reservoir models.

In addition,it should be pointed out that,the current topological theory seems to be applicable only to disjoint fractures (i.e.,individual fractures),and a more general theory for intersecting fractures needs to be studied,which is deserving of significant future work.If topological theory,such as that discussed in this study,is still applied,it will be difficult to ensure that fluxes agree with the actual physical process of multiple fracture intersections.

Funding Statement: This work was supported by the National Natural Science Foundation of China (No.52104017),and National Key Research and Development Program of China (Grant No.2019YFA0705501),and State Center for Research and Development of Oil Shale Exploitation,and Cooperative Innovation Center of Unconventional Oil and Gas (Ministry of Education&Hubei Province),Yangtze University (No.UOG2020-17).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.