US20060277013A1

US20060277013A1 - Method for simulating fluid flows within a reservoir by means of a chimera type discretization

Info

Publication number: US20060277013A1
Application number: US11/444,474
Authority: US
Inventors: Chakib Bennis; Isabelle Faille; Lionel Ladouari; Roland Masson; Bertrand Moraisin
Original assignee: IFP Energies Nouvelles IFPEN
Current assignee: IFP Energies Nouvelles IFPEN
Priority date: 2005-06-02
Filing date: 2006-06-01
Publication date: 2006-12-07
Also published as: CA2548638A1; EP1729153B1; FR2886743B1; FR2886743A1; US7617083B2; NO20062498L; EP1729153A1

Abstract

Method for simulating flows in reservoirs by means of chimera type grids. The method first consists in generating a chimera grid made up of a set of elementary grids capable of taking account of the various characteristics of the flows and of the geometry (radial grid in the well zone and CPG grid for the reservoir for example) in such a way that their union covers the entire simulation domain and that there is a minimum overlap between two neighbouring elementary grids. Interpolation parameters allowing to connect the flows calculated on this grid to the flows calculated on the neighbouring grids are then determined on the border of each elementary grid. Finally, flows are simulated within the reservoir and the well, and the variables describing the flow on the border of the overlap zones are interpolated using the interpolation parameters. Application: simulators for hydrocarbon reservoirs for example.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention
The present invention relates to a method for simulating fluid flows within a reservoir from methods referred to as chimera methods; which is an alternative solution for reservoir simulation, allowing accounting for the radial directions of flow around wells and thus improving the calculation accuracy. Furthermore, this modular approach allows testing in an efficient and interactive way several development scenarios during modelling.
2. Description of the Prior Art
An oil reservoir is a porous and permeable subsoil formation associated with an impermeable cover that traps hydrocarbons. The order of magnitude of the lateral extension of a reservoir is one kilometer whereas the depth thereof is rather about ten meters. It is crossed by a variable number of geometric discontinuities such as wells or faults.
In oil exploration, survey of an oil reservoir requires simulation of the fluid flows within the reservoir and of the geometric discontinuities running therethrough in order to define the economic interest and the production modes of the field. These simulations are referred to as “reservoir simulations”.
It is therefore necessary to simulate fluid flows over the entire domain covered by the reservoir, by taking into account all the geometric discontinuities running therethrough. The simulation domain is then defined as the set of the reservoir and the geometric discontinuities. This simulation requires, on the one hand, physical modelling of the flows by a set of equations referred to as flow model and, on the other hand, discretization of these equations to solve the problem in an approximate way. The simulation domain is represented by means of a grid allowing to discretize it in space, thus allowing solution of the discrete equations and to provide an approximation of the physical solution. Simulation of fluid flows within a reservoir therefore requires two stages:

- discretizing the structure of the domain by means of a grid, and
- discretizing the flow equations according to the grid selected.

Grid generation is a crucial element for new-generation oil reservoir simulators. A grid allows description of the geometry of the geologic structure studied by means of a representation in discrete elements wherein the simulation is carried out. Better understanding of the physical phenomena involved requires 3D simulation of the multiphase flows in increasingly complex geologic structures, in the neighborhood of different types of singularities such as stratifications, faults, onlaps, channels and complex wells. All this complexity has to be taken into account first by the grid that has to reproduce as accurately as possible the geologic information in its heterogeneous character. Furthermore, good apprehension of the physical phenomena occurring in these complex structures requires a hierarchical and modular grid approach.
The following documents mentioned in the course of the description hereafter illustrate the state of the art:

- Alexander S. Williamson and John E. Chappelear, [
  ] “Representing Wells in Numerical Reservoir Simulation: Part 1 Theory”, SPE 7697, 1981
- Sophie Balaven-Clermidy, “Génération de Maillages H{dot over (y)}brides Pour la Modélsation de Réservoirs Pétroliers”, Thèse de doctorat, IFP, Ecole des mines de Paris, 2001.
- Chesshire et Henshaw, “Composite Overlapping Meshes for the Solution of Partial Differential Equations”, Journal of Computational Physics 90, 1-64, 1990.
- Shih, “Overset grids: Fundamentals and Practical Issues”, AIAA Applied Aerodynamics Conference. 3259, 2002.
- Lin, C. W., Smith, G. D., and Fisher, S. C.—Application of a Multiblock Grid Generation Approach to Ship Configuration, 3rd International Conference on Numerical Grid Generation in CFD, Spain, June 1991.

