The present invention relates to the field of oil exploration, and more particularly the study of fluid flows within a subterranean formation. The method in particular relates to the study of the displacements of fluids such as hydrocarbons in a deposit, or subterranean reservoir, traversed by faults which is carried out with basin modeling, reservoir simulation, and/or simulations of CO2 injection into a subterranean formation.
The goal of basin modeling is to reconstruct the geological history of a sedimentary basin and its oil systems to help to locate hydrocarbon traps known as reservoirs, to estimate their quantity and quality, and finally, to evaluate the risks of encountering pressure excesses during drilling. Reservoir simulation studies the evolution over time of the proportions of water, gas and petroleum in the reservoir in order to assess profitability, to validate or to optimize the position of the wells ensuring operation of the reservoir exploitation. In an era where durable development and the protection of the environment are becoming paramount, a third study related to oil exploration performs simulations for the injection of CO2 into a porous subterranean medium.
Basin modeling, reservoir simulation and CO2 simulation are techniques based on the simulation of flow in a porous medium. A representation of the subterranean medium, also called a reservoir model, constitutes a mockup of the sub-soil, representative both of the structure and of its behavior. Generally this type of mockup is represented on a computer. A flow simulator is software making possible, inter alia, modelling the production of a deposit as a function of time, on the basis of measurements describing the deposit, on the basis of a representation of the subterranean medium. These simulations are carried out by a system of partial differential equations by finite volume procedures on a meshed model of the subterranean medium concerned.
Today, exploration is concerned with zones with complex geometry where faults are numerous and their positions random. The automatic construction of a 3D mesh capable of representing this geometric complexity is the first indispensable step in the formulation of the simulation software, the flow simulator, for such a medium, on the basis of the horizons which delimit the various geological layers, and of the faults which cut the horizons. The horizons and the faults are provided as surfaces, triangulated on the basis of a net of points. These points generally result from seismic surveys. Having regard for the flow simulators used in the industry, it is necessary that the mesh comprises as many hexahedral elements as possible so as to allow a better result of the simulation.
For example, three basin simulation software packages have been developed at IFP Energies nouvelles, namely Temis 2D®, Temis 3D® and Ceres 2D®. The first two software packages make possible modeling of the basin but neglect the presence of faults, whereas the last is suitable for taking faults into account in the two-dimensional case.
To obtain a three-dimensional mesh with as many hexahedral elements as possible, it is necessary that the two-dimensional grids making possible construction of the three-dimensional mesh comprise a maximum of elements of optimal shape, that is to say substantially rectangular, and that the other elements of the grid not be degenerate, that is with angles far removed from right angles.
Furthermore, a procedure has been developed, referred to as the constrained grid procedure and described in particular in patent application FR 2 948 215, for constructing a suitable mesh, representing a subterranean medium comprising several faults. This procedure is based on carrying out the following steps:
The meshes constructed with the constrained grid procedure are very well suited to the problems of basin studies, since the procedure processes a large number of mesh cells with complex geometry. Moreover, the procedure makes it possible as it were to project each horizon topology onto another horizon.
The procedure for fitting the faults, proposed by this procedure, advances over a fault trace from one extremity to the other while fitting the closest nodes on the traces of the faults while avoiding creating degenerate quadrilaterals. The procedure for fitting the faults uses geometric criteria, that is criteria based essentially on the closest distance by incrementing as one advances over the nodes belonging to the fault trace. This requires a uniform dispersion of the quadrilaterals neighboring this trace. However, the procedure for generating the regular two-dimensional grid is suitable only for the generation of the meshes of quadrilaterals for geometric maps and it actually complies with the non-degeneracy of the mesh cells. The dispersion is not taken into account at all and the quality of the deformed mesh cells in proximity to the faults is not optimal. This can be illustrated by the case where there are more mesh cells on one side of a trace than on the other. Moreover, after the fitting, the deformation procedure is employed only on the mesh cells neighboring the trace with this optimization being insufficient. Furthermore, in the zones of folds of the regular two-dimensional grid, the mesh cells do not possess a shape which is optimized for the simulations. Thus, the procedure described in this prior art document does not make possible construction of an accurate model, suitable for basin simulation.
The method according to the invention relates to a method for exploiting a determined subterranean medium on a basis of a meshed representation of the subterranean medium. The mesh quality of the representation of the medium is optimized by a deformation of a two-dimensional grid while minimizing the gradient of displacement of the nodes of the grid. The method according to the invention can be used as a supplement to the procedure described in patent FR 2 948 215 in order to form a more accurate model, suitable for simulations.
The invention relates to a method for exploiting a subterranean medium, in which the subterranean medium is exploited according to an exploitation scheme defined on the basis of a representation of the medium, on the basis of a two-dimensional grid representing the subterranean medium. The quality of the mesh cells of the grid is optimized by displacing the nodes of the grid by a procedure for minimizing the deformation of the grid which comprises:
According to the invention, the position of the node corresponding to the rigid node in the grid to be optimized is determined in accordance with geological or geometric considerations.
Advantageously, a rigid node is a node belonging to the boundary of grid to be optimized or a node corresponding to a trace of a fault of the subterranean medium in the grid to be optimized.
Furthermore, the invention relates to a method for exploiting a subterranean medium according to an exploitation scheme defined on a basis of a representation of the medium, the representation of the medium comprising a three-dimensional mesh, the subterranean medium comprising at least one sedimentary layer traversed by at least one fault, the layer being delimited vertically by two geological horizons. For this method, the following steps are carried out:
According to the invention, the procedure for minimizing the deformation of the grid furthermore comprises a final step of testing validity of the mesh and if the mesh is not valid, steps (1) to (3) are repeated.
Advantageously, the validity test detects the folded two-dimensional surfaces of the optimized two-dimensional grid and if folded surfaces are detected, then the displacement field is divided by two while repeating steps (1) to (3).
Preferably, the procedure for minimizing the gradient of the displacement field furthermore comprises a step of converting the grid into a triangular mesh for the steps of displacing the nodes.
According to one embodiment of the invention, the reference grid is generated by carrying out the following steps:
Preferably, the regular grid is generated by the following steps:
In an advantageous manner, for step v) before displacing the closest end, a test is carried out to verify that the displacement does not produce a quadrilateral having at least one angle which is greater than a fixed angle threshold and if such is the case, the other end of the edge which is intersected is displaced.
Preferably, the regular grid is transformed into a three-dimensional gridded surface by carrying out the following steps:
According to another embodiment, links are created between the nodes by linking each node of each three-dimensional gridded surface having like coordinates i, j, and if a fault intersects this link, a node is linked with the fault by considering a direction of a neighboring node.
The invention also relates to a computer program product downloadable from a communication network and/or recorded on a support readable by computer and/or executable by a processor; the computer program product comprises program code instructions for the implementation of the method such as described above, when the program is executed on a computer.
Other characteristics and advantages of the method according to the invention will become apparent on reading the description hereinafter of nonlimiting examples of embodiments, while referring to the appended figures described hereinafter.
The invention relates to a method for exploiting a subterranean medium, in which the subterranean medium is exploited according to an exploitation scheme defined on the basis of a representation of the medium, established on the basis of a two-dimensional grid representing the subterranean medium. The optimization of the quality of the mesh cells of a regular two-dimensional grid representing a subterranean medium is obtained by displacement of the nodes of the grid by a procedure for minimizing the deformation of the grid which comprises the following steps:
a) generation of a reference grid;
b) rigid node displacement;
c) construction of the optimized two-dimensional grid.
The method according to the invention can furthermore comprise a step d), for which a test of validity of the optimized mesh obtained is performed.
The ability of a mesh cell to give coherent results during simulations is called the quality of a mesh cell. A rectangular mesh cell is the type of mesh cell whose quality is highest; conversely a so-called degenerate mesh cell (with small angles) possesses poor quality. This quality can be determined by computing the ratio of the value of the smallest angle of the mesh cell to the value of a right angle (90° or
$$
Thus, a rectangular mesh cell will have a quality value of 1.
Step a) generation of a reference grid
The goal of the method according to the invention is to obtain a mesh optimized in terms of quality. Accordingly, a two-dimensional grid representing the subterranean medium is deformed by displacing certain nodes of the grid in an imposed manner. The other nodes are displaced as a function of their location in the grid. To obtain this result, a reference grid covering the existing two-dimensional grid is generated initially. This reference grid comprises an equal number of mesh cells to that of the grid to be optimized. Furthermore, the reference grid is formed solely of rectangles and the dimensions of the reference grid are greater than the dimensions of the existing two-dimensional grid. Consequently, the reference grid is of good quality and covers the entirety of the grid to be optimized.
A plane having two families of perpendicular straight lines delimiting rectangles, preferably squares, is called a grid or a framework (set of principal lines and points of a figure).
In one embodiment of the invention, the reference grid is generated by carrying out the following steps:
The goal of this process is to reduce the initial displacement field and therefore to improve the accuracy as well as the computation time.
Preferably, the quadrilaterals of the reference grid are transformed into triangles to facilitate the construction of the optimized grid. Here this entails separating each quadrilateral into two triangles by passing through one of their diagonal. The reference grid of
Step b) displacement of the rigid nodes
For the method for optimizing the mesh according to the invention, the displacement of certain nodes of the reference grid, called rigid nodes, is imposed in accordance with geological or geometric considerations. The rigid nodes can be in particular the nodes of the boundaries of the grid or the nodes corresponding to the traces of faults, whose displacements are known.
These imposed displacements displace the nodes of the reference grid toward the location of the corresponding nodes in the two-dimensional grid. The imposing of the displacement field on the boundaries makes possible retention and compliance with the topology of the surface whose grid is to be optimized. For
Step c) construction of the optimized two-dimensional grid
This step deforms the reference grid while minimizing deformation. The minimum deformation makes it possible to preserve the quality of the reference grid while preserving the topology of the grid to be optimized representative of the subterranean medium. The optimized grid is therefore at the same time representative of the medium and of good quality. Indeed, through this deformation, the quality of the reference grid is diffused and distributed to the grid to be optimized. The minimum deformation according to the invention corresponds to the minimization of the gradient of the displacement field of the nodes of the reference grid toward the surface to be optimized.
During this step, all the nodes of the grid whose displacement is not imposed (the non-rigid nodes) are therefore displaced, so as to obtain a good quality mesh representative of the subterranean medium. This optimizes not only the mesh cells close to the boundaries and traces but also all the mesh cells of the grid.
For the preferential embodiment, for which the quadrilaterals have been triangulated, the inverse transformation is operated during this step to obtain quadrilaterals.
Given the points M=(x,y), M′=(x′,y′)∈R^{2}, a displacement field {right arrow over (u)}=(u_{1}(dx,y),u_{2}(x,y)) from M to M′ is defined by the following transformation:
$$
According to the invention, the displacement field {right arrow over (u)} is imposed on the nodes of the traces of faults and those of the boundary of the surface. The choice of this displacement field is defined on the basis of the deformation from a reference surface S_{R }formed by a regular mesh of very good quality to the surface S_{o }whose mesh is to be optimized. Thereafter, this deformation which is diffused inside the surface is minimized by smoothing all the internal mesh cells. This smoothing must comply with a certain graduation of the initially deformed quantity which is attenuated on moving further away from the neighborhood of the traces of the faults and the boundary.
The embodiment of the computation of the displacement field and of its minimization is detailed as an annex to the description.
Step d) test of validity of the mesh
After having deformed the reference grid, it is possible to perform a validity test on the grid. If the grid obtained is not valid, steps a) to c) are repeated. For example, the test can detect the resulting folded surfaces of the optimization method. Indeed, during the step of displacing the nodes, it is possible that folds may be generated in the grid. These folds are not suitable for the simulations and also for the construction of a three-dimensional mesh. If, by this test, folds are detected, then steps a) to c) are repeated while dividing the displacement field by two.
The invention also relates to a method for exploiting a subterranean medium according to an exploitation scheme defined on the basis of a representation of said medium. The representation of the medium comprises a three-dimensional mesh. The subterranean medium comprises at least one sedimentary layer traversed by at least one fault with the layer being delimited vertically by two geological horizons. According to the invention, the mesh of the representation of the subterranean medium is generated by carrying out the following steps:
Steps a), b), c), e), g) and h) are known, in particular from patent document FR 2,948,215. Consequently for these steps, only the essential characteristics will be discussed.
Step a) discretization of the geological horizons
During this step, the two geological horizons delimiting each sedimentary layer are discretized by two triangulated three-dimensional surfaces.
Step b) transformation into triangulated two-dimensional surface
This step unfolds in an isometric manner the triangulated 3D surfaces to obtain a triangulated 2D surface and makes possible stitching together the tears in the faults present on the triangulated 3D surfaces.
A transformation of a 3D surface into a 2D surface (an operation which flattens out a surface) is called “unfolding”. This transformation is “isometric” when it preserves measures, such as the lengths of edges, and therefore surfaces, such as the areas of triangles.
Tools and procedures are known for carrying out an isometric unfolding of a surface. The method according to the invention does not depend on the type of unfolding procedure.
Step c) generation of a regular two-dimensional grid
This step generates a regular 2D grid, on the basis of the periphery of the triangulated 2D surface.
This regular grid is generated by the following steps:
For all the grids, a corner is matched and the boundary curves are oriented in such a way that the domain is situated on the left. Here domain is a mathematical term meaning the geometric shape that it is desired to mesh. This step serves to orient all the horizons in the same manner to facilitate the final step which is the matching of the horizons. Next, the boundary curves are subdivided into N or M segments of constant length with N and M being the same as the grid beneath.
Given four parametric curves f1(u), f2(u), g1(v), g2(v) (0≤u, v≤1) defining four connected boundaries, the Coons formula computes the tightest surface which passes through the four boundaries (which interpolates these boundaries):
$$
The vertices P(i,j) are the four corners of the surface.
It is then possible to obtain a regular mesh N×M of this surface by simple sampling, taking as vertices the points S(i,j) corresponding to the ui=i/(N−1), vj=j/(M−1), i=0 to N−1, j=0 to M−1. The formula restricted to the points of the mesh becomes:
$$
Step d) optimization of the quality of the mesh cells
With the grid optimization steps described hereinabove, a two-dimensional grid is generated and optimized in terms of quality of the mesh. For the application of the grid optimization steps described above, the rigid nodes correspond to the nodes of the boundaries of the grid to be optimized. The displacement of these rigid nodes is therefore imposed on the nodes of the boundaries of the grid to be optimized.
Step e) fitting of the faults
Taking a fault into account is by deforming the mesh in such a way that it conforms with the reality of the subterranean medium to be meshed. Indeed, if a fault passes through a medium, generally the sedimentary layers of this medium are broken and deformed. The mesh must therefore account for these deformations induced by the faults.
Case of a single fault
In the course of this step, the fault is taken into account within the regular grid, by carrying out the following steps:
A difficulty may arise due to the fact that degenerate quadrilaterals might be produced by moving this end. This problem is illustrated in
Consider three nodes such as:
If (i0=i1 or j0=j1) then (i0, j0) and (i1, j1) form an edge (two nodes moved successively form either a diagonal, or an edge of a quadrilateral). Likewise, (i1, j1) and (i2, j2) form a second edge if (i1=i2 or j1=j2). Consequently, these two constructed edges will be in the same quadrilateral if |i0−i2|=1 and |j0−j2|=1 and will be considered to be aligned if the angle between them exceeds a certain alignment threshold. This situation is not acceptable since the hexahedrons having this quadrilateral as a face will be twisted. To surmount this difficulty, before displacing the closest end, a check is carried out to verify that this displacement does not give rise to a degenerate quadrilateral. A quadrilateral is called degenerate if at least one of its angles is greater than a fixed angle threshold (if an angle of the quadrilateral is above this threshold, this signifies that another angle of this quadrilateral will be small and consequently the value of the quality will be low). If, in actual fact, displacing the closest end gives rise to such a quadrilateral, then the other end of the intersected edge is displaced. According to the example of
The procedure is implemented by a computer and the algorithm relating to this step is described. The algorithm commences with the identification of the quadrilateral which contains the first point of the fault. Next, the intersection between the first segment of the fault and the edges of this quadrilateral is computed. If no intersection is detected, the algorithm proceeds to the next segment of the fault and loops until an intersection is found. Thereafter, the end of the intersected edge closest to the point of intersection is moved toward the latter. If this results in the formation of a degenerate quadrilateral, it is the other end of the edge which is moved. To compute the next intersection, the algorithm advances along the quadrilaterals having the last node paired in common and excludes those which have the last two nodes paired like a diagonal or an edge. If there is no longer any intersection between the current segment and the quadrilaterals to be visited, the quadrilaterals to be visited are reassigned with that which contains the extremity of the current segment, and the algorithm proceeds to the next segment to compute the intersection, and so on and so forth. On completion, the two extremities of the fault are fitted with the closest nodes.
Case of multiple faults
In the presence of multiple faults, not all the nodes of the grid are free to move, since a node which has already been displaced to account for a first fault must not be modified to account for a second fault. The fitting algorithm must therefore account for this new more constraining context, and a refinement procedure can be performed to free oneself from this constraint.
According to the invention, multiple faults are taken into account within the regular grid, by carrying out the following steps:
The refinement of the mesh can be carried out by the following steps:
This refinement is also applied to all the horizons situated lower down, if the horizons are processed from the lowest to the highest.
The procedure is implemented by a computer with the algorithm relating to this step being described as follows. To describe whether a node is constrained by one or more faults, two notions are introduced: the first termed “Faulted degree of a node” and the second termed “Passing faults of a node”.
Faulted Degree of a Node (Fd):
A node S is said to be of faulted degree n if it is the point of intersection of n faults. If n is zero, then no fault passes through this node. The faulted degrees of the nodes are stored as a property of the grid, and are useful for the phase of taking the faults into account, to verify whether a node is constrained by faults, for the grid optimization phase (fd=0: move the node freely; fd=1: project onto the passing fault; fd>=2: do not move), and finally for the grid 3D mapping phase, so as to decide when and how a node is split.
Passing Faults of a Node:
This is the list of faults which pass through a node.
The global process of fitting the multiple faults is iterated from bottom to top on all the horizons, to generate regular grids and fit the faults thereto. For each horizon, a grid of the same dimension as those situated lower down (geologically older) are generated first, and then all the points of intersection between the faults and the ends of the faults are fitted. The global process iterates over all the faults to account for them in the mesh. When a refinement is necessary on a grid of a horizon to pair a fault, the grids of the horizons already processed undergo the same refinement in an exhaustive manner, to retain a single dimension for all the grids.
Concerning the structures of the data, having regard to the fact that the number of the nodes of a grid on a horizon may increase subsequent to a refinement, a chained list is used to represent the grid, to facilitate the addition of a new node at an arbitrary position. Hence, the nodes traversed by a fault are recorded in a linked list to allow fast insertion.
To fit a given fault, a segment-by-segment search is conducted for the points of intersection between the fault and the grid, and the nodes to be moved to these points of intersection are decided. A point of intersection is computed in the following manner. The first node in the list of nodes traversed by a fault which represents the start of the first segment of the fault is the beginning, the four neighboring cells of this node are recovered, and a search is conducted for the intersection of the first segment with these four cells. If an intersection is found, an end of the intersected edge is chosen and it is moved to the point of intersection. The neighboring cells of the node representing this intersection become the cells to be visited to compute the next intersection with the same fault segment, with the exclusion of those which have the new paired node and the one before as edge or diagonal. If there is no longer any intersection between the current segment and the cells to be visited, the cells to be visited are reassigned with that which contains the extremity of the current segment, and the next segment is processed to compute the intersection, and so on and so forth.
While scanning across the fault, already paired points may be encountered, which represent for example the intersections of the fault with other faults present on the same horizon. It is therefore necessary to verify, for the last paired node, whether there exists a previous node in the list of nodes traversed by a fault. If so, if the two nodes form an edge or a diagonal of a quadrilateral, the previous node becomes the current node for computing the cells to be visited, and computation proceeds to the fault segment which corresponds to this node to continue the intersection computation.
Thereafter, the refinement process can be triggered by the addition as a new node of a point of intersection between a fault and a grid.
Step f) optimization of the quality of mesh cells
With the aid of the grid optimization steps described hereinabove, a two-dimensional grid is generated and optimized in terms of quality of the mesh. For the application of these grid optimization steps, the rigid nodes correspond to the nodes of the boundaries of the grid to be optimized and to the nodes on the traces of the faults of the grid to be optimized. The displacement of these rigid nodes is therefore imposed on the nodes of the boundaries and the nodes of the traces of faults of the grid to be optimized.
Step g) transformation into three-dimensional gridded surface
This step maps the regular grid, which is deformed to account for the faults and optimized, into a 3D real horizon.
According to the invention, it is possible to transform the regular grid into a three-dimensional gridded surface by performing a change of reference frame, from the reference frame of the regular grid to the reference frame of the triangulated 3D surface, of the nodes of the grid. The coordinates of the off-fault nodes are determined on the basis of their barycentric coordinates in a reference frame defined by the triangle of the 2D triangulated surface to which they belong. The coordinates of the nodes situated on a fault are determined on the basis of their curvilinear abscissae on the fault. The connectivities are thereafter established to maximize the number of quadrilaterals with only the quadrilaterals for which a fault passes through one of its diagonals being divided into two triangles.
In particular it is possible to carry out the following steps:
This procedure is implemented by a computer with the algorithm relating to this procedure is described as follows. The principal structure represents a quasi-regular 3D grid, with certain nodes (I, J) split into two or more.
The principal structure is described by:
For a node (I, J), the following information is stored:
For a pseudo-quadrilateral cell, the following information is stored:
The algorithm commences by scanning through the nodes of the 2D regular grid paired with the faults. For each node (I, J), the passing faults are first recovered. If there are no passing faults, a triangle containing the node in the unfolded triangulation is located, and the corresponding node on the 3D real horizon is computed using the barycentric coordinates of the triangle of the same number on the 3D original triangulation. If the node (I, J) is situated on a single fault, the segment is first recovered for the fault where this node is situated in unfolded space. The curvilinear abscissa of the node on this segment is thereafter computed. Next, the two corresponding split nodes are computed using this curvilinear abscissa on the same segment on two sides of the fault in real space. For a node situated at the intersection of the faults, its number in the triangulation is recovered and the references on the sides of faults which intersect thereat are recorded.
Apart from the coordinates, the information relating to the faults or the containing triangle is stored to facilitate the generation of the volume mesh by matching the horizons.
Once the nodes of the unfolded grid have been mapped into a real horizon, the connectivity between these nodes is established by forming pseudo-quadrilateral cells. These cells are constructed one by one by scanning first the direction I and then J. The type of the current cell is first determined. To do this, it is necessary just to verify whether there is a passing fault on one of two diagonals of the quadrilateral in the unfolded grid. If so, the type is 1 or 2; otherwise, the type is 0. The number of vertices for the current cell is then obtained as a function of its type, as well as the index (I, J) of each vertex. With this index, a test is carried out to verify whether a vertex is or is not split, by examining the list of nodes (I, J) which have been previously filled. If all the vertices are split, the current cell is added to a temporary list, which saves the cells with all its split vertices for processing at the end, and processing passes to the next cell. Otherwise, the non-split node IJ is immediately mapped, since there is a single corresponding 3D node. For each already mapped node of the current cell, the adjacent and opposite nodes split in the cell are computed in the following manner: If the non-split node and the adjacent one is on the same fault, the node on the same side of the fault is chosen; otherwise, the hundredth point on the edge of the non-split node and the adjacent one is taken. It is located in the unfolded triangulation, and it is mapped into a real triangulation by using the same triangle and the same barycentric coordinates. The distance is thereafter computed between the mapped point and all the split nodes of the adjacent node IJ. The closest node is chosen as the corresponding vertex of the current cell. For the node opposite the current node (already mapped) of the diagonal, the split node on the other side of the fault is chosen. This procedure loops back until all the vertices are mapped.
For a cell with all its vertices split, for each of its edges, the facing edge in the neighboring cell is tagged, the reference of the fault where this facing edge is situated is recovered, and the current edge is therefore on the other side of the fault with a different reference. The corresponding split nodes are thus determined according to this reference.
Step h) generation of the mesh of the subterranean medium
Steps a) to f) are carried out on each horizon while preserving the same number of quadrilateral for each horizon. It is therefore possible to reconstruct a three-dimensional mesh between the various horizons. This step constructs a three-dimensional mesh of the medium, by directly linking the nodes having the same coordinates I, J on two neighboring horizons to form three-dimensional mesh cells, and in optionally cutting these mesh cells if they are traversed by the fault. If a fault intersects this link, then a node is linked with the fault by considering a direction of a neighboring node. This is illustrated in
When a quadrilateral on one of two horizons is split into two triangles via a fault which passes through the diagonal, the corresponding quadrilateral on the other horizon is also cut virtually in the same manner and two prisms will be formed. Nevertheless, one of the two prisms is twisted and the interface of the fault is not complied with. In
A problematic configuration is produced when the horizons on the two sides of the fault slip on the latter, as illustrated by
The malformation of the elements corresponding to the two configurations hereinabove is due to an incorrect choice of the direction of tying between two horizons. The solution thus resides in the detection and the correction of wrong directions. To detect the presence of a wrong direction during the creation of a 3D element, a test is first carried out to verify whether there exists a segment linking the same I, J which passes through a fault. If so, the 3D element actually comprises directions to be corrected and it must be replaced with two 3D elements which follow the corrected directions. Otherwise, a test is carried out to verify whether there exists a fault which passes through two edges of the 3D element, one on the horizon at the top and another on that at the bottom. If so and if the local indices of these two edges are different, the element also has to be corrected. Let E be the 3D element in question, F be the traversing fault, E1 and E2 be the two 3D elements to be created to replace E. The correct directions are computed in the following manner:
By following the corrected directions in E2, all the as yet unprocessed nodes in H0 are projected onto the fault surface and are tied to its corresponding images. All the unprocessed nodes in H1 are projected onto the fault surface following the corrected directions in E1, as illustrated by
A mesh constructed according to the invention is particularly suitable for simulating the flows within a subterranean medium in a zone with complex geometry. The invention thus provides an accurate tool for carrying out a numerical representation of the subterranean medium (for example basin modeling). This representation of the medium allows in particular reservoir simulation, or simulations of injection of CO2 into a subterranean formation. A procedure for simulating the flows within a subterranean medium in a zone with complex geometry comprises the generation of a hexa-dominant mesh on the basis of the procedure according to the invention, and then the carrying out of simulations by suitable software (flow simulator, basin simulator or reservoir simulator) underpinned by the generated mesh.
These simulations allow testing several production schemes and then to optimizing oil field exploration, geological reservoir exploration or exploitation, or the injection of gas into subterranean media by carrying out the optimal production scheme (for example according to a criterion of profitability, volume of oil recovered, etc.). A simulated production scheme can correspond to the choice of the location of a new well (producer or injector), to the choice of tools (for drilling, exploration, etc.), to the choice of fluids used or recovered, to the choice of exploitation conditions (injection flowrate, etc.). The exploitation of the subterranean medium then carries out the choice that was made, such as, for example, the drilling of a new well, the injection of fluid, the modification of the exploitation and/or drilling conditions, etc.
Annex: Computation and Minimization of the Gradient of the Displacement Field
In order to better explain the procedure for computing and minimizing the gradient of the displacement field of the nodes of the grid, the following notations and definitions are introduced:
In this case, it may be written:
$$
where C_{n }is a connected component.
$$
It may be assumed that the surface Ω_{0 }is connected. If this surface is not connected, it suffices to perform an optimization on each of its connected components C_{n}.
Optimizing the mesh of the surface Ω_{0 }minimizes the deformation of the displacement field from the surface Ω_{R }to the surface Ω_{0}. Therefore, the following optimization problem should be solved:
$$
where:
$$
Denote by T″ the triangulation of the reference surface Ω_{R}. It may therefore be written:
$$
where T is a triangle of Ω_{R }of area equal to ST.
Denote by N_{i}, N_{j }and N_{k }the three nodes of the triangle T and let M be an arbitrary point inside the triangle T. Therefore, it may be written:
{right arrow over (OM)}=α{right arrow over (ON)}_{i}+β{right arrow over (ON)}_{j}+γ{right arrow over (ON)}_{k} (3)
With:
$$
The displacement field {right arrow over (ν)}(M)=(ν_{1}(M),ν_{2}(M)) at the point M is such that:
ν_{m}(M)=αu_{m}(N_{i})+βu_{m}(N_{j})+γu_{m}(N_{k})m=1,2 (5)
Where {right arrow over (u)}(N_{t})=(u_{1}(N_{t}),u_{2}(N_{t}))t=i, j, k represents the displacement field at the node N_{t}. Thus, the gradient of this displacement field is obtained:
∇ν_{m}(M)=∇(α)u_{m}(N_{i})+∇(β)u_{m}(N_{j})+∇(γ)u_{m}(N_{k}) (6)
From this the basis of equation (4) is deduced as:
$$
According to (7), ∇ν_{m}(M) is independent of x and y, therefore this gradient is constant over T. It can then be deduced:
$$
We thus obtain a quadratic form on the triangle T that may be written:
$$
where for t=i, j, k u_{m}′=u_{m}(N_{t}) and
$$
Ultimately, it is obtained from equation (9) and by definition of the adjacent mesh cells:
$$
where a summation over each node i of the triangulation and then over each node j connected to the node i by the edge [N_{i},N_{j}] is performed. This connection is denoted by j↔i as is illustrated in
$$
Thus, the function to be minimized is a quadratic form whose variables are the displacement fields u_{m}^{i }of each node N_{i }of the reference surface. It is denoted by:
$$
where Nb is the number of nodes of the reference surface. This function is a minimum if its gradient is zero, that is:
$$
The minimization therefore amounts to solving the following system:
$$
that may also be written in the form:
$$
with:
$$
According to the invention, the approach involves three steps:
1. The first solves the following system by imposing the displacement field solely on the boundary of Ω_{R}:
$$
that may be written in the form:
L^{(1)}U=B (17)
where L^{(1) }is a square matrix of order N equal to the number of nodes not belonging to the boundary of Ω_{R }such that:
$$
with the right-hand side B of size N such that:
$$
2. The second step fits the rigid nodes other than those of the boundaries (for example the traces of faults) on the optimized mesh resulting from step 1.
3. After fitting the fault traces on the mesh, the third step solves the following system by imposing the displacement field not only on the boundary but also on the other rigid nodes (for example the traces of faults):
$$
that may be written in the form:
L^{(2)}U=T (21)
where L^{(2) }is a square matrix of order M equal to the number of nodes not belonging to the boundary and to the rigid nodes (for example the traces of faults) such that:
$$
with the right-hand side T of size M such that:
$$
The right-hand sides B and T imply that the displacement field is imposed on the boundary in step 1 and on the boundary and the traces of faults in step 3. The matrices L^{(1) }and L^{(2) }occurring in the systems (17) and (21) are two sparse symmetric matrices since, in the ith row, all the coefficients are zero except for those j such that j is connected to i by an edge. Thus, band storage can be used for the storage of these two matrices. The adaptive Gauss Seidel procedure can be used to solve the two systems.
It is generally possible to find formula (11) for deformations in 2D or 3D and also for conformal unfoldings. The articles cited hereinbelow use them for this purpose. There are also other ways of approximating the Laplacian such as taking for example
$$
where N_{i }is the number of neighbors of node i.
Oscar Kin-Chung Au, Chiew-Lan Tai, Ligang Liu, and Hongbo Fu, “Dual Laplacian Editing for Meshes”, IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. 12, NO. 3, pages 386-395, Hangzhou-China, MAY/JUNE 2006.
H. Masuda*, Y. Yoshioka, Y. Furukawa, “Interactive Mesh Deformation Using Equality-Constrained Least Squares”, Tokyo, Japan, May 2006.
Mathieu Desbrun, Mark Meyer, Pierre Alliez, “Intrinsic Parameterizations of Surface Meshes”, EUROGRAPHICS 2002, Volume 21, Number 2, Oxford UK, 2002.
Yanzhen Wang*, Kai Xu, Yueshan Xiong and Zhi-Quan Cheng, “2D shape deformation based on rigid square matching”, COMPUTER ANIMATION AND VIRTUAL WORLDS, pages 411-420, Changsha-China, August 2008.
The algorithm can be tested on triangulated surfaces. A good distribution of the mesh cells after fitting is noted for an optimization performed with the Matlab® (Mathworks, USA) computation software. The reference horizon surface used is a Coons mesh generated starting from the boundary. The improvement in quality makes itself felt, especially when far from the boundary. A posteriori, the two parts separated by the traces of faults considered at the start to be connected components do not actually have the same uniformity of meshing. No folding is detected for this case because the deformations are not sufficiently imposing.
Number | Date | Country | Kind |
---|---|---|---|
12 01451 | May 2012 | FR | national |
Number | Name | Date | Kind |
---|---|---|---|
8504300 | Dorn | Aug 2013 | B2 |
8711140 | Mallet | Apr 2014 | B1 |
8743115 | Mallet | Jun 2014 | B1 |
8965745 | Lepage | Feb 2015 | B2 |
9959662 | Mitchell | May 2018 | B2 |
20010007451 | Aono | Jul 2001 | A1 |
20040194051 | Croft | Sep 2004 | A1 |
20060267978 | Litke | Nov 2006 | A1 |
20080243454 | Mallet | Oct 2008 | A1 |
20100017181 | Mouton | Jan 2010 | A1 |
20110015910 | Ran | Jan 2011 | A1 |
20110106507 | Lepage | May 2011 | A1 |
20120022837 | Asbury | Jan 2012 | A1 |
20130124161 | Poudret | May 2013 | A1 |
20130204598 | Mallet | Aug 2013 | A1 |
20130218539 | Souche | Aug 2013 | A1 |
20160125555 | Branets | May 2016 | A1 |
20160245951 | Kartasheva | Aug 2016 | A1 |
20180031721 | Etiene Queiroz | Feb 2018 | A1 |
20180232950 | Brewer | Aug 2018 | A1 |
Number | Date | Country |
---|---|---|
2 948 215 | Jan 2011 | FR |
Entry |
---|
Bennis et al, 3D line-support grid flattening for more accurate geostatistical reservoir population with petrophysical properties, Engineering with Computers (2014) 30:403-421 (published online Jan. 24, 2013). |
Horna et al, Isometric Unfolding of Stratigraphic Grid Units for Accurate Property Populating: Mathematical Concepts, EAGE, European Conference on the Mathematics of Oil Recovery, Sep. 2010, Oxford, United Kingdom, 14 pages. |
Frantz Maerten, Chaper 11 of Geomechanics to solve geological structure issues: forward, inverse and restoration modeling, Geophysics, Universite Montpellier II Sciences et Techniques du Languedoc, pp. 331-357, 2010. |
Jean-François Rainaud, Vincent Clochard, Thomas Crabié, and Houman Borouchaki, Using a ChronoStratigraphic Unfolding Workflow to build an a priori model for Stratigraphic Inversion with accurate Horizon and Fault Fitting, SEG Technical Program Expanded Abstracts 2015: pp. 1927-1931. |
Ran, Longmin, Hex-dominant mesh generation for subterranean formation modeling, Engineering with Computers, 28, pp. 255-268, 2011. |
Brahim Yahiaoui, Houman Borouchaki, Abdallah Benali, Hex-Dominant Mesh Improving Quality to Tracking Hydrocarbons in Dynamic Basins, Oil & Gas Science and Technology—Revue d'IFP Energies nouvelles, Institut Francais du Petrole, 2014, 69 (4), pp. 565-572. |
Brahim Yahiaoui, Section 5.2 of Maillage dynamique tridimensionnel pour la simulation de l'ecoulement dans un bassin sedimentaire, Dissertation, pp. 103-112, Dec. 17, 2013. |
Poudret, Mathieu et al, A Volume Flattening Methodology for Geostatistical Properties Estimation, Proceedings of the 20th International Meshing Roundtable, Springer, pp. 569-585, Jan. 2012. |
Mello et al. Tulsa: AAPG/Datapages Discovery Series 7 (2003): 27 pages. Obtained from https://www.researchgate.net/profile/Paulo_Cavalcanti/publication/228592962_A_Topologically-based_Framework_for_Three-dimensional_Basin_Modeling/links/0912f50d57a58188dd000000.pdf on May 10, 2018 (Year: 2003). |
D'Azevedo, Eduardo F., and R. Bruce Simpson. “On optimal triangular meshes for minimizing the gradient error.” Numerische Mathennatik 59, No. 1 (1991): 321-348. (Year: 1991). |
Canann, et al. “An Approach to Combined Laplacian and Optimization-Based Smoothing for Triangular, Quadrilateral, and Quad-Dominant Meshes.” In IMR, 16 pages. 1998. Obtained from http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.88.9616&rep=rep1&type=pdf on May 10, 2018 (Year: 1998). |
Oscar Kin-Chung Au et al: “Dual Laplacian Editing for Meshes”, IEEE Transactions on Visualization and Computer Graphics, vol. 12, No. 3, pp. 386-395, Hangzhou-China, May/Jun. 2006. |
H. Masuda et al: “Interactive Mesh Deformation Using Equality-Constrained Least Squares”, Tokyo, Japan, May 2006. (24 pages). |
Mathieu Desbrun et al: “Intrinsic Parameterizations of Surface Meshes”, Eurographics 2002, vol. 21, No. 2, Oxford UK, 2002. (10 pages). |
Yanzhen Wang et al: “2D Shape Deformation Based on Rigid Square Matching”, Computer Animation and Virtual Worlds, pp. 411-420, Changsha-China, Aug. 2008. |
Number | Date | Country | |
---|---|---|---|
20140052427 A1 | Feb 2014 | US |