METHOD OF MODELLING A SEDIMENTARY BASIN USING A HEX-DOMINANT MESH REPRESENTATION

Abstract

The present invention relates to a method of modelling a sedimentary basin by means of a numerical basin simulation solving at least a balance equation of poromechanics according to a face-based smoothed finite-element method for determining at least a stress field and a deformation field. The method according to the invention notably comprises the following steps: subdividing the hexahedral cells of a mesh representation of a state of the basin into at least eight hexahedral subcells, determining a transition relation between the degrees of freedom of the nodes of the subcells and the degrees of freedom of the nodes of the cell to which the subcells belong, and determining a stiffness and nodal forces from at least this transition relation and a strain-displacement relation determined for the subcells.

Claims

1. A computer-implemented method of modelling a sedimentary basin, the sedimentary basin having undergone a plurality of geological events defining a sequence of states of the basin, by means of a computer-executed numerical basin simulation, the numerical basin simulation solving at least one balance equation of poromechanics according to a face-based smoothed finite-element method for determining at least a stress field and a strain field, characterized in that the method comprises carrying out at least the following steps: A. performing physical quantity measurements relative to the basin by means of sensors and constructing a mesh representative of the basin for each of the states of the basin, the meshes representative of the basin for each of the states predominantly consisting of hexahedral cells, B. by means of the numerical basin simulation and of the meshes for each of the states, determining at least a strain field and a stress field for each of the states by carrying out at least the following steps for each of the meshes representative of the states: a) subdividing each of the hexahedral cells of the mesh into at least eight hexahedral subcells, b) for each face of each of the subcells, determining a smoothing domain according to the face-based smoothed finite-element method applied to the subcells, c) for each of the smoothing domains, determining a strain-displacement relation according to the face-based smoothed finite-element method applied to the smoothing domains, d) for each of the smoothing domains, determining a transition relation between the degrees of freedom of the nodes of the subcells containing the smoothing domain and the degrees of freedom of the nodes of the at least one cell to which the subcells belong, e) determining a stiffness and nodal forces for each of the smoothing domains from at least the transition relation determined for the smoothing domain and from the strain-displacement relation determined for the smoothing domain, f) determining a stiffness and nodal forces relative to the mesh from at least the stiffness and the nodal forces determined for each of the smoothing domains, g) modelling the sedimentary basin by determining at least the displacement field and the stress field for the mesh, by means of the numerical basin simulation and at least the stiffness and the nodal forces determined for the mesh.

2. A method as claimed in claim 1, wherein each of the hexahedral cells of the mesh is subdivided into eight hexahedral subcells, the subdivision of one of the cells being performed in such a way that each face of the cell consists of four faces of four subcells among the eight subcells.

3. A method as claimed in claim 1 wherein, in step b), the smoothing domain relative to one of the faces belonging to at least one of the subcells is determined by connecting each node of the face with at least one point located at the barycenter of the subcell(s) to which the face belongs.

4. A method as claimed in claim 3, wherein the strain-displacement relation is a transformation matrix relating a displacement vector to a strain vector according to a formula of the type: $.Math.B=\left[B_1.Math.B_2.Math.\mathrm{.Math.}.Math..Math.B_{n-1}.Math.B_n\right]$ $.Math.\mathrm{with}.Math..Math.B_i=\frac{1}{\text{?}}.Math.\int \text{?}.Math.{\mathrm{nN}}_i.Math.\mathrm{dT}$ $\text{?}.Math.\text{indicates text missing or illegible when filed}$ where B.sub.i is a matrix respectively associated with node i, with i=1 . . . n, of the subcell(s) to which the smoothing domain Ω.sub.k.sup.s belongs, A.sub.k.sup.s is the volume of the smoothing domain Ω.sub.k.sup.s, ∂Ω.sub.k.sup.s is a boundary of the smoothing domain, n is an outward normal of the smoothing domain Ω.sub.k.sup.s and N.sub.i is a shape function associated with node i.

5. A method as claimed in claim 4, wherein the transition relation for a smoothing domain relative to one of the faces of at least one of the subcells is a transition matrix determined according to at least the shape function matrices of the finite element assigned to the subcell(s) to which the smoothing domain belongs, the shape function matrices being evaluated at each node of the subcell(s) to which the smoothing domain belongs.

6. A method as claimed in claim 5, wherein the stiffness and the nodal forces for each of the smoothing domains take the form of a stiffness matrix and a nodal force vector respectively, and the stiffness matrix K.sub.Ωvirt and the nodal force vector f.sub.Ωvirt for one of the smoothing domains Ω.sub.virt are respectively determined according to formulas of the type:
K.sub.Ω.sub.virt=∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.TDB.sub.Ω.sub.virtT.sub.Ω.sub.virtdV
f.sub.Ω.sub.virt=−∫.sub.Ω.sub.virtN.sup.TbdV−∫.sub.BΩ.sub.virtN.sup.TtdA−∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.TDε.sub.0dV+∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.Tσ.sub.0dV wherein T.sub.Ωvirt is the transition matrix for the smoothing domain Ω.sub.virt, B.sub.Ωvirt is the transformation matrix for the smoothing domain Ω.sub.virt, N is a matrix of the shape functions, D is a matrix representative of the material stiffness, b is a body force vector, t is a surface force vector, ε.sub.0 is an initial strain and σ.sub.0 is an initial residual stress.