Wells are generally thin tubes having a radius of about ten centimetres and comprising several perforations (at the most permeable points). Scale problems in reservoir modelling are therefore flagrant: the reservoir/well flow models can therefore only be coupled models. This scale difference also has to be taken into account if the discontinuity is a fault, which is why it is crucial to obtain a fine description of the flows around the discontinuities.
There are several known approaches based on the use of a complete grid representing all of the simulation domain (reservoir+discontinuities). A commonly used method consists in gridding the reservoir part of the domain in a single block (several in case of faults) and to integrate the presence of other discontinuities such as wells in the equations. The presence of wells is considered as a source term defined by means of a productivity index Ip. This method is for example described in the document by Williamson (1981).
However, the necessity for a finer description of the flows around discontinuities such as wells leads to the elaboration of a modular approach where the various objects (reservoir and well for example) considered can be taken into account individually. Flow modelling around wells for example is carried out by means of the “well models”. These models allow simulation of multiphase flows in the well drainage area by adopting a small-scale description of the characteristics of the porous medium in the area around the well where the flows are fast.
Thus, in a more recent method developed by the assignee, the reservoir, in the area close to the discontinuity are gridded and the two grids are combined by creating a transition zone. The result is referred to as hybrid grid, as described in the document by S. Balaven-Clermidy (2001).
This method requires a considerable effort for generating the grid of the transition zone notably in three dimensions.
An alternative to this approach is the use of multi-block grids wherein the calculation domain is represented by a set of subdomains, gridded independently of one another. They can overlap one another or not. There are for example multi-block grids with non-coincident overlaps referred to as “chimera grids” or “Composite Overlapping Meshes”. This type of grid is described for example by Chesshire and Henshaw (1990). These chimera grids have been used for about ten years in the spheres of aeronautics (Shih, 2002) and hydrodynamics (Lin, Smith, G. D., and Fisher (1991) for complex geometries and in particular those with mobile bodies (missiles, helicopter blades, fluid structure coupling, . . . ). This type of grid allows adequate description of the flows on the border of discontinuities, but control of the subdomains as regards flows requires particular methods called chimera methods. However, the current chimera grid creation techniques and the associated chimera methods are not applied for reservoir simulation.
The method according to the invention allows simulation of fluid flows within the reservoir and its geometric discontinuities by modelling the structure of the domain by means of a chimera grid and by adjusting the fluid flow model to this grid by means of appropriate chimera methods.

SUMMARY OF THE INVENTION

The invention relates to a method for simulating fluid flow parameters within a heterogeneous medium forming a first subdomain crossed by at least one geometric discontinuity forming a second subdomain. The method comprises the following stages:
generating a chimera grid to construct a discretized representation of the subdomains, by forming at least a first grid for gridding the first subdomain and by forming at least a second grid for gridding the second subdomain, the two grids having at least one overlap zone;
determining interpolation parameters to estimate flow parameters in overlap zone while imposing continuity of the flow parameters between each subdomain;
simulating the flow parameters by solving the equations of at least one flow model in the subdomains and by interpolating the flow parameters at the level of the overlap zone using the interpolation parameters.
According to the method, generation of the chimera grid can comprise the following stages:
defining for the first subdomain an interpolation method allowing estimation of fluid flow parameters for each cell of the first grid included in the overlap zone;
defining an overlap parameter defining the size of the overlap zone;
determining a limit boundary in accordance with the interpolation method and the overlap parameter; and
determining the cells of the first grid to be interpolated according to the limit boundary, by selecting the cells all the points to be interpolated of which are located within the limit boundary.
Generation of the chimera grid can also comprise the following stages:
determining and deactivating cells, among the cells of the first grid to be interpolated, that will not be used a priori to interpolate flow parameters for cells of the second grid, by selecting cells all the neighboring cells of which are to be interpolated;
defining for the second subdomain an interpolation method allowing estimation of fluid flow parameters for each cell of the second grid included in the overlap zone;
checking that all the cells of the second grid to be interpolated are interpolable and adding the deactivated cells required for these interpolations.
The interpolation methods can be selected from among at least one of the following methods:
interpolation of the flow parameters at the center of the cells;
interpolation of the flow parameters at the center of the faces of the cells; and
interpolation of the flows on the faces of the cells.
If the geometric discontinuity is a well, the second grid can be a radial grid generated from at least one of the following geometric methods:
a geometric method based on the initial and the final radius, the number of cells along r and the desired geometric progression q;
a geometric method based on the initial and the final radius, the number of cells along r and the length of the border cells along r;
an automatic method locally based on the average length of the cell edges of the first grid.
Cells necessary for interpolation of the flow parameters at the level of the overlap zones can be determined by determining the convex envelope of cell centers containing a given point and by determining a cell containing the given point. The convex envelope can be determined by means of a sorting technique referred to as bucket sorting, and a cell containing the given point can be determined by means of a virtual grid.

BRIEF DESCRIPTION OF THE FIGURES

Other features and advantages of the method according to the invention will be clear from reading the description hereafter of embodiments given by way of non limitative example, with reference to the accompanying figures wherein:
FIG. 1 shows a radial grid;
FIG. 2 illustrates the definition of a cell of a radial grid;
FIG. 3 illustrates the definition of a limit circle with an interpolation at the center of the cells;
FIG. 4A illustrates part of a chimera grid with an overlap of 1;
FIG. 4B illustrates part of a chimera grid with an overlap of 2;
FIG. 5 shows an example of a 3D chimera grid constructed according to the method wherein an inclined well is integrated in a CPG grid;
FIG. 6 illustrates the interpolation method at the center of a cell;
FIG. 7 illustrates the interpolation method at the center of an edge; and
FIG. 8 illustrates the flow interpolation method at the center of an edge.

DETAILED DESCRIPTION

According to the invention, the method for simulating fluid flows within a reservoir can be broken up into three main stages:
1. Construction of a chimera grid representing the simulation domain, that is the set of the reservoir subdomain and of the subdomain of the discontinuities running through the reservoir;
2. Definition of interpolation formulas at the level of the overlap zone to provide continuity of the solution on the borders of the subdomains;
3. Simulation of the fluid flows within the simulation domain.
To simplify the description, the discontinuities are assumed to be wells, which allows to illustrate radial flows being taken into account, without limiting the invention to this type of discontinuity.
1. Construction of a Chimera Grid Representing the Simulation Domain
A chimera grid is a multi-block grid with overlap. These grids are a set of elementary grids covering the domain studied, having as only constraints a minimum overlap and similar cell sizes in the overlap zones. The main asset of this type of grid lies in the flexibility and the modularity of the definition of the grid, of the models and of the numerical methods of each block. On the other hand, they allow using structured elementary grids and therefore fast solvers in the blocks, and they are therefore suited for parallel calculation.
According to the invention, construction of a chimera grid can be broken up into three stages:

- Construction of a grid representing the reservoir;
- Construction of grids representing the wells that cross the reservoir;
- Definition of the overlap zone for the reservoir and well grids.

1.1 Construction of a Grid Representing the Reservoir
The reservoir can be discretized by regular Cartesian grids (identical discretization intervals in all the directions), irregular Cartesian grids (elongate in a direction) or Corner-Point geometry (CPG) grids:
A Cartesian grid is a grid whose vertices are of the form (x_i,y_j,z_k) with i, j and k natural numbers. A Cartesian grid can be generated by discretizing each direction of the bounding box identically or not according to the type of Cartesian grid required (regular or irregular).
A CPG grid can be defined starting from 12 curves defining each an edge of the grid. Discretization of each one of the curves then gives the desired grid. Such grids are for example described in French Patent 2,747,490 and corresponding U.S. Pat. No. 5,844,564] filed by the assignee.
The reservoir grid is thus defined by:
its bounding box, that is a Cartesian volume,
its discretization intervals in each direction.
1.2 Construction of Grids Representing the Wells Running Through the Reservoir
The wells can be discretized by a circular radial type structured grid allowing the drainage area of the wells to be modelled. FIG. 1 illustrates a radial grid. Such a grid is defined by:
the well trajectory (at least two points P₁and P₂defining an axis {right arrow over (a)}),
the initial (R_i) and the final (R_f) radius.
A cell can then be defined by (FIG. 2):
the minimum and maximum radius of this cell r1 and r2,
the minimum and maximum angle of this cell θ1 and θ2.
A point in this grid is then determined by its polar coordinates r and θ (FIG. 2):
r: distance to the axis of the well (in relation to the origin O),
θ: rotation around the axis of the well.
In order to better meet the requirements of the simulation according to the invention, three radial grid generation methods are presented:
1. A geometric method based on the initial and the final radius, the number of cells along r and the desired geometric progression q.
2. A geometric method also based on the initial and the final radius, the number cells along r and the length of the border cells along r (initial value of the geometric sequence).
3. An automatic method based on the initial and the final radius.
With Geometric Progression
The first method is the simplest one since it is based on the definition of the geometric sequence only.
With an initial radius R_iand a final radius R_f, we can find a geometric sequence of fixed common ratio q passing through n+1 points (R₀=R_iand R_n=R_f). Each point of this sequence is defined by: $R_{k} = R_{i} + h \cdot \frac{1 - q^{k}}{1 - q}$
The length of the first cell (from R₀to R₁), h, then just has to be determined as a function of the other parameters: $R_{f} = R_{i} + h \cdot \frac{1 - q^{n}}{1 - q} \Rightarrow h = \frac{(R_{f} - R_{i}) \cdot (1 - q)}{1 - q^{n}}$
It is now possible to evaluate the first formula for k=0 . . . n to know the radius of each discretization level of the well.
With Geometric Progression and Fixed Cell Length
The second method is still based on geometric sequences, but this time the value of h is fixed from the start. This method can be useful if one wants to check that the border cells of the grid generated have a given length. The missing parameter in this case is the common ratio q, which verifies a polynomial equation of degree n to be solved as a function of the other parameters.
With Fixed Cell Length
The last method is no longer based on the definition of the geometric sequences. It is the most automatic method: the parameters are calculated from the cells of the reservoir grid.
Let H be the average length of the edges of a reservoir cell. One decides to fix this length H as the length of the outer edge of the well border cells. The goal is here to generate a radial grid whose border cells have a size that is more or less equivalent to that of the reservoir cells.
H being fixed, we determine the discretization interval n_θ of the well according to the angle of rotation θ around axis {right arrow over (a)} of the well (see FIG. 2): $n_{θ} = \frac{2 π \cdot R_{f}}{H}$
Since the discretization interval has to be whole, the integer is taken that is the closest to n_θ, denoted by ñ_θ, then define {tilde over (H)} and {tilde over (h)} can be defined, the lengths along r of the outer and inner border cells of the well: $\tilde{H} = \frac{2 π \cdot R_{f}}{{\tilde{n}}_{θ}}$ $and$ $\tilde{h} = \frac{2 π \cdot R_{i}}{{\tilde{n}}_{θ}}$
These two values represent the initial and the final value of the progression of the cells along r.
The next stage determines the number of cells to be inserted for these dimensions to be best respected (an approximation of the size of the reservoir cells is sought).
Let the geometric interpolation function □s) defined by: $β (s) = \tilde{h} \cdot {(\frac{\tilde{H}}{\tilde{h}})}^{s}$
The value of n is deduced: $n = \frac{d}{\tilde{h}} \cdot [\frac{\tilde{h} - \tilde{H}}{\tilde{H} \cdot \ln (\frac{\tilde{h}}{\tilde{H}})}]$
Here again, the integer is taken that is the closest to n, denoted by ñ, since the number of cells is determined to be inserted along the radius of the well.
Then, the position of each discretization point in relation to the initial and the final radius has to be determined.
Therefore n=N is fixed and solve the equation for t: $N \cdot \frac{n}{\tilde{n}} = d \cdot \int_{0}^{t} \frac{1}{β (s)} ⅆ s$
for N=0 . . . ñ, ñ being the integer that is the closest to n. Multiplication by n/ñ simply allows readjustment of the value of N as a function of the rounding performed on n.
The value found for t allows location of discretization point N along the radius, from R_i.
1.3 Definition of the Overlap Zone of the Chimera Grid
Definition of the overlap zone of a chimera grid from the reservoir and well grids essentially consists in defining the cells for which no flow parameter simulation is necessary. After creating the reservoir grid and the radial grid, there are three types of cells in the chimera grid under creation:
the cells of the reservoir grid outside the overlap zone,
the cells of the reservoir grid of the overlap zone (which require an interpolation),
the cells of the well grid necessarily in the overlap zone (which require an interpolation).
Definition of the overlap zone comprises:
determining and marking the cells of the reservoir grid whose flow parameters have to be interpolated by means of the flow parameters of the cells of the well grid of the overlap zone; and
determining and removing the cells of the reservoir grid for which no flow calculation will be carried out.
The chimera grid thus constructed therefore has four types of cells:
the cells of the reservoir grid outside the overlap zone,
the cells of the reservoir grid of the overlap zone that have not been deactivated,
the cells of the well grid of the overlap zone,
the cells of the well grid outside the overlap zone (that require no more interpolation).
The number of cells requiring calculation of the flow parameters or interpolation is reduced while still ensuring a fine description of the flow phenomena at the level of the discontinuities.
For clarity reasons, the method is first presented in two dimensions. The three-dimensional generalization will be described later.
Definition of the overlap zone comprises 4 main stages:
stage 1: definition of a limit circle centered on the axis of the well,
stage 2: determination and marking, as a function of the limit circle, of the cells of the reservoir grid whose flow parameters will require interpolation,
stage 3: determination and deactivation of the interpolable cells of the reservoir grid that will not be used a priori to interpolate flow parameters of cells of the well grid and for which no flow calculation will be carried out,
stage 4: addition of some cells to the reservoir grid in order to guarantee interpolation on the radial (well) grid.
These various stages are described hereafter in detail:
Stage 1: Definition of a Limit Circle
Definition of the limit circle depends on two parameters:
a type of interpolation that it is desired to use on the subdomain of the reservoir, that is the type of interpolation used to define the flow parameters of the cells of the reservoir grid to be interpolated. Three types of interpolation can be considered: at the centre of the cells, at the center of the cell faces and at the flows,
an overlap parameter defining the size of the overlap. This parameter corresponds to the number of outer layers of the radial grid that are covered by the reservoir subdomain. This definition corresponds to a discrete overlap notion, that is corresponding to the grid.
The limit circle is defined as a function of the interpolation type as follows:
If, in the reservoir, interpolation is carried out with the first type of interpolation (at the center of the cells), the limit circle is the largest circle centered on the well containing cell centers of the radial grid (well grid). FIG. 3 illustrates this type of limit circle where MR represents the reservoir grid, MP the well grid and CL the limit circle. The points represent the border cells of the reservoir grid and the crosses represent the deactivated cells.
If, still in the reservoir, interpolation is carried out with the second type of interpolation (at the center of the faces), the limit circle is still the largest circle centered on the well, containing cell centers of the radial grid.
Finally, if interpolation of the reservoir is carried out at the flows, the limit circle is the largest circle centered on the well contained in the radial grid.
Then, this initial limit circle is redefined as a function of the second parameter, the overlap parameter. Definition of the overlap size of the chimera grid ensures that a certain number of layers of the radial grid are covered by the reservoir grid once the cells of the reservoir grid deactivated (stage 3). One therefore checks that a certain layer of the radial grid is interpolable by cells of the reservoir grid. A layer is referred to as interpolable if it comprises only interpolable cells. Thus, this layer and all the layers of greater radius are covered by reservoir cells. In fact, if a layer of the radial grid is covered by the reservoir grid, considering the geometry, it will also be the case for the layers whose radius is greater. Similarly, if a layer of the radial grid is interpolable, it is also the case for all the layers whose radius is greater.
For overlap of a single layer, one therefore takes the outer layer of the well, that is the layer corresponding to the initial limit circle as defined above (see FIG. 4A). On the other hand, for a greater overlap, one takes the layer such that the number of layers that are external thereto plus itself is equal to the desired overlap parameter, thus redefining the limit circle. In the example of FIG. 4B, the overlap is two.
It is thus possible to define, whatever the interpolation type, the “size” of the overlap which will be the number of circles beyond the initial limit circle.
Stage 2: Determination and Marking of the Interpolable Reservoir Grid Cells
The limit circle thus defined allows to identify the centers or the edges of the reservoir grid cells to be interpolated by well grid cells, as well as the centers or the edges of the cells in line with the axis of the well.
If, in the reservoir, interpolation is carried out with the first interpolation type (at the centre of the cells), the interpolable cells are the cells whose centers are within this limit circle.
If interpolation is carried out with the second interpolation type (at the center of the faces), the interpolable cells are the cells all the faces of which are interpolable, an interpolable face being a face whose center is within the limit circle.
Finally, if the reservoir is interpolated at the flows, an interpolable cell will be, as in the previous case, a cell all the faces of which are interpolable and therefore the middle of each edge of which is located within the limit circle.
All these interpolable reservoir grid cells are thus determined, then marked.
Stage 3: Rough Deactivation of Most of the Interpolable Cells
Most of the interpolable reservoir grid cells are roughly deactivated (a “hole” is made in the reservoir grid): in case of an interpolation at the center of the cells, all the cells adjacent to four interpolable cells (therefore all the interpolable cells, except those located the furthest from the center) are deactivated and, in case of an interpolation at the flows or at the centre of the edges, they are all deactivated. The minimum number of active cells for which interpolation on the reservoir grid is possible is thus obtained. In fact, in the first case, the cells on the border of the reservoir subdomain are cells interpolable through the well and, in the second case, the border of the reservoir subdomain consists of faces interpolable through the well.
Stage 4: Addition of Cells for Interpolation on the Well Border
At this stage, interpolation on the well grid cannot be guaranteed. In fact, some reservoir grid cells necessary for interpolation of the radial grid cells may have been deactivated. In this case, the reservoir cells necessary for interpolation on the well grid are reactivated by scanning one by one the radial grid cells in contact with the reservoir grid cells (referred to as well border cells) and the reservoir grid cells necessary for interpolation are reactivated.
It is therefore necessary to know the type of interpolation to be carried out to define the flow parameters of the well grid cells to be interpolated. It can be noted that, on the one hand, the latter stage depends on the type of interpolation performed on the well grid and, on the other hand, that the interpolation type is not necessarily identical to the interpolation type used on the reservoir subdomain.
The main stages of the algorithm used to define the overlap zone of the chimera grid can then be described as follows:
stage 1: definition of a limit circle, from the interpolation type to be used on the reservoir subdomain and an overlap parameter,
stage 2: determination and marking of the cells all the points of which to be interpolated (centers or edge centers) are located within the limit circle,
stage 3: determination and deactivation of the cells all the neighboring cells of which are interpolable (necessary for the well grid to be taken into account),
stage 4: checking that all the border cells of the well are interpolable and adding the deactivated reservoir grid cells necessary for these interpolations.
1.4 Methods of Detecting Interpolative Cells in 2D
The grid thus defined therefore is overlap zone of a set of reservoir cells to be interpolated and a set of well cells to be interpolated. The interpolable cells have been clearly identified at this stage. It is then necessary to define the methods allowing determination of the interpolative cells, that is the cells that will be used to interpolate the interpolable cells. In particular, the methods used to:
find the cells whose convex envelope of the centers contains a given point,
know in which cell a given point is located
have to be described.
To determine the convex envelope of cell centers containing a given point, the vertex that is the closest to the point in question has to be sought. Once this vertex found, the cells sharing this vertex just have to be scanned to recover the centres whose convex envelope contains the given point.
To find the closest vertex, the cell barycentre that is the closest to the point first has to be sought. The method can therefore use an efficient sorting technique referred to as bucket sorting. The vertex that is the closest among the 4 vertices (8 in 3D) is then determined.
To determine which cell contains a given point, the barycenter that is the closest to this point is first sought. Once the barycentre obtained, a virtual grid of the cell found, based on triangles, is created. These triangles have the barycenter and two vertices neighboring the cell as the vertices.
If a grid is not regular, it is possible to rapidly advance from cell to cell until the cell containing the point is found. In fact, the barycenter that is the closest to the point may not be that of the cell containing the point.
In this case, the list of the triangles of the virtual grid just has to be properly organized so as to know, when leaving a given triangle, towards which neighboring cell one is heading (top, bottom, left or right in 2D). For example, the triangles can be stored in this order. When leaving the cell through the last triangle, the cell located to the right of the current cell therefore has to be analized.
In order to know if a point is in a given triangle, its position in relation to each side of the triangle just has to be tested using the equations of the lines passing through these sides.
Creation of a virtual grid is very simple. The triangles just have to be created, for any cell, in the same order, for example top, bottom, left, right. The main thing is to know, for a given triangle, towards which neighbouring cell one is heading.
1.5 Extension to 2.5D
To carry out extension of the method to a 2.5D case, we consider a Cartesian reservoir grid and a radial well grid. The discretization interval of these two grids along z is the same. The ends of the well grid correspond to the limits of the reservoir grid. The method described above can then be applied to the upper face of the two grids (2D case) and the result is extended to the whole. In this case, an eliminated cell becomes a whole column of eliminated cells.
It is also possible to project the initial 2D result onto a surface of the form z=f(x,y) instead of remaining in a plane. This projection is then repeated along axis z to obtain another 2.5D grid type.
1.6 Generalization to 3D
The definition of the overlap zone of the chimera grid according to the invention can be directly extended to 3D, the limit circle becoming then a limit boundary (sphere, ellipsoid, . . . ). However, it is necessary to add a hemisphere to each end of the well so as to model the flows correctly, the latter being no longer radial but spherical at the ends of the well.
The method provided can therefore give correct results, whether the well is vertical or inclined in a Cartesian or CPG reservoir.
Detection of the interpolative cells in 2D can be directly extended to 3D.
In case if interpolation at the center of the cells or at the center of the faces, the eight cells the convex envelope of the centers of which contains the point to be interpolated have to be determined. These are, as in 2D, cells sharing the vertex that is the closest to the point.
As for interpolation at the flows, the cell containing the point to be interpolated just has to be determined. The method is the same, but the virtual grid of the cell then is made of tetrahedrons.
To define the tetrahedrons of the virtual grid, each face is divided into two triangles. Each face being defined by four vertices, the diagonal is taken passing through the vertex having the smallest index. Thus, the triangles are defined independently of the cells to which they may belong.
We can then define twelve tetrahedrons having the center of the cell as the vertex and a triangle of a face as the base. The tetrahedrons are then stored in a precise order allowing to know, when leaving one of the tetrahedrons, in which neighboring cell to carry out the new tests (same principle as in two dimensions). The same method as in 2D is used to know if a point is in a tetrahedron.
An example of a 3D chimera grid constructed according to the method is illustrated by FIG. 5, wherein an inclined well is integrated in a CPG grid.
2. Definition of Interpolation Formulas for the Overlap Zone
The main asset of chimera grids lies in the flexibility and the modularity of the definition of the grid, the models and the numerical methods of each subdomain. Besides, they allow to use structured elementary grids and therefore fast solvers in the subdomains, and they are naturally suited to parallel calculation. On the other hand, they require, in the overlap zone, the definition of points and of interpolation formulas for information transfer between the grids.
The continuous problem is written as follows:
find a solution to the equations of flow in the reservoir,
find a solution to the equations of flow in the well, and
find these solutions in such a way that these solutions and their derivatives link up continuously.
When discretizing, that is once the chimera grid is constructed, what is provided:
a system of equations for the reservoir (one equation per cell of the reservoir grid),
a system of equations for the well (one equation per cell of the well grid),
linking conditions which, in the chimera case, are interpolation formulas connecting unknowns of a grid to the unknowns of the other. These formulas are referred to as chimera methods.
Methods referred to as chimera methods (methods of modelling from a chimera grid) aim to solve partial derivative equations on subdomains whose grids overlap but are independent. In particular, no vertex, no cell center and no edge is common to the two grids, and edge conditions (flows or values at the points) have to be defined to solve the equations in each subdomain. Interpolation formulas are therefore selected to set the system.
It is crucial to find good interpolation formulas since the good behavior of the numerical scheme will depend on the validity thereof. This definition is essential for keeping the order of approximation of the schemes used in each subdomain. One of the goals is also to keep certain properties of the conventional schemes such as the stability (solution that remains bounded throughout the simulation time) and, in some cases, the conservativity (conservation of the mass of hydrocarbons in the reservoir) of the numerical schemes obtained on the elementary grids. The subdomain grids being made in such a way that solution of the problem is very efficient, the quality loss of the scheme (flow consistency, scheme stability and mass conservativity) will be directly linked with a defective interpolation method.
Three methods were retained, the first two methods for Dirichlet conditions (value of the solution imposed on the border) and the last one for Neuman conditions (flow value imposed on the border).
The interpolation formulas for the flow parameters on a subdomain are linear combinations of the flow parameters of the other subdomain. In general terms, these formulas allow interpolation of a value at the level of a point referred to as “point to be interpolated” which represents, for example, the center of a cell of one of the grids, from points referred to as “interpolative points” which represent, for example, the centers of cells of the other grid.
Interpolation at the Center of a Cell
In this case, what is referred to as the “point to be interpolated” is the center of a cell of one of the grids, and the “interpolative points” are the centers of cells of the other grid, used for interpolating the point to be interpolated.
The value at the center of a cell is interpolated by the values of the four cell centers whose convex envelope contains the point to be interpolated (8 centres in 3D).
Calculation of the interpolation coefficients is described below for interpolation within the reservoir grid, then within the well grid.
At the reservoir level: let A, B, C, D, E, F, G and H be the interpolative points and P the point to be interpolated. ABCDEFGH is assumed to be rectangular. We want λ_A, λ_B, λ_C, λ_D, λ_E, λ_F, λ_G, and λ_Hto be such that:
u _p=λ_A u _A+λ_B u _B+λ_C u _C+λ_D u _D+λ_E u _E+λ_F u _F+λ_G u _G+λ_H u _H
with u the value to be interpolated (typically a pressure).
We put:
{right arrow over (u)}={right arrow over (AB)}, {right arrow over (v)}={right arrow over (AF)} and {right arrow over (w)}={right arrow over (AD)}
$\hat{x} = \frac{\det (\vec{AP}, \overset{->}{v}, \overset{->}{w})}{\det (\overset{\overset{->}{->}}{u}, \overset{->}{v}, \overset{->}{w})}, \hat{y} = \frac{\det (\overset{->}{u}, \vec{AP}, \overset{->}{w})}{\det (\overset{\overset{->}{->}}{u}, \overset{->}{v}, \overset{->}{w})}$ $and$ $\hat{z} = \frac{\det (\overset{->}{u}, \overset{->}{v}, \vec{AP})}{\det (\overset{\overset{->}{->}}{u}, \overset{->}{v}, \overset{->}{w})}$
Here, the determinants 3×3 represent the volumes of the parallelepipeds constructed on the three vectors. We thus calculate the quotient of the volume corresponding to P in each direction on the total volume of the cell.
We thereafter define the interpolation coefficients as follows:
λ_A=(1−{circumflex over (x)})(1−ŷ)(1−{circumflex over (z)})
λ_B ={circumflex over (x)}(1−ŷ)(1−{circumflex over (z)})
λ_C ={circumflex over (x)}(1−ŷ)−{circumflex over (z)}
λ_D=(1−{circumflex over (x)})(1−ŷ){circumflex over (z)}
λ_E ={circumflex over (x)}ŷ(1−{circumflex over (z)})
λ_F=(1−{circumflex over (x)})ŷ(1−{circumflex over (z)})
λ_G =(1−{circumflex over (x)})ŷ{circumflex over (z)}
λ_H={circumflex over (x)}ŷ{circumflex over (z)}
At the well level: considering the radial grid used, it is necessary to work in polar coordinates. For each point, we have to determine r which is the distance thereof to the axis of the well, and θ which is the rotation around the axis of the well (see FIG. 2).
To find r, we simply determine the distance between the point and the axis of the well by projecting this point onto the axis. If the point is in a cap (intersection between a ball and a half space) at one end of the well, we take the distance to the vertex of the corresponding axis.
To determine θ, a local reference mark is used at the cell containing the point, as illustrated by FIG. 2. By projecting the point of coordinates (r, □) onto axis {right arrow over (u)}, we find point O, then abscissa x. Angle θ is deduced as follows: $θ = θ_{1} + ar \cos (\frac{x}{r})$
θ₁being known upon creation of the cells because the values of θ for each vertex of the well are stored when initializing the well.
The formulas used are then the same as in the Cartesian case of the reservoir.
FIG. 6 illustrates this type of interpolation for the reservoir and for the well, with:

- MR: reservoir grid
- MP: well grid
- MRI: interpolated reservoir grid cell
- MPI: interpolated well grid cell
- MRPI: reservoir grid cells used for interpolation of MPI
- MPPI: well grid cells used for interpolation of MRI.

Interpolation of the Unknown at the Center of an Edge
In this case, what is referred to as “point to be interpolated” is the center of the edge of a cell of one of the grids, and “interpolative points” are the centers of the cells of the other grid used for interpolating the point to be interpolated.
This interpolation is very similar to the previous one: the value of the unknown is interpolated with the same formulas, but this time at the center of the edge. Consequently, the following equation is also written in the border cells:
u _p=λ_A u _A+λ_B u _B+λ_C u _C+λ_D u _D+λ_E u _E+λ_F u _F+λ_G u _G+λ_H u _H
But, in this case, P is the center of a border edge and the Dirichlet conditions are applied thereto (value of the solution assumed on the border of the subdomain), which concretely amounts to introducing a border unknown (u_P) and an equation (the previous equation) in the system of equations to be solved to find the approximate solution.
The value at the center of an edge is interpolated by the values of the four cell centres whose convex envelope contains the point to be interpolated (8 centres in 3D), and calculation of the interpolation coefficients is identical to that carried out in the case of interpolation at the cell centers.
FIG. 7 illustrates this type of interpolation for the reservoir and the well, with the same notations as in FIG. 6.
Flow Interpolation at the Centre of an Edge
In this case, what is referred to as “point to be interpolated” is the center of an edge of a cell of one of the grids, and the “interpolative points” are the centers of the edges of the cell of the other grid used for interpolating the point to be interpolated.
The flow is here to be directly interpolated on the border of the subdomains. This interpolation is performed in the middle of the edge by a formula of mixed finite elements. The difference is that the latter no longer depends on the value of the unknowns at the centre of the neighbouring cells, but on the value of the flow at the edges of the interpolative cell.
To calculate interpolation coefficients α_i, we want to find here the flow on a face as a function of the flows of the faces of the interpolative cell: $F_{P} = \sum_{i = 1}^{6} α_{i} F_{i}$
It is then set: $F_{P} = \sum_{i = 1}^{6} α_{i} F_{i}$
where A is the surface area of the face to be interpolated, Ai the surface area of the interpolative face i, {overscore (n)} the normal external to the face to be interpolated and {right arrow over (v)}_ia base function defined on the interpolative cell as follows:
{right arrow over (v)} ₀=−(1−{circumflex over (x)})·{right arrow over (u)}
{right arrow over (v)} ₁=−(1−ŷ)·{right arrow over (v)}
{right arrow over (v)} ₂=−(1−{circumflex over (z)})·{right arrow over (w)}
{right arrow over (v)} ₃ ={circumflex over (x)}·{right arrow over (u)}
{right arrow over (v)} ₄ =ŷ·{right arrow over (v)}
{right arrow over (v)} ₅ ={circumflex over (z)}·{right arrow over (w)}
Vectors {right arrow over (u)}, {right arrow over (v)}, {right arrow over (w)} and values {circumflex over (x)}, ŷ, {circumflex over (z)} are the same as in the other two interpolation types.
By defining the {right arrow over (v)}_ithis way, we ensure continuity of the flows from one cell to the next because the flow on a face does not depend on the other faces. This can be simply explained by the fact that:
{right arrow over (v)} _i ·{right arrow over (n)} _j=δ_ij=0 if i≠j
hence: $α_{i} = \frac{A}{A_{i}} \overset{⟶}{v_{i}} \cdot \vec{n}$
If we evaluate for example the flow of face 0 as a function of faces 0, 1, 2, 3, 4 and 5, we will find F₀because only {right arrow over (v₀)},{right arrow over (n₀)} will be non-zero.
FIG. 8 illustrates this type of interpolation for the reservoir and the well, with the same notations as in FIGS. 6 and 7.
3. Simulation of Fluid Flows Within the Simulation Domain
Simulation of the fluid flows consists in assigning to each cell of the previously constructed chimera grid flow parameter values by defining a numerical scheme and by determining the interpolations inherent in the chimera nature of the grid.
The physical model used for modelling fluid flows in the reservoir and in the radial zone around the well can be a compressible and non-compositional two-phase model well known to the man skilled in the art and described for example in the following document:
Peter Bastian, “Numerical Computation of Multiphase Flows in Porous media”, Habilitationsschrift, Technische Fakultät, Christian-Albrechrs-Universität Kiel, 1999.
The numerical scheme used to model these equations can be an upstream off-centered finite-volume type scheme known to the man skilled in the art and described for example in the following document:
Thierry Gallouët, Raphaëlle Herbin and Robert Eymard, “The Finite Volume Methods”, in Philippe Ciarlet and J-L Lions, editors, Handbook of Numerical analysis, pages 713-1020, North Holland, 2000.
The physical model used for modelling fluid flows within the well can be taken into account as an edge condition (imposed pressure or flow).
As for the overlap zone, we use the interpolation formulas from the selected chimera methods described above.

