Well fracture modelling

Abstract

The present disclosure relates to numerical simulation of hydrocarbon reservoirs, in particular a method for modelling the fluid flow in fractures close to a well bore. The present invention describes a fast solver for multi-component multi-phase flow and its embodiment on a programmable electronic automaton.

Claims

1. Method for calculating flow conditions in fractures in the vicinity of a well bore, the well fracture being defined by a two-dimensional network of nodes and links describing flow paths intersecting a geological formation, the geological formation defined in terms of three-dimensional grid cells and allowing a fluid flow between said formation grid cells and the fracture, and from the fracture to the well, the method includes: a) iterative solution of the well flows contingent on boundary conditions from the formation, the well bottom hole pressure being varied from one iteration to the next until all of a plurality of well constraints are satisfied to a predetermined tolerance or a chosen iteration number; b) at each iteration, this involves a non-iterative solution of any fracture flows contingent on the well pressures and the boundary conditions from the formation grid cells, the pressure solution of the linear system of equations for any voidage flows in the fracture defines the individual drawdowns for each of a plurality of cell-fracture connections determining the component rates between the formation and the fracture.

2. Method according to claim 1, including application of non-linear models to affect the boundary conditions mentioned in claim 1, step b.

3. Method according to claim 1, wherein the adjustments of well pressure values in claim 1, step a involves a Newton iteration process.

4. Method according to claim 1, wherein said given precision is in the range of, and preferably within a threshold defined by a relative deviation of 10.sup.−6.

5. Method according to claim 1, wherein said well constraints are defined by data including individual pressure and fluid composition data for each grid cell.

6. Method according to claim 1, wherein the iteration is stopped after a predetermined maximum number of iterations, preferably in the order of 10.sup.2.

7. Method according to claim 1, wherein said linear set of equations is based on the single phase Darcy flow between neighboring network nodes.

8. Method according to claim 1, wherein the nodes include external nodes with links to internal nodes describing the flow from the formation into the network, the method including a step of, for each external node associated with a reservoir cell having a positive drawdown, calculating the phase/component flows to the adjacent internal node from the drawdowns determined in claim 1, step b and the said cell properties.

9. Method according to claim 1, wherein the nodes include external nodes with links to internal nodes describing the flow from the network into the formation, the method including a step of, for each external node associated with a reservoir cell having a negative drawdown, the phase/component flows are based on a stock tank model where the fractions for the well fluid as a whole are determined by averaging the fractions of all inflows, and then that split is applied to all outflows.

10. Method according to claim 1, wherein the fracture component flows into the well through each wellbore node are found by calculating a sum over all flows obtained in the previous two steps and attributing fractions according to the relative drawdowns of the wellbore nodes.

11. Method of displaying the static and dynamic properties of well fractures evaluated by application of method described in claim 1, including three-dimensional visualisation of the reservoir and wells and well fractures, and plotting of time-dependent phase and component rates for each well fracture.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) For a more complete understanding of the invention, reference is made to the following description and accompanying drawings, which illustrate the invention by way of examples:

(2) FIG. 1 illustrates the process used according to the invention; and

(3) FIG. 2 illustrates the node grid used according to the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

(4) According to the present invention, the fluid flows in the reservoir are solved iteratively. The details of the numerical method in general can be found in the book by D. Peaceman on Fundamentals of Numerical Reservoir Simulation [4] and the innovation presented here lies in the calculation of the fracture flows and their coupling to the formation and well flows. To evolve the total system from time t to time t+Δt, an iteration over the following instructions is performed until a solution within the chosen overall error tolerance is found or a maximum number of iterations has been reached. In the latter case, the time step size Δt is considered unsolvable and the procedure is repeated with a new time step Δt′<Δt. The iteration instructions are: 1. Calculate the net component and associated fluid volume changes for each reservoir cell. In- and outflows are from and to other reservoir cells or well connections, and their calculation may entail 2. Iterative solution of each well given the boundary conditions set by the state of all connected reservoir cells. Each well solution is obtained by varying the well pressure until all well control targets are met within a given tolerance or a maximum number of iterations has been reached. Solving each well may entail 3. A single direct (non-iterative) calculation of the internal well fracture pressures under the boundary conditions set by the well pressure at the first well iteration and the state of all reservoir cells connected to the fracture. The found solution is scaled to match the well pressure in each subsequent well iteration. This solution determines the component rates from the reservoir through the fracture into the well. The following section describes this method in detail.