7. A method as claimed in claim 6 wherein, in step f), and if all of the cells of the mesh are hexahedral, the stiffness and the nodal forces relative to the mesh are determined according to respective formulas of the type:
K=Σ.sub.i=1.sup.nvK.sub.Ω.sub.virt t
f=Σ.sub.i=1.sup.nvf.sub.Ω.sub.virt t where nv is the total number of the smoothing domains.

8. A computer program product downloadable from a communication network and/or recorded on a computer-readable medium and/or processor executable, comprising program code instructions for implementing the method as claimed in claim 1, when the program is executed on a computer.

9. A method for exploiting hydrocarbons present in a sedimentary basin, the method comprising at least implementing the method for modelling the basin as claimed in claim 1, and wherein, from at least the modelling of the sedimentary basin, an exploitation scheme is determined for the basin, comprising at least one site for at least one injection well and/or at least one production well, and the hydrocarbons of the basin are exploited at least by drilling the wells of the site and by providing them with exploitation infrastructures.

Description

BRIEF DESCRIPTION OF THE FIGURES

[0053] 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:

[0054] FIG. 1 shows an illustrative example of a sedimentary basin,

[0055] FIG. 2 shows an example of a sedimentary basin (left) and an example of a mesh representation (right) of this basin,

[0056] FIG. 3 shows an example of structural reconstruction of a sedimentary basin according to an embodiment of the invention, represented by three strain states taken in three different geological states,

[0057] FIG. 4 shows the strain in a beam subjected to bending, calculated with the FEM method and the FS-FEM method, in the case of flattened hexahedral cells,

[0058] FIG. 5 shows an example of subdivision of a hexahedral cell into eight hexahedral subcells, and

[0059] FIG. 6 illustrates the construction of a smoothing domain according to the invention for a face belonging to two subcells.

DETAILED DESCRIPTION OF THE INVENTION

[0060] According to a first aspect, the invention relates to a computer-implemented method of modelling a sedimentary basin, by means at least of a numerical basin simulation solving at least a balance equation of poromechanics according to a face-based smoothed finite-element method known as FS-FEM method.

[0061] In general terms, the balance equation of poromechanics is written (see for example document (Coussy, 1991)):

∇.Math.σ′+(ρ+m)g=∇.Math.(pB)

where σ′ is the effective stress tensor, ρ the homogenized density of the porous medium, m the fluid mass exchange, g the acceleration of gravity, p the fluid pressure and B the Biot tensor. In a numerical simulation such as a basin simulation, this equation is numerically solved. In general, this equation is usually discretized in a mesh representation (or mesh) by means of the finite-element method (FEM) or one of its variants such as, for example, the smoothed finite-element method (S-FEM).

[0062] Different variations of the S-FEM method exist in the literature in the field of mechanical simulation, each with its convergence and accuracy properties. The variation known as FS-FEM (Face-based Smoothed Finite-Element Method) is known as the most accurate variation of the S-FEM method. However, in the literature, the FS-FEM method is described to date only for tetrahedral cells, and all the application cases found in the literature in the field of mechanical simulation use a tetrahedral mesh only.

[0063] The applicant has observed that the FS-FEM method is in practice unstable in the case of a mesh having very flattened hexahedral cells. It therefore appears that the FS-FEM method cannot be applied as is in the case of basin models whose mesh representation comprises cells representative of thin geological layers (for example of the order of a few meters in the vertical direction) and/or when these layers comprise stratigraphic pinchouts (characterized by an angle less than 1°).

[0064] The instability of the FS-FEM method in the presence of flattened hexahedral cells (i.e. whose ratio of the greatest horizontal edge to the smallest vertical edge is at least 5, preferably 10) materializes in two different ways in numerical simulations: either the computation stops with a fatal error because the stiffness matrix of the model is singular (zero determinant), which means that it cannot be inverted, or the computed solution is incorrect and it locally exhibits oscillations and/or peaks in the computed displacements, stresses and strains. An example of the latter case is shown in FIG. 4, within the context of a beam subjected to bending, the beam being represented by cells consisting of flattened hexahedra (aspect ratio 1:10). More precisely, FIG. 4, top, shows the geometry of the deformed beam computed with the FEM method. This deformed shape is correct because it is coherent with the known analytical solution to the problem. FIG. 4, bottom, shows the deformed beam geometry obtained with the FS-FEM method: the strain is clearly incorrect and the accordion shape is the typical sign of numerical instability (see also paragraph 4.6 of document (Cook, 1995)).