Claims

1. A method for simulating fluid flow parameters within a heterogeneous medium forming a first subdomain, crossed by at least one geometric discontinuity forming a second subdomain, comprising:

generating a chimera grid to construct a discretized representation of the subdomains, by forming at least a first grid for gridding the first subdomain and by forming at least a second grid for gridding the second subdomain, the two grids having at least one overlap zone;

determining interpolation parameters to estimate flow parameters in the overlap zone while imposing continuity of the flow parameters between each subdomain; and

simulating the flow parameters by solving the equations of at least one flow model in the subdomains and by interpolating the flow parameters at the level of the overlap zone using the interpolation parameters.

2. A method as claimed in claim 1, wherein generation of the chimera grid comprises:

defining for the first subdomain an interpolation method allowing estimation of fluid flow parameters for each cell of the first grid included in the overlap zone;

defining an overlap parameter defining the size of said overlap zone;

determining a limit boundary in accordance with the interpolation method and the overlap parameter; and

determining the cells of the first grid to be interpolated according to the limit boundary, by selecting the cells all the points to be interpolated of which are located within the limit boundary.

3. A method as claimed in claim 2, wherein generation of the chimera grid comprises:

determining and deactivating cells, among the cells of the first grid to be interpolated, that will not be used a priori to interpolate flow parameters for cells of the second grid, by selecting cells all the neighboring cells of which are to be interpolated;