(5) The well fracture is modelled as a network of nodes i bearing pressures p.sub.i and links (i,j) bearing voidage flows {dot over (V)}.sub.ij with unknown pressures p.sup.in on the internal nodes and unknown flows {dot over (V)}.sup.ex between internal and external nodes. See FIG. 2 for an illustration of the different node types.

(6) Solving a well fracture happens in three stages: 1. We find the internal node pressures by modelling single phase Darcy flow between neighbouring network nodes,
{dot over (V)}.sub.ji≡{dot over (V)}.sub.j.fwdarw.i=μ.sub.ji[p.sub.j−p.sub.i−ρ.sub.ji(d.sub.j−d.sub.i)],  (1) where μ.sub.ij=μ.sub.ji denotes the voidage conductance for the link between nodes i and j, and ρ.sub.ij=ρ.sub.ji is the average fluid density in that link. 2. For each external node j associated with a reservoir cell that has a positive drawdown, we calculate the component flows to the adjacent internal node i from the drawdowns determined in the first step,
q.sub.ji.sup.c≡q.sub.j.fwdarw.i.sup.c=ccf.sub.j.Math.m.sub.j.sup.c[p.sub.j.sup.e−p.sub.i.sup.i−ρ.sub.ji(d.sub.j−d.sub.i)],  (2) with ccf.sub.j the geometrical connectivity factor and m.sub.j.sup.c the component mobility for node j. The former is a purely geometrical number describing the intersection between fracture plane and reservoir grid cell j while the latter encodes the fluid properties in cell j. 3. For each external node with a negative drawdown, we base the phase/component flows on a stock tank model. That is, the fractions for the well fluid as a whole may be determined by averaging the fractions of all inflows, and then that split is applied to all outflows. 4. To obtain the fracture component flows into the well through each wellbore iso node, we sum all flows obtained in the previous two steps and attribute fractions according to the relative drawdowns of the wellbore nodes.
The conductances between internal and external nodes are calculated time-implicitly, while those between two internal nodes are updated time-explicitly.

(7) As usual, we allow for hydrostatic adjustment h.sub.ij ≡ρ.sub.ij(d.sub.i−d.sub.j)=−h.sub.ji between node depths d.sub.i and d.sub.j. The well fluid density is an explicit quantity in the dynamic reservoir model. Hence, it lags behind one time step.

(8) In step 2, we use the bulk component-in-phase mobilities of the respective cells in calculating the inflows into the fracture. Thus the solution can be considered first order multi-component multi-phase. For outflows we use the wellbore mobilities once all producing well connections have been summed up.

Solution Procedure

(9) In order to accomplish step 1, observe that mass conservation (or rather voidage conservation under the incompressibility assumption) requires

(10) $\begin{array}{cc}\underset{j\in N\left(i\right)}{.Math.}{\stackrel{.}{V}}_{\mathrm{ji}}=\left\{0,-{\stackrel{.}{V}}_i^e\right\}& \left(3\right)\end{array}$
to hold at internal and external nodes, respectively, where N(i) denotes the topological neighbourhood of node i. Inserting FIG. 1 into the l.h.s. of FIG. 3 gives

(11) $-\left\{0,-{\stackrel{.}{V}}_i^e\right\}=\underset{j}{.Math.}{\mu }_{\mathrm{ij}}\left(p_i-p_j\right)-\underset{j}{.Math.}{\mu }_{\mathrm{ij}}{h}_{\mathrm{ij}}=p_i\underset{j}{.Math.}{\mu }_{\mathrm{ij}}-\underset{j}{.Math.}{\mu }_{\mathrm{ij}}{p}_j-\underset{j}{.Math.}{\mu }_{\mathrm{ij}}{h}_{\mathrm{ji}}^T.$