[0065] The present invention relates to an improvement to the FS-FEM method in order to make it applicable to the hexahedral cells used in basin models. More specifically, in a particularly shrewd manner, the applicant has shown that improving the FS-FEM method according to the invention allows to stabilize the numerical simulation, including in the case of flattened hexahedral cells.

[0066] Thus, the numerical basin simulation according to the invention aims to solve numerically at least the balance equation of poromechanics according to the face-based smoothed finite-element method (FS-FEM) in the case of a hex-dominant mesh representation.

[0067] In general terms, discretization of the balance equation of poromechanics leads to write a system of algebraic equations of the following form (see for example Equation 2.23 of the document (Zienkiewicz and Taylor, 2000)):

Ka+f=r  (1)

where K is the stiffness matrix of the system, a is the vector of the degrees of freedom of the nodes of the mesh used for discretization (nodal displacements, see Equation 2.1 of the document (Zienkiewicz and Taylor, 2000)), f is a vector of the nodal forces (equivalent nodal forces, see Equation 2.24B of the document (Zienkiewicz and Taylor, 2000)), and r is the vector of the external concentrated forces (see Equation 2.14 of the document (Zienkiewicz and Taylor, 2000)). Details of the solution of this system of equations by means of a face-based smoothed finite-element method specifically suited to the case of hexahedral cells are given in step 2 of the method according to the invention.

[0068] According to the invention, it is assumed that the sedimentary basin has undergone a plurality of geological events defining a sequence of states of the basin, each of said states extending between two successive geological events. Preferably, the sequence of states can cover a period of time covering at least the production of hydrocarbons, notably by maturation of an organic matter present in a source rock of the basin, deformation of the basin due to 3D mechanical stresses, displacement of these produced hydrocarbons into at least one reservoir rock of the basin over geological times. By way of non-limitative example hereafter, Ai denotes a state of the sequence of states of the basin, i being an integer ranging from 1 to n, An representing the state of the basin at the current time. According to the invention, n is at least 2. In other words, the sequence of states according to the invention comprises the state of the basin at the current time and at least one state of said basin at an earlier geological time.

[0069] According to a second aspect, the invention relates to a method of exploiting hydrocarbons present in a sedimentary basin, the method according to the second aspect comprising implementing the method of modelling a sedimentary basin according to the first aspect of the invention.

[0070] The method according to the first aspect of the invention comprises at least steps 1) and 2) described hereafter.

[0071] The method according to the second aspect of the invention comprises at least steps 1) to 3) described hereafter.

1) Constructing Mesh Representations for Each of the States

[0072] This step consists in measuring physical quantities relative to the basin by means of sensors and in constructing a mesh representative of the basin at the current time, then meshes for each state of the basin. According to the invention, the mesh representations of the basin for each of said states predominantly consist of hexahedral cells. This is conventionally known as hex-dominant meshing.

[0073] During this step, at least the three following substeps are carried out:

[0074] 1.1) Measuring Physical Quantities Relative to the Basin

[0075] This substep consists in acquiring physical quantity measurements relative to the basin studied, by means of sensors.

[0076] By way of non-limitative example, the sensors can be logging tools, seismic sources and receivers, fluid samplers and analyzers, etc.

[0077] Thus, the measurements according to the invention can consist of outcrop surveys, seismic acquisition surveys, measurements in wells (logging for example), petrophysical and/or geochemical analyses of core samples taken in situ.

[0078] These measurements allow to deduce petrophysical properties associated with the basin studied, such as facies (lithology), porosity, permeability, or the organic matter content at measuring points of the basin. Information relative to the properties of the fluids present in the basin can also be obtained, such as saturation values of the various fluids present in the basin. Temperatures can also be measured at different points of the basin (notably bottomhole temperatures).

[0079] 1.2) Constructing a Mesh Representative of the Basin at the Current Time

[0080] This substep consists in constructing a hex-dominant mesh representative of the basin at the current time, from the physical quantity measurements performed in the previous substep.

[0081] More precisely, construction of a mesh representation of a basin consists in discretizing in three dimensions the architecture of the basin and in assigning properties to each of the cells of this mesh. The physical quantity measurements performed at various points of the basin as described above are therefore notably exploited, extrapolated and/or interpolated, in the various cells of the mesh, according to more or less restrictive hypotheses.

[0082] Most often, the spatial discretization of a sedimentary basin is organized in cell layers representing each the various geological layers of the basin studied. FIG. 2 illustrates, on the left side, an example of a sedimentary basin at the present time and, on the right side, an example of a mesh of this basin.

[0083] According to an implementation of the invention, the mesh constructed for the current state An of the basin studied notably comprises in each cell information on the lithology, a porosity value, a permeability value, an organic matter content, and properties relative to the fluids present in the cell, such as saturation.

[0084] According to an implementation of the invention, the mesh constructed for the current state An of the basin predominantly consists of hexahedral cells. The OpenFlow® software (IFP Energies nouvelles, France) or the GOCAD® software (Emerson-Paradigm, USA) can be used.