defining for the second subdomain an interpolation method allowing estimation of fluid flow parameters for each cell of the second grid included in the overlap zone; and

checking that all the cells of the second grid to be interpolated are interpolable and adding the deactivated cells required for these interpolations.

4. A method as claimed in claim 2, wherein said interpolation methods are selected from among at least one of the following methods:

interpolation of the flow parameters at the center of the cells;

interpolation of the flow parameters at the center of the faces of the cells; and

interpolation of the flows on the faces of the cells.

5. A method as claimed in claim 1 wherein, the geometric discontinuity being a well, said second grid is a radial grid.

6. A method as claimed in claim 5, wherein the radial grid is generated from at least one of the following geometric methods:

a geometric method based on the initial and the final radius, the number of cells along r and the desired geometric progression q;

a geometric method based on the initial and the final radius, the number of cells along r and the length of the border cells along r; and

an automatic method locally based on the average length of the cell edges of the first grid.

7. A method as claimed in claim 1, wherein cells required for interpolation of the flow parameters at the level of the overlap zones are determined by determining a convex envelope of cell centers containing a given point and by determining a cell containing the given point.

8. A method as claimed in claim 7, wherein the convex envelope is determined by means of a sorting technique referred to as bucket sorting.

9. A method as claimed in claim 7, wherein a cell containing the given point is determined by means of a virtual grid.

10. A method as claimed in claim 3, wherein said interpolation methods are selected from among at least one of the following methods:

interpolation of the flow parameters at the center of the cells;

interpolation of the flows on the faces of the cells.