(12) Let the N.sup.i internal and N.sup.e external node indices be sorted as {1, . . . N.sup.i, . . . N.sup.i+N.sup.e}. In shorthand notation, the combined equations (1) and (3) read

(13) $\begin{array}{cc}\left(\begin{array}{c}{0}^{\left(N^i\right)}\\ {\stackrel{.}{V}}^{\mathrm{ex}}\end{array}\right)+\mathrm{diag}\left(M.Math.H^T\right)=\left[-M+\mathrm{diag}\left(M.Math.e\right)\right]\left(\begin{array}{c}{p}^{\mathrm{in}}\\ p^{\mathrm{ex}}\end{array}\right).& \left(4\right)\end{array}$
Here, H denotes the matrix of hydrostatic head differences and M denotes the conductance matrix, i. e. it evaluates to μ.sub.ij for neighbouring nodes i, j and zero otherwise. The term diag(M.Math.e) is the diagonal matrix of row sums of M.

(14) The system (4) is semi-decoupled in that the first N.sup.i dimensions can be solved 170 independently from the last N.sup.e dimensions, and the solution p.sup.in of the former can then be inserted into the latter system to directly calculate {dot over (V)}.sup.ex. So the computational cost is determined by N.sup.i.

(15) To be precise, consider the following split into fracture and reservoir parts:

(16) $\begin{array}{cc}-M+\mathrm{diag}\left(M.Math.e\right)=\left(\begin{array}{cc}{F}^U& R^U\\ F^L& R^L\end{array}\right),& \left(5\right)\end{array}$
where F.sup.U is an N.sup.i×N.sup.i matrix, and the other are consequently of dimensions dim(R.sup.U)=N.sup.i×N.sup.e, dim(F.sup.L)=N.sup.e×N.sup.i, and dim(R.sup.L)=N.sup.e×N.sup.e. Note that R.sup.L is diagonal and F.sup.L=(R.sup.U).sup.T. Since R.sup.Up.sup.e does not contain any unknowns, the first N.sup.i equations of the system (4) read
diag(M.Math.H.sup.T)|.sub.1.sup.N.sup.i−R.sup.Up.sup.e=F.sup.Up.sup.in.  (6)
Internal Pressures
We solve for the internal node pressure vector p.sup.in as follows: 1. Solve
Lz=diag(M.Math.H.sup.T)|.sub.1.sup.N.sup.i−R.sup.Up.sup.e  (7) for z by forward-substitution. 2. Solve
L.sup.Tp.sup.in=z  (8) by backward-substitution.
In equations (7) & (8), L denotes the lower triangular matrix defined by the Cholesky decomposition
F.sup.U=L.Math.L.sup.T.  (9)
Mathematical Prerequisite
The main algorithmic task lies in decomposing matrix F.sup.U. The voidage mobility of any real fluid being positive, it must hold that μ.sub.ij=μ.sub.ji≥0 and μ.sub.ii=0∀i. Hence, F.sup.U follows to be a positive definite matrix which we can always decompose by an efficient Cholesky algorithm.
Alternative Calculation of Voidage Rates for Pure Injectors
In the case of non-crossflowing injectors, only the external network voidage flows q.sup.e are required. So the internal pressures need not be calculated explicitly. From FIG. 4 we find
q.sup.e=F.sup.Lp+R.sup.Lp.sup.e−diag(M.Math.H.sup.T)|.sub.N.sub.i.sub.+1.sup.N.sup.e  (10)
We see that the internal pressures p can be eliminated in order to find the linear functional q.sup.e[p.sup.e]

(17) $\begin{array}{cc}{q}^e={F^L\left(F^U\right)}^{-1}\left(\mathrm{diag}\left(M.Math.H^T\right){|}_1^{N^i}-R^U{p}^e\right)+R^L{p}^e-\mathrm{diag}\left(M.Math.H^T\right){|}_{N^i+1}^{N^e}=\left[R^L-{F^L\left(F^U\right)}^{-1}{R}^U\right]p^e+\left(\begin{array}{cc}{F^L\left(F^U\right)}^{-1}& 0\\ 0& -{N^e}_{}\end{array}\right)\mathrm{diag}\left(M.Math.H^T\right).& \left(11\right)\end{array}$