[0085] 1.3) Structural Reconstruction of the Basin Architecture for the Various States

[0086] This substep consists in reconstructing the past architectures of the basin for the various states Ai, with i ranging from 1 to n−1. The mesh constructed in the previous substep, which represents the basin at the current time, is therefore deformed in order to represent the anti-chronological evolution of the subsoil architecture over geological times, and for the various states Ai. At the end of this substep, a mesh is available for each state Ai, with i ranging from 1 to n.

[0087] According to a first embodiment of the present invention, structural reconstruction can be particularly simple if it is based on the assumption that its deformation only results from a combination of vertical movements by compaction of the sediment or by uplift or downwarping of its basement. This technique, known as backstripping, is described in (Steckler and Watts, 1978) for example.

[0088] According to another embodiment of the present invention, in the case of basins whose tectonic history is complex, notably in the cases of basins with faults, less restrictive hypotheses, such as structural restoration, need to be used. Such a structural restoration is described for example in document FR-2,930,350 A1 (US-2009/0,265,152 A1). Structural restoration consists in computing the successive deformations undergone by the basin by integrating the deformations due to compaction and those resulting from tectonic forces.

[0089] In the example of FIG. 3, three states are used to represent the subsoil deformation over geological times. The mesh on the left represents the current state, where a slip interface (a fault here) can be observed. The mesh on the right represents the same sedimentary basin for a state Ai, prior to the current state. For this state Ai, the sedimentary layers are not fractured yet. The central mesh represents an intermediate state, i.e. the sedimentary basin in a state Ai′ between state Ai and the current state. It is observed that the slip along the fault has started to modify the basin architecture.

2) Determining a Strain and Stress Field by Numerical Basin Simulation for Each State

[0090] In this step, by means of numerical basin simulation and of the meshes determined for each of the states in the previous step, at least one strain field and one stress field are determined for each state, by carrying out at least substeps 2.1 to 2.5 and substeps 2.7 to 2.8 described below for each of the states. In other words, at least substeps 2.1 to 2.5 and substeps 2.7 to 2.8 described below are applied for each mesh representative of a state of the basin studied.

[0091] According to the invention, substeps 2.1 to 2.5 are only applied to the hexahedral cells of the meshes of each state. According to an implementation of the invention, when the mesh representing any state comprises both hexahedral cells and cells of another type, substep 2.6 can be applied to the cells of another type than the hexahedral type. According to the invention, substeps 2.7 and 2.8 are applicable whatever the cell type of the mesh representative of a state.

[0092] 2.1) Subdividing the Hexahedral Cells of the Mesh

[0093] This substep consists in subdividing each of the hexahedral cells (also referred to as parent cells hereafter) of the mesh considered into at least eight hexahedral subcells (also referred to as virtual cells hereafter). It is clear that the non-hexahedral cells of the mesh considered are not subdivided as described in this substep. Processing of the other cell types of the mesh than the hexahedral type is detailed hereafter in substep 2.6.

[0094] According to a preferred embodiment of the invention, each hexahedral cell of the mesh can be subdivided into eight hexahedral subcells, subdivision being performed in such a way that each face of the parent cell consists of four faces of four subcells among the eight subcells.

[0095] According to this preferred embodiment of the invention, subdividing a hexahedral cell into eight hexahedral subcells can be performed as follows: [0096] dividing each edge of the parent cell into two edges of any length. The point in common between the two edges thus formed is referred to as virtual node and the other end of each edge, which coincides with one of the corners (real node) of the parent cell, is referred to as real node. An illustrative example of subdividing edges of a hexahedral cell is shown in FIG. 5, top left, where a black point represents a real node and a grey point represents a virtual node, [0097] subdividing each face of the parent cell into quadrangles, from the virtual nodes created by subdividing the edges and by adding an additional virtual node inside the face considered (or, in other words, not on an edge). An illustrative example of subdividing a face of a hexahedral cell into quadrangles is shown in FIG. 5, top right, where a quadrangle is formed by a part (between a real node and a virtual node) of two edges of the parent cell and by two edges (represented by a dotted grey line between two virtual nodes), one of these virtual nodes being located inside the face, [0098] from the quadrangles thus created and by adding an additional virtual node inside the cell (or, in other words, not on a face), from a given quadrangle (also referred to as base face), constructing a subcell (or virtual cell) using the quadrangle selected, the adjacent quadrangles that are not coplanar with the base face and the faces of three other subcells (or virtual cells). An illustrative example of subcell construction is shown in FIG. 5, bottom left. The subcells are thus created on the basis of both real and virtual nodes. It is observed that each face of the parent cell consists of four faces of four subcells.

[0099] FIG. 5, bottom right, shows the final result of the subdivision of a parent cell into eight subcells (or virtual cells) according to this preferred embodiment of the invention. The subcells are represented here in a disjoint manner to improve the clarity of the figure, it is however clear that two neighbouring subcells share a common face.

[0100] This substep is repeated for each hexahedral cell of the mesh representative of a state.

[0101] 2.2) Determining a Virtual Smoothing Domain for Each Face of Each Subcell

[0102] In this substep, for each face of each of the subcells as determined in the previous substep, a smoothing domain is determined according to the face-based smoothed finite-element method applied to the subcells. These domains are hereafter referred to as “virtual smoothing domains” (and denoted by Ω.sub.virt) because they are constructed from the virtual cells alone. “Applied to the subcells” means that the face-based smoothed finite-element method is applied as described in general terms in the literature, but by considering the hexahedral subcells as determined above, whereas the smoothed finite-element method is usually described for an application to cells. The difference between a subcell and a cell is not geometrical: both are hexahedral cells, but the nodes defining a subcell comprise both virtual and real nodes, whereas those defining a cell are real nodes exclusively. The difference between a real node and a virtual node is that the degrees of freedom of the problem (vector a in Equation (1) above) are associated only with the real nodes. No degree of freedom being associated with a virtual node, the virtual smoothing domains according to the invention require additional processing in relation to what is done for a smoothing domain of a real cell (parent cell in this case), and therefore in relation to what is described in the literature for the smoothed finite-element method. This additional processing is described in substep 2.4 below.

[0103] In general terms, the face-based smoothed finite-element method (FS-FEM) and, more generally, the smoothed finite-element method (S-FEM) use, in the mechanical balance formulation, a smoothed strain calculated by a weighted average of the strain according to an equation of the type (see for example Equation 4.19 of document (Liu and Nguyen Thoi Trung, 2010)):

[00002]$\begin{array}{cc}\stackrel{\_}{\mathrm{.Math.}}\left(x_k\right)=\frac{1}{A_k^s}.Math.\underset{{\Omega }_k^s}{\int }.Math.\widetilde{\mathrm{.Math.}}\left(x\right).Math.d.Math..Math.\Omega ,& \left(2\right)\end{array}$

where ε is the smoothed strain of the S-FEM method or of the FS-FEM method, {tilde over (ε)} is the strain (used by the FEM method) and A.sub.k.sup.s is the volume of a smoothing domain Ω.sub.k.sup.s at on which the strain is smoothed.

[0104] According to these methods, whose description can be found in document (Liu and Nguyen Thoi Trung, 2010), Equation (2) above is calculated in each smoothing domain Ω.sub.k.sup.s. The smoothing domains for the S-FEM method, the FS-FEM method or the smoothing domains according to the invention are determined in such a way that: [0105] the smoothing domains cover the mesh entirely, [0106] there is no overlap between the domains.

[0107] According to an implementation of the invention, at least one smoothing domain associated with a face belonging to one or two subcells is constructed by connecting each node (real or virtual) of the face considered with at least one point located at the barycenter of the subcell(s) to which the face belongs.

[0108] FIG. 6 illustrates the construction of such a smoothing domain for a face belonging to two subcells. Thus, the smoothing domain associated with face ABCD, which is common to two subcells, is constructed by connecting each node A, B, C, D of the face with the barycenter G and I of each subcell. The resulting smoothing domain, in grey in the figure, consists of two pyramids A BCDG and ABCDH.

[0109] It is clear that, for a face of a hexahedral subcell located on the borders of a mesh, i.e. a face belonging to a single subcell (in other words, a free face), the smoothing domains consists of a single pyramid because the nodes of a face can only be connected to one barycenter.

[0110] 2.3) Determining a Strain-Displacement Relation for Each Virtual Smoothing Domain

[0111] In this step, for each of the virtual smoothing domains as determined in the previous substep, a strain-displacement relation is determined according to the face-based smoothed finite-element method applied to the virtual smoothing domains as determined above. “Applied to the smoothing domains” means that the face-based smoothed finite-element method is applied as described in general terms in the literature, but by considering the smoothing domains as determined above, i.e. for hexahedral subcells.

[0112] According to an implementation of the invention, said strain-displacement relation takes the form of a matrix relating the displacement vector to the strain vector. This matrix is also referred to as strain-displacement transformation matrix.

[0113] According to an implementation of the invention, using Equation (2) above for the definition of the smoothed strain leads to the following expression for the strain-displacement transformation matrix, denoted by B (see also Equations 4.29, 4.30 and 4.31 in document (Liu and Nguyen Thoi Trung, 2010)):

[00003]$\begin{array}{cc}.Math.B=\left[B_1.Math.B_2.Math.\mathrm{.Math.}.Math..Math.B_{n-1}.Math.B_n\right]& \left(3\right)\\ .Math.\mathrm{with}& \\ .Math.B_i=\frac{1}{\text{?}}.Math.\int \text{?}.Math.{\mathrm{nN}}_i.Math.\mathrm{dT}& \left(4\right)\\ \text{?}.Math.\text{indicates text missing or illegible when filed}& \end{array}$