(18) We can avoid explicit inversion of F.sup.U altogether in order to evaluate FIG. 11 by calculating the following terms separately using the identity (L.Math.L.sup.T).sup.−1=(L.sup.T).sup.−1.Math.L.sup.−1=(L.sup.−1).sup.T.Math.L.sup.−1: 1. Y=L.sup.−1.Math.R.sup.U by substitution in the column vector equations Ly.sub.j=r.sub.j.sup.U. And since F.sup.L=(R.sup.U).sup.T, we find
F.sup.L(F.sup.U).sup.−1R.sup.U=Y.sup.T.Math.Y.  (12) 2. z=L.sup.−1 diag(M.Math.H.sup.T)|.sub.1.sup.N.sup.i by substitution in the system Lz=diag(M.Math.H.sup.T)|.sub.1.sup.N.sup.i. Hence,
F.sup.L(F.sup.U).sup.−1 diag(M.Math.H.sup.T)|.sub.1.sup.N.sup.i=Y.sup.Tz.  (13)
With these terms at hand, FIG. 11 simplifies to
q.sup.e=(R.sup.L−Y.sup.T.Math.y)p.sup.e+Y.sup.Tz−diag(M.Math.H.sup.T)|.sub.N.sub.i.sub.+1.sup.N.sup.e.  (14)
Multi-Phase and Multi-Component Flows
Once the internal node pressures are known, we determine the component rates directly from the set of drawdowns p.sub.j.sup.e−p.sub.i−ρ.sub.ji(d.sub.j−d.sub.i) utilising the two-point multi-component multi-phase flow expressions applied throughout the well solver. In order to calculate the flow derivatives with respect to the well bhp, we merely need the drawdown derivatives which are simply given by ∂p.sub.i/∂p.sub.j.sup.e for any of the external nodes j in the wellbore (since the latter only differ by constant heads). Taking the derivative ∂/∂p.sub.j.sup.e of FIG. 6 reads in index notation:

(19) $\begin{array}{cc}{F}_{\mathrm{ki}}^U\frac{\partial p_i}{\partial p_j^e}=-R_{\mathrm{kj}}^U.& \left(15\right)\end{array}$
This equation is solved by the same forward and backward substitution as in equations (7) & (8), i. e.