where B.sub.1 . . . B.sub.n are the matrices respectively associated with nodes i=1 . . . n of the subcell(s) to which smoothing domain Ω.sub.k.sup.s belongs, A.sub.k.sup.s is the volume of smoothing domain Ω.sub.k.sup.s, ∂Ω.sub.k.sup.s is the boundary of the smoothing domain, n is the outward normal of smoothing domain Ω.sub.k.sup.s and N.sub.i is a shape function (see paragraph 2.2.1 of document (Zienkiewicz and Taylor, 2000)) associated with node i.

[0114] 2.4) Determining a Transition Relation Between the Degrees of Freedom of the Subcell Nodes and of the Cell Nodes for Each Virtual Smoothing Domain

[0115] This substep, which is applied for each of the smoothing domains determined as described above, consists in determining a transition relation between the degrees of freedom of the nodes (located at the ends of the subcell edges), real or virtual, of the subcell(s) containing the smoothing domain considered and the degrees of freedom of the nodes of the (parent) cell(s) to which the subcell(s) considered belong.

[0116] According to an implementation of the invention, said transition relation can take the form of a transformation matrix or a transition matrix. According to this implementation of the invention, this transition matrix allows to express the degrees of freedom of the nodes (real or virtual) of the virtual cells in relation to the degrees of freedom of the real nodes of the parent cells of the mesh.

[0117] According to an implementation of the invention, the transition matrix relative to a smoothing domain is determined according to at least shape function matrices of the finite element assigned to the cell(s) to which said smoothing domain belongs, the shape function matrices being evaluated at each node of the subcell(s) to which said smoothing domain belongs.

[0118] According to an implementation of the invention, this substep is carried out for a virtual smoothing domain Ω.sub.virt as follows:

[0119] a) determining the subcell(s) i containing the smoothing domain considered. If the face on which the smoothing domain considered was generated is a free face (i.e. if the face is located on the border of the mesh), the smoothing domain is contained in a single subcell (and we note i=1 hereafter). If the face is common to two subcells, the smoothing domain is contained in two subcells (and we note i=1, 2 hereafter),

[0120] b) for each subcell i (i=1 or i=1, 2) thus identified, calculating a transition matrix T.sub.i by means of the shape functions N of the parent cell that contains the subcell considered according to a formula of the type (see Equation 2.1 of document (Zienkiewicz and Taylor, 2000)):

[00004]$\begin{array}{cc}{T}_i=\left[\begin{array}{c}{N}_i\left(x_1\right)\\ \mathrm{.Math.}\\ N_i\left(x_m\right)\end{array}\right]& \left(5\right)\end{array}$

where N.sub.i(x.sub.1) is the shape function matrix (see paragraph 2.2.1 of document (Zienkiewicz and Taylor, 2000)) of the finite element assigned to the parent cell that contains subcell i (i=1 or i=1, 2) evaluated at node j of subcell i, with j=1 . . . n, and n being the number of nodes (real and virtual) of subcell i,

[0121] c) in the case of a free face (i=1), transition matrix T.sub.Ωvirt relative to virtual smoothing domain Ω.sub.virt is equal to matrix T.sub.i as determined above; otherwise, transition matrix T.sub.Ωvirt relative to virtual smoothing domain Ω.sub.virt is constructed from the two matrices T.sub.i calculated above according to a formula of the type:

[00005]$\begin{array}{cc}{T}_{{\Omega }_{\mathrm{virt}}}=\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]& \left(6\right)\end{array}$

[0122] 2.5) Determining a Stiffness and Nodal Forces for Each Virtual Smoothing Domain

[0123] This substep consists in determining a stiffness and nodal forces for each of the smoothing domains determined as described above, from at least the transition relation determined for the smoothing domain considered as described above and the strain-displacement relation determined for this smoothing domain as described above.

[0124] According to an implementation of the invention, said stiffness and said nodal forces can take the form of a stiffness matrix and of a nodal force vector respectively.

[0125] According to an implementation of the invention, stiffness matrix K.sub.Ωvirt and nodal force vector f.sub.Ωvirt are determined for a virtual smoothing domain Ω.sub.virt as follows:

K.sub.Ω.sub.virt=∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.TDB.sub.Ω.sub.virtT.sub.Ω.sub.virtdV  (7)

f.sub.Ω.sub.virt=−∫.sub.Ω.sub.virtN.sup.TbdV−∫.sub.BΩ.sub.virtN.sup.TtdA−∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.TDε.sub.0dV+∫.sub.Ω.sub.virtT.sub.Ω.sub.virt.sup.TB.sub.Ω.sub.virt.sup.Tσ.sub.0dV  (8)

wherein transformation matrix T.sub.Ω.sub.virt is the transition matrix determined as described in substep 2.4, B.sub.Ω.sub.virt is a strain-displacement transformation matrix determined as described in substep 2.3, N is the shape function matrix of the finite element of the parent cell as described above, D is a matrix representative of the material stiffness, b is the body force vector, t the surface force vector (distributed external loading), ε.sub.0 the initial strain (see paragraph 2.2.3 of document (Zienkiewicz and Taylor, 2000)) and σ.sub.0 the initial residual stress (see paragraph 2.2.3 of document (Zienkiewicz and Taylor, 2000)).

[0126] 2.6) Determining a Stiffness and Nodal Forces for Non-Hexahedral Cells

[0127] This substep is optional and it can be applied in cases where the mesh representative of the state of the basin considered is a hex-dominant mesh, predominantly comprising hexahedral cells, but also comprising a minority of other types of cell (such as tetrahedral, pyramidal, pentahedral cells, etc.).

[0128] In this case, the face-based smoothed finite-element method as described in document (Liu and Nguyen Thoi Trung, 2010) can be applied to the other cell types of the mesh. Applying this method to the other cell types leads to define smoothing domains for each face of these cells (without cell subdivision) and to determine a stiffness matrix (denoted by K.sub.Ω hereafter) and a nodal force vector (denoted by f.sub.Ω hereafter) for each smoothing domain relative to each of the faces of the other cell types of the mesh.

[0129] According to an implementation of the invention, a smoothing domain for a face belonging to one or two cells of another type than the hexahedral type can be constructed by connecting each node of the face considered with at least one additional point located at the barycenter of the cell(s) to which the face belongs.

[0130] According to an implementation of the invention, when a face is common to a hexahedral subcell and a cell of another type, this face is considered as a free face for the construction of smoothing domains. In other words, for such a face, two smoothing domains are constructed, a first one on the side of the hexahedral subcell and a second on the side of the cell of another type, each smoothing domain being constructed by connecting each node of the face considered to the barycenter of the cell to which the node belongs.

[0131] According to an implementation of the invention, stiffness matrix K.sub.Ω and nodal force vector f.sub.Ω relative to a smoothing domain constructed for non-hexahedral cells can be determined by means of Equations 2.24a and 2.24b of document (Zienkiewicz and Taylor, 2000) and by using the definition of the strain-displacement transformation matrix denoted by B, described in Equation 3.67 of document (Liu and Nguyen Thoi Trung, 2010).

[0132] 2.7) Determining a Stiffness and Nodal Forces for the Mesh

[0133] This substep consists in determining a stiffness and nodal forces relative to the mesh representative of the state considered from at least the stiffness and the nodal forces determined for each virtual smoothing domain defined as described above for this mesh.

[0134] According to an implementation of the invention wherein the stiffness and/or the nodal forces can take the form of a stiffness matrix and a nodal force vector respectively, stiffness matrices K.sub.Ωvirt and vectors f.sub.Ωvirt determined for each virtual smoothing domain Ω.sub.virt in the previous substep are assembled by means of a standard assembly used in the finite-element method (see paragraph 1.3 of document (Zienkiewicz and Taylor, 2000)), which can be written as follows:

K=Σ.sub.i=1.sup.nvK.sub.Ω.sub.virt t  (9)

f=Σ.sub.i=1.sup.nvf.sub.Ω.sub.virt t  (10)

where nv is the total number of virtual smoothing domains Ω.sub.virt.

[0135] According to an implementation of the invention wherein the mesh representative of the state considered predominantly comprises hexahedral cells and a minority of other types of cell, a stiffness matrix and a nodal force vector are determined for the mesh according to a formula of the type:

K=Σ.sub.i=1.sup.nvK.sub.Ω.sub.virt t+Σ.sub.i=1.sup.nK.sub.Ω.sub.t  (11)

f=Σ.sub.i=1.sup.nvf.sub.Ω.sub.virt t+Σ.sub.i=1.sup.nf.sub.Ω.sub.t  (12)

where nv is the total number of virtual smoothing domains Ω.sub.virt constructed for the hexahedral cells and n is the total number of smoothing domains Ω constructed for the other cell types as described in substep 2.6.

[0136] 2.8) Determining the Displacement Field and the Stress Field for the Mesh

[0137] This substep consists in modelling the basin studied by determining at least the displacement field and the stress field for the mesh considered, by means of the numerical basin simulation according to the invention and of at least the stiffness and the nodal forces determined for the mesh considered and determined as described in the previous substep.

[0138] According to an implementation of the invention, the discretized problem defined by Equation (1) above, implemented in the basin simulator according to the invention and used with stiffness matrix K determined as described above (and computed with Equations (9) or (11) above for example) and nodal force vector f determined as described above (and computed with Equations (10) or (12) above for example), is solved. Determining vector r of Equation (2), which is the boundary conditions vector, can be done using the conventional FEM method, as described for example in paragraph 1.3 of document (Zienkiewicz and Taylor, 2000).

[0139] According to the invention, at least substeps 2.1 to 2.5 and substeps 2.6 and 2.7 are repeated for each mesh representative of each state of the basin, substep 2.6 being also repeated in the case of a mesh comprising hexahedral cells and other types of cell.