(20) $\begin{array}{cc}\mathrm{Lz}=-r_j{L}^T\frac{\partial p}{\partial p_j^e}=z,& \left(16\right)\end{array}$
where r.sub.j denotes the j-th column vector of R.sup.U.
Near Wellbore Modelling
The flow pattern far from the wellbore is assumed to be characterised by streamlines being perpendicular to the fracture plane. However, in particular for fractures along vertical wells, this assumption breaks down near the wellbore. Given the wellbore 220 radius an analytic approximation can be used to evaluate the deviation of the streamlines from the normal, which then enters (2) via the connectivity factors ccf.sub.j of the reservoir nodes according to their distance from the well. Likewise, the radial inflow directly into the wellbore (not through the fracture) is affected by the deviation of the streamlines from their unperturbed radial direction at the well radius. This effect can be taken into account by modifying the radial ccf appropriately.
Incorporating Non-Linear Effects
One can argue that multi-phase and non-Darcy effects are most prominent in the near wellbore region where the pressure changes appreciably, while the flow inside the fracture can be well approximated by the linear calculation above. To a good approximation, for a horizontal fracture, the flow pattern close to the well is radially symmetric no matter what the fracture looks like. Hence, the non-linear problem boils down to solving the flow in a disc, given boundary conditions at some outer effective well radius and the wellbore pressure at the inner radius. I. e. there exists some non-linear function
q=q(p.sub.bh,p.sub.well.sup.e),  (17)
containing one of the external pressures (the node at the wellbore of the linear problem in the previous section. This modelling is most relevant for horizontal wells where there is really only one network node close to the well—in contrast to the sketch of a vertical well in FIG. 2. The well fracture is then solved by satisfying the condition
q.sub.well.sup.eq(p.sub.bh,p.sub.well.sup.e),  (18)
which can be done by Newton's method or bisection. It is sped up by the fact that q.sub.well.sup.e linearly depends on p.sub.well.sup.e and thus requires only straightforward calculation.
Coupling with External Fracture Modelling Tools
The flexible design of this well fracture modelling facility both in terms of the user interface and the internal algorithm allows coupling with external fracture modelling tools. These could be used to supplement it with high-resolution fracture data such as orientation, height, width and permeability along the fracture. In return, results from the flow simulation could be fed back into the external tools to dynamically adjust the fracture properties in order to further enhance the accuracy of the combined model.

SUMMARY

(21) To summarize, the present invention relates to a method for calculating flow conditions in fractures in the vicinity of a well bore. The well fracture being defined by a two-dimensional network of network nodes and links describing the network of flow paths intersecting a geological formation. The nodes include internal and external nodes with links between internal nodes describing the flow within the network and links between internal and external nodes describing the flow into or out of the network. The well fracture thus allowing a fluid flow between said formation grid cells and the fracture, and from the fracture to the well. The geological formation may be defined in terms of three-dimensional grid cells including a set of fracture properties, the fracture properties are defining all available well fracture properties, involving geometrical data and fluid characteristics. The information may for example be provided through a graphical user interface or a textual definition file

(22) The method includes the steps of: a) iterative solution of the well flows contingent on boundary conditions from the formation, the well bottom hole pressure being adjusted from one iteration to the next until all well constraints are satisfied to a given precision or iteration number. b) at each iteration, this involves a non-iterative solution of the fracture flows contingent on the well pressures and the boundary conditions from the formation grid cells, the pressure solution, being based on a linear system of equations for the voidage flows in the fracture defined by the single phase Darcy flow between neighboring network nodes, defines the individual drawdowns for each cell-fracture connection determining the component rates between the formation and the fracture.

(23) Preferably the method utilizes the application of non-linear models to affect the boundary conditions mentioned in step, and the adjustments of well pressure values in a) may involve a Newton iteration process.

(24) The well constraint data may include individual pressure fluid composition data for each grid cell. The given precision of the well constraints is preferably within a threshold defined by a relative deviation of 10.sup.−6 or the iteration may be stopped after a number of iterations in the order of 10.sup.2.

(25) The network nodes include external nodes with links to internal nodes describing the flow from the reservoir to the two dimensional network, wherein for each external node being associated with a reservoir cell having a positive drawdown, calculating the phase/component flows to the adjacent internal node from the drawdowns determined in step and the said cell properties.

(26) For each external node associated with a reservoir cell having a negative drawdown, the phase/component flows are based on a stock tank model where the fractions for the well fluid as a whole may be determined by averaging the fractions of all inflows, and then that split is applied to all outflows.

(27) The fracture component flows into the well through each wellbore node, may be found by calculating a sum all flows obtained in the previous two steps and attribute fractions according to the relative drawdowns of the external wellbore nodes.

(28) The static and dynamic properties of well fractures evaluated by application of the method according to the invention may include three-dimensional visualization of the reservoir and wells and well fractures, and plotting of time-dependent phase and component rates for each well fracture.

REFERENCES TO RELATED APPLICATIONS

(29) [1] S. Kennon, S. Canann, S. Ward, and F. Eaton, “Method and system for modeling and predicting hydraulic fracture performance in hydrocarbon reservoirs,” Apr. 12, 2011. U.S. Pat. No. 7,925,482. 2 [2] G. Bowen, T. Stone, D. Bradley, and N. Morozov, “Multisegment fractures,” Jul. 12, 2016. U.S. Pat. No. 9,390,204. 2 [3] G. Bowen and T. Stone, “Multiphase flow in a wellbore and connected hydraulic fracture,” Mar. 25, 2014. U.S. Pat. No. 8,682,628. 2 [4] D. Peaceman, Fundamentals of Numerical Reservoir Simulation. Developments in Petroleum Science, Elsevier Science, 2000. 5