[0140] 3) Exploiting the Hydrocarbons of the Basin

[0141] This step is carried out within the context of the second aspect of the invention, which concerns a method of exploiting the hydrocarbons present in a sedimentary basin.

[0142] After carrying out the previous steps, basin simulation results are available. The basin simulation according to the invention allows at least to determine the basin stress and strain field for various states of the basin, in a stable manner in the case of a hex-dominant mesh some cells of which have a poor aspect ratio.

[0143] Implicitly, and as is conventional in basin simulation, the amount of hydrocarbons present in each cell of the mesh representative of the basin for the past and current times is also known.

[0144] Furthermore, depending on the basin simulator used for implementing the invention, it is possible for example to have information on:

[0145] i. the development of sedimentary layers

[0146] ii. their compaction under the effect of the weight of the overlying sediments

[0147] iii. their temperature evolution during burial

[0148] iv. the changes in fluid pressure resulting from this burial

[0149] iv. the thermogenic formation of hydrocarbons.

[0150] From such information, the specialist can then determine cells of the mesh representative of the basin at the current time comprising hydrocarbons, as well as the proportion, the nature and the pressure of the hydrocarbons trapped therein. The zones of the basin studied having the best petroleum potential can then be selected. These zones are then identified as reservoirs (or hydrocarbon reservoirs) of the sedimentary basin studied.

[0151] This step consists in determining at least one scheme for exploiting the hydrocarbons contained in the sedimentary basin studied. Generally, an exploitation scheme comprises a number, a geometry and a site (position and spacing) for injection and production wells to be drilled in the basin. An exploitation scheme can further comprise a type of enhanced recovery for the hydrocarbons contained in the reservoir(s) of the basin, such as enhanced recovery through injection of a solution containing one or more polymers, CO.sub.2 foam, etc. A hydrocarbon reservoir exploitation scheme must for example enable a high rate of recovery of the hydrocarbons trapped in this reservoir, over a long exploitation time, and requiring a limited number of wells. In other words, the specialist predefines evaluation criteria according to which a scheme for exploiting the hydrocarbons present in a sedimentary basin is considered sufficiently efficient to be implemented.

[0152] According to an embodiment of the invention, a plurality of exploitation schemes is defined for the hydrocarbons contained in one or more geological reservoirs of the basin studied, and at least one evaluation criterion is assessed for these exploitation schemes, by means of a reservoir simulator (such as the PumaFlow® software (IFP Energies nouvelles, France)). These evaluation criteria can comprise the amount of hydrocarbons produced for each of the various exploitation schemes, the curve representative of the production evolution over time for each well considered, the gas-oil ratio (GOR) for each well considered, etc. The scheme according to which the hydrocarbons contained in the reservoir(s) of the basin studied are really exploited can then correspond to the one meeting at least one of the evaluation criteria of the various exploitation schemes.

[0153] Then, once an exploitation scheme determined, the hydrocarbons trapped in the petroleum reservoir(s) of the sedimentary basin studied are exploited according to this exploitation scheme, notably at least by drilling the injection and production wells of the exploitation scheme thus determined, and by installing the production infrastructures necessary to the development of this or these reservoirs. Moreover, in cases where the exploitation scheme has been determined by estimating the production of a reservoir associated with different enhanced recovery types, the selected additive type(s) (polymers, surfactants, CO.sub.2 foam) are injected into the injection well.

[0154] It is understood that a scheme for exploiting hydrocarbons in a basin can evolve during the exploitation of the hydrocarbons of this basin, for example according to additional basin-related knowledge acquired during this exploitation and to improvements in the various technical fields involved in the exploitation of a hydrocarbon reservoir (advancements in the field of drilling, of enhanced oil recovery for example).

[0155] Equipment and Computer Program Product

[0156] The method according to the invention is implemented by means of an equipment (a computer workstation for example) comprising data processing means (a processor) and data storage means (a memory, in particular a hard drive), as well as an input/output interface for data input and method results output.

[0157] The data processing means are configured for carrying out in particular step 2 described above.

[0158] Furthermore, the invention concerns a computer program product downloadable from a communication network and/or recorded on a computer-readable medium and/or processor executable, comprising program code instructions for implementing the method as described above, and notably step 2), when said program is executed on a computer.

[0159] Thus, the method according to the first aspect of the invention allows to stabilize a numerical basin simulation solving at least one balance equation of poromechanics according to a face-based smoothed finite-element method, in the case of a hex-dominant mesh, even in the case of fine and/or pinched-out geological layers in the sedimentary basin.

[0160] Besides, the method according to the second aspect of the invention allows to predict the petroleum potential of complex sedimentary basins, which may have undergone for example complex tectonic movements and exhibit fine layers and/or stratigraphic pinchouts, which contributes to improving the exploitation of hydrocarbons in this type of sedimentary basins.