# Finite Element Simulation of Oil and Gas Reservoir In Situ Stress Based on a 3D Corner-point Grid Model

Authors: Liu Yuyang; Liu Shiqi; Pan Mao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A three-dimensional (3D) corner-point grid model gives a relatively accurate description of the structural properties and spatial distribution of oil and gas reservoirs than Cartesian grids. The finite element simulation of the stress field provides a relatively probable presentation of the in situ stress distribution. Both methods are of great importance to the exploration and development of oil and gas fields. Implementing the finite element simulation of in situ stress on a 3D corner-point grid model not only retains the structural attributes of a reservoir but also allows the accurate simulation of the 3D stress distribution. In this paper, we present a method for implementing the finite element simulation of in situ stress based on a 3D corner-point grid model. We first established a fine 3D reservoir model with corner-point grids and then converted the grids into corresponding 3D finite element grid models using a grid conversion algorithm. Next, we simulated the in situ stress distribution with the finite element method. The stress model is then resampled to corresponding corner-point grid geological models using the reverse algorithm. The grid conversion algorithm is to provide data support for the subsequent numerical simulation and other research efforts, thereby guaranteeing procedure continuity and data consistency. Finally, we simulated the stress distribution of a real oil field, theX region. Comparing the simulated result with the measured result, the high agreement validated the effectiveness and accuracy of the proposed method.


---

## Body

## 1. Introduction

The in situ stress, also known as the crustal stress, is a force existing in rock in its natural state that includes the gravitational force, tectonic stress, and residual stress [1]. Regional tectonic movements, fluid overpressures, and a range of other physical and chemical processes over the long geological history can cause continuous readjustment and evolution in the in situ stress field, giving rise to spatially and temporally complex distribution patterns of the in situ stress [1, 2]. The development of unconventional reservoirs, mainly refers to the shale reservoirs, tight gas and oil, and coalbed methane, often involves mass hydraulic fracturing operations, and the magnitude and direction of the in situ stress control the formation and distribution of artificial fractures in oil and gas field development [2–4]. The stress field distribution also plays a critical role in the well-bore stability of drilling operations [5]. Needless to say, investigating the in situ stress distribution is of vital importance to oil and gas field exploration and development [6, 7].The corner-point grid is a non-Cartesian but orthogonal areal grid using corner-point geometry, specifying the corners of each grid block in grid building to represent complicated reservoir geometries [8–10]. It is a technique introduced by Ding and Lemonnier [8] and applied to the numerical simulation of reservoirs by Ponting [11], which is widely used in geological modeling and numerical simulation of reservoirs [6, 12]. Over the past few years, studies on the corner-point grid have become more refined with the development of technology [13–15]. The corner-point grid is superior to its conventional counterparts in describing the heterogeneity of a reservoir since it avoids the inflexibility of the conventional orthogonal grids while providing a precise description of the interfaces, fault planes, and pinch-outs [9, 10, 16], an accurate presentation of the spatial distribution, and the attributes of the reservoir [17, 18]. Unfortunately, however, because of the data organization and the grid block attribute treatment [8, 11], little has been reported on three-dimensional (3D) in situ stress simulation with corner-point grids. Rather, most studies so far have relied on the finite element method (FEM) to simulate 3D stress distributions with finite element grids because of extensive applicability and practical significance of FEM in solving continues media and field problems [19].FEM is a numerical simulation procedure that is able to analyze and solve continuous fields [20–22]. FEM is among the most widely applied approaches for the numerical simulation of in situ stresses [23, 24]. It is possible to establish a relatively accurate 3D in situ stress model based on finite element grids using FEM simulation by establishing a finite element model with boundary loads and boundary displacements constraints [25]. Many authors worldwide have succeeded in simulating stress fields in geologically complicated areas including coal mines [26, 27] and tunnels [28, 29], but most of these studies have directly relied on finite element software to establish their models. For oil and gas reservoirs, which have complicated structures, great burial depths, and high heterogeneous attributes, modeling is always related to some geological statics algorithms; it may not establish a sufficiently accurate 3D model with FEM software for the reservoir scale [25, 30, 31]. Some studies also report these problems [32, 33], and approaches have been implemented to resolve these obstacles [34].Nevertheless, as the grid generation approaches vary among different types of grids, a geological model based on corner-point grids cannot be directly used for finite element simulation. Because of limitations in the modeling features of finite element simulation programs, it may not use these programs to establish a model that provides a fine description of the reservoir attributes for the reservoir scale; geological modeling programs often use corner-point grids to provide a fine description and characterization of reservoirs. Many studies have been carried out to resolve this obstacle [35]. In this paper, we present a method for implementing a finite element simulation of in situ stress based on a 3D corner-point grid model combining the merits of both the corner-point and FEM grids. Using this method, we were able to build fine 3D geological models based on corner-point grids and simulate the 3D stress distribution with FEM. This both effectively retains the reservoir features and allows a relatively accurate simulation of the 3D stress distribution of the study area.After establishing fine 3D geological models based on a corner-point grid, including tectonic framework models, facies models, and attribute models, we designed and implemented a grid conversion algorithm and, while retaining the structural and attribute distributions of the reservoir, directly converted the 3D models based on a corner-point grid into 3D finite element models that can be applied to finite element simulation. After that, FEM simulators are employed to simulate the in situ stress distribution based on the 3D finite element models. Finally, through a postprocessing algorithm for grid conversion, we resampled the in situ stress back into the fine 3D geological models based on corner-point grids. Thus, the finite element simulation of in situ stress based on a 3D corner-point grid model is fulfilled.In this paper, we first designed and outlined an overall procedure for implementing the finite element simulation of in situ stress based on the 3D corner-point grid model. Then, we performed a detailed analysis of the data structures of corner-point grids and finite element grids and included a preliminary analysis of the grid conversion. Next, we elaborated on the grid conversion algorithm and the simulation procedure; finally, we simulated the stress distribution in oil fields in regionX and compared our simulated result with the measured result. The result validates the effectiveness and accuracy of the proposed method.


## 2. Overall Procedural Design

The finite element numerical simulation of an in situ stress based on a 3D corner-point grid geological model typically involves the following procedure: first, structural models and fine 3D geological attribute models are established based on corner-point grids, which form the basis for simulating the in situ stress. Specifically, this step consists of building 3D geological grid models through 3D geological structural modeling with logging, seismic, and experimental data and building rock mechanical parameter field models through attribute modeling interpolation algorithms.Then, the resulted corner-point grid models are converted into corresponding finite element grid models using a grid conversion algorithm to allow the simulation of the stress distribution with the FEM. The conversion algorithm extracts a series of information of all nodes of the corner-point model, including node number, node coordinate, and node attribute. The mapping relationships of nodes and grid blocks are also extracted and saved. The algorithm then sorts all nodes, and a serious modification (duplicated nodes operation, mapping relationship construction, and so on) is operated to set up a new node model based on finite element grid. Meanwhile, considering about the aforementioned mapping relationship, a new mapping relationship between the finite element grid nodes and blocks is established. Furthermore, a mapping relationship between the corner-point grid block attributes and finite element grid block material properties has been established. The attribute of each grid block of corner-point grid is also connected to the material properties of the corresponding cells of finite element grid through the mapping relationship. Thus, the corner-point grid is converted into a finite element grid. The finite element models contain the finite element node model (stored the node information), the finite element cell model (stored the cells information and the mapping relationship between nodes and cell), and the finite element cell material properties models (stored each cells’ material properties related the rock mechanic attributes, including Young’s Modulus, Poisson ratio, and rock density).Finally, the in situ stress is simulated by using an FEM simulator and the simulated finite element in situ stress models are converted back into corresponding corner-point grid models using a reverse algorithm for grid conversion to facilitate subsequent analysis and further simulation efforts such as sweet spot area evaluation or well pattern optimization. This conversion algorithm does not introduce any instability and reduces accuracy result in the model. The detailed reverse algorithm is almost similar to that of the conversion from corner-point grid to finite grid while it is an opposite process. The in situ stress attributes including magnitude and direction are also linked to each grid blocks. The overall procedure for the design is shown in Figure1.Figure 1
Overall procedure design.As both corner-point grids and finite element grids will be used during the simulation, the continuity of the simulation procedure is heavily challenged by the difficulty in converting and discriminating between these two types of grids due to their highly different forms of definition and data structures. The whole procedure is divided into two different parts: the preprocessing procedure and postprocessing procedure.By designing and implementing a grid conversion algorithm, the main part of the preprocessing procedure, the 3D corner-point grid models are converted into finite element models. Also, the rock mechanical parameters of the attribute models of each cell correspond to the material parameters of the finite element grid cells. Dealing with complicated geological modeling with multilevel faults, the FEM simulator often cannot establish relatively accurate model than geological modeling software due to its complicated morphology and topology, and corner-point grid-based modeling is widely employed among most of the geological software to build up and represent complicated geological models [32]. Meanwhile, in order to establish a relatively accurate heterogeneous attribute model, some geological statistics-based modeling approaches and algorithms are needed [30, 32]. The current existing FEM simulator cannot directly set up a fine heterogeneous attribute model-based geological statistics methods to represent heterogeneous attribute distribution [30, 32]. This approach not only allows the accurate preservation of the features of the fine geological models and attribute models established by the geological modeling program but also mitigates the low modeling and grid subdivision capabilities that are challenges in the use of finite element programs [8–10, 19].By defining the related constraints with respect to displacements or boundary loads, it is possible to use FEM to simulate the 3D in situ stress distribution [32, 33] and, through multiple inversions of the constraints at known target points through specific mathematical models [1] or artificial neural network, obtain the more probable distribution of the 3D in situ stress. The finite element in situ stress model is usually used for the numerical simulation of reservoirs, pattern arrangement of wells, and hydraulic fracture simulation. A reverse algorithm for grid conversion is then designed and implemented to resample the stress attributes of the 3D finite element stress field back into the fine 3D corner-point grid geological models. The corner-point grid-based model contains both the magnitude and direction of the in situ stress and can be used in subsequent analysis to guarantee procedure continuity and data consistency.


## 3. Grid Data Structure

### 3.1. Corner-Point Grid

The corner-point grid and finite element grid are grid models based on the grid generation algorithm. Grid generation is mainly used for numerical simulation and computation. Different means of grid generation algorithm generate different types of grids [36]. Grids can be categorized into orthogonal and nonorthogonal according to whether or not they are orthogonal, the Cartesian grid and tetrahedron grid, into structured, nonstructured, and hybrid structured-and-nonstructured grids according to whether or not their nodes are well arranged, the corner-point grid, tetrahedron grid, and perpendicular bisection grid, and into adaptive and nonadaptive grids according to whether or not they adapt to time intervals during simulation [37, 38].A corner-point grid is categorized as an orthogonal, structured, and nonadaptive grid [10, 36]. In a structured grid, all nodes in the grid region share the same number of neighboring grid blocks, each horizontal or longitudinal row shares the same number of nodes, each node at the model’s boundary has two cells around it, and each internal node has four cells around it [11, 14]. The structural properties of the corner-point grid make it easy to fit the regional boundaries of a reservoir and to simulate the fluid in a finite difference simulator for describing the microtectonic morphology of oil and gas reservoirs, the boundaries of reservoirs, the type of flow, and horizontal wells, directional wells, and faults. Structured grids are also applauded for their simple grid generation algorithms, high generation speed, good generation quality, and simple and well-organized data structures [39].By definition, a corner-point grid is an irregular hexahedral grid [11]. The morphology of its cell blocks are described by four coordinate lines (defined by the top and bottom regular topologic control planes) and the coordinates of eight grid nodes (corner points), as shown in Figures 2(a) and 2(b). Two Nx+1×Ny+1 regular topologic control planes control the generation of individual structural cells; the control plane and intermediate points define the top and bottom boundary coordinates of each cell with Nz−1 gliding lines, which logically constitute a regular topological model with Nx×Ny×Nz cells, as shown in Figures 2(c) and 2(d). For the depth coordinates of a corner point, a grid describing the distribution of Cartesian coordinates in the computational domain is directly generated by the interpolation algorithm in formula (1) according to the coordinate system of the upper and lower boundaries of the computational area. The coordinates are XA,YA,ZA for point A, XE,YE,ZE for point E, XI,YI,ZI for point I, and XJ,YJ,ZJ for point J. Logically, there are Nx+1, Ny+1, and Nz+1 coordinate lines in the X, Y, and Z (depth) directions, respectively, and the grid model is divided into Nx, Ny, and Nz grid cells:(1)Xk=XI+Zk−ZI×XJ−XIZJ−ZI,k=A,E,Yk=YI+Zk−ZI×YJ−YIZJ−ZI,k=A,E.Figure 2
Schematic of a corner-point grid: coordinates (a), cells (b), transverse index (c), and vertical index (d).
(a)
(b)
(c)
(d)


### 3.2. Finite Element Grid

FEM is a scheme that discretizes the equation of a target region such that the problem is converted into a number of finite elements that are interconnected by nodes and solved in these individual finite elements [22, 32, 40]. A finite element grid is a structural unit formed through a grid subdivision algorithm during the application of FEM. It can be a structured grid or a nonstructured one. The former grids are typically hexahedral, and the latter are typically tetrahedral; both are widely applied [41]. Figure 3(a) shows the conversion between a structured grid (hexahedral) and a nonstructured grid (tetrahedral). Here, two grid cells, I and II, are shown. The former contains finite element nodes 1-2-4-5-7-8-10-11, and the latter contains finite element nodes 2-3-5-6-8-9-11-12. Grid cell I can be converted into four corresponding tetrahedral structural units: i (nodes 1-10-8-7), ii (nodes 1-8-5-2), iii (nodes 1-5-10-4), and iv (nodes 1-8-10-5).Figure 3
Schematic of a finite element grid: (a) conversion from a structured grid to a nonstructured finite element grid and (b) data structure of a corner-point grid versus that of a finite element grid.
(a)
(b)A finite element grid is connected by nodes. It is distinguished from a corner-point grid in the way it deals with the nodes between adjacent cells. For finite element grids, to ensure the continuous transmission of forces and other variables, the same grid nodes are shared between adjacent grid blocks; for corner-point grids, the nodes between adjacent grids are relatively independent and each grid has its own grid nodes. Figure3(b) shows a simplified grid system that contains only grid cells I and II. For finite element grids, the entire grid system contains two grid cells (I and II) and 12 nodes (1–12), where nodes 2-5-11-8 are shared between cells I and II. For corner-point grids, however, the grid system still contains two grid cells but 16 nodes—cell I contains nodes 1-2-4-5-7-8-10-11 and cell II contains nodes 2′-3-6-5′-8′-9-12-11’. Among these nodes, 2-2′, 5-5′, 8-8′, and 11-11′ are nodes that are included in each of the adjacent grids. Despite the identical position information, they belong to different grid cells.


## 3.1. Corner-Point Grid

The corner-point grid and finite element grid are grid models based on the grid generation algorithm. Grid generation is mainly used for numerical simulation and computation. Different means of grid generation algorithm generate different types of grids [36]. Grids can be categorized into orthogonal and nonorthogonal according to whether or not they are orthogonal, the Cartesian grid and tetrahedron grid, into structured, nonstructured, and hybrid structured-and-nonstructured grids according to whether or not their nodes are well arranged, the corner-point grid, tetrahedron grid, and perpendicular bisection grid, and into adaptive and nonadaptive grids according to whether or not they adapt to time intervals during simulation [37, 38].A corner-point grid is categorized as an orthogonal, structured, and nonadaptive grid [10, 36]. In a structured grid, all nodes in the grid region share the same number of neighboring grid blocks, each horizontal or longitudinal row shares the same number of nodes, each node at the model’s boundary has two cells around it, and each internal node has four cells around it [11, 14]. The structural properties of the corner-point grid make it easy to fit the regional boundaries of a reservoir and to simulate the fluid in a finite difference simulator for describing the microtectonic morphology of oil and gas reservoirs, the boundaries of reservoirs, the type of flow, and horizontal wells, directional wells, and faults. Structured grids are also applauded for their simple grid generation algorithms, high generation speed, good generation quality, and simple and well-organized data structures [39].By definition, a corner-point grid is an irregular hexahedral grid [11]. The morphology of its cell blocks are described by four coordinate lines (defined by the top and bottom regular topologic control planes) and the coordinates of eight grid nodes (corner points), as shown in Figures 2(a) and 2(b). Two Nx+1×Ny+1 regular topologic control planes control the generation of individual structural cells; the control plane and intermediate points define the top and bottom boundary coordinates of each cell with Nz−1 gliding lines, which logically constitute a regular topological model with Nx×Ny×Nz cells, as shown in Figures 2(c) and 2(d). For the depth coordinates of a corner point, a grid describing the distribution of Cartesian coordinates in the computational domain is directly generated by the interpolation algorithm in formula (1) according to the coordinate system of the upper and lower boundaries of the computational area. The coordinates are XA,YA,ZA for point A, XE,YE,ZE for point E, XI,YI,ZI for point I, and XJ,YJ,ZJ for point J. Logically, there are Nx+1, Ny+1, and Nz+1 coordinate lines in the X, Y, and Z (depth) directions, respectively, and the grid model is divided into Nx, Ny, and Nz grid cells:(1)Xk=XI+Zk−ZI×XJ−XIZJ−ZI,k=A,E,Yk=YI+Zk−ZI×YJ−YIZJ−ZI,k=A,E.Figure 2
Schematic of a corner-point grid: coordinates (a), cells (b), transverse index (c), and vertical index (d).
(a)
(b)
(c)
(d)


## 3.2. Finite Element Grid

FEM is a scheme that discretizes the equation of a target region such that the problem is converted into a number of finite elements that are interconnected by nodes and solved in these individual finite elements [22, 32, 40]. A finite element grid is a structural unit formed through a grid subdivision algorithm during the application of FEM. It can be a structured grid or a nonstructured one. The former grids are typically hexahedral, and the latter are typically tetrahedral; both are widely applied [41]. Figure 3(a) shows the conversion between a structured grid (hexahedral) and a nonstructured grid (tetrahedral). Here, two grid cells, I and II, are shown. The former contains finite element nodes 1-2-4-5-7-8-10-11, and the latter contains finite element nodes 2-3-5-6-8-9-11-12. Grid cell I can be converted into four corresponding tetrahedral structural units: i (nodes 1-10-8-7), ii (nodes 1-8-5-2), iii (nodes 1-5-10-4), and iv (nodes 1-8-10-5).Figure 3
Schematic of a finite element grid: (a) conversion from a structured grid to a nonstructured finite element grid and (b) data structure of a corner-point grid versus that of a finite element grid.
(a)
(b)A finite element grid is connected by nodes. It is distinguished from a corner-point grid in the way it deals with the nodes between adjacent cells. For finite element grids, to ensure the continuous transmission of forces and other variables, the same grid nodes are shared between adjacent grid blocks; for corner-point grids, the nodes between adjacent grids are relatively independent and each grid has its own grid nodes. Figure3(b) shows a simplified grid system that contains only grid cells I and II. For finite element grids, the entire grid system contains two grid cells (I and II) and 12 nodes (1–12), where nodes 2-5-11-8 are shared between cells I and II. For corner-point grids, however, the grid system still contains two grid cells but 16 nodes—cell I contains nodes 1-2-4-5-7-8-10-11 and cell II contains nodes 2′-3-6-5′-8′-9-12-11’. Among these nodes, 2-2′, 5-5′, 8-8′, and 11-11′ are nodes that are included in each of the adjacent grids. Despite the identical position information, they belong to different grid cells.


## 4. Grid Conversion Algorithm and Implementation of Simulation Procedure

### 4.1. Grid Conversion Algorithm

We designed a grid conversion algorithm that makes conversion between the two grids possible and guarantees the procedure continuity and data consistency. First, we defined two structural bodies, GPoint{} and GVolume{}, to store the node information and grid information, respectively. GPoint{}, which is used to store the position information of each point in the corner-point grid, mainly contains theX, Y, and Z coordinate information; GVolume{}, which is used to store cell-related information, mainly contains the numbering of the grids, the corner points contained in the grids, and the attribute information of the grids (including the grid effectiveness and density and attributes of the rock mechanical parameters).The following functions are contained in the bidirectional grid conversion algorithm (Figure4): PreProcessAlgorithm() is used to optimize the input corner-point grid model files and check their validity. InputModelOperation() is used to process the input corner-point grid models and record information relating to them. GenerateModelNode() is used to convert corner-point grid models into finite element grid nodes. GenerateFiniteGrid() is used to establish 3D finite element grid models with the resulting finite element nodes. GenerateMaterialNumber() and GenerateMateiralAttribute() are used to number the 3D finite element grids and define the material attribute parameters of each finite element grid against attribute information such as the grid density and rock mechanical parameters at the corresponding corner points. GenerateLoad() and GenerateDisplacement() are used to define the grid boundary loads and boundary constraints. StressOperation() involves multiple reversions of the model boundary loads and boundary constraints at known target points and the execution of the finite element simulation algorithm to obtain a 3D finite element in situ stress model. OutputCtrl() is used to automatically store the coordinates as well as the stresses and strains of the nodes in each finite element grid and to process and store the magnitude and direction of the stress in each finite element grid. Finally, PostProcessAlgorithm() reads the stored 3D finite element in situ stress models and converts them into 3D in situ stress models based on the corner-point grid.Figure 4
Algorithm functions and their uses.


### 4.2. Preprocessing Procedure

First, PreProcessAlgorithm() is used to check the validity of the input model. Then, InputModelOperation() is used to read the original 3D corner-point grid models, record all node and grid information, and extract the node information from the model. Next, GenerateModelNode() is used to deal with the repetitive grid nodes and convert them into corresponding finite element model nodes. After obtaining a finite element node model, GenerateFiniteGrid() is used to read the grid validity information in the grid model, correlate the nodes corresponding to these valid grids, and obtain a 3D finite element grid model. Finally, GenerateMateiralNumber() and GenerateMateiralAttribute() are used to correlate the attribute information in the grids to the material attributes of the finite element grids, number the grid materials in the individual grids, assign an attribute to each grid material (including Young’s modulus, Poisson ratio, and density), and obtain a complete 3D finite element model, as shown in Figure5.Figure 5
Schematic flowchart showing the preprocessing procedure.We tested the 3D corner-point grids without inactive grids (grids are useful in data organization while are not useful in model description and analysis) and with inactive grids using the grid conversion algorithm, read each of them into the corner-point grid model, and converted the grid nodes into the corresponding finite element node models. According to the correspondences between the nodes and grids and the activity of the grids, we connected the nodes into a finite element grid, thereby successfully converting a corner-point grid (Figure6 shows in the corner-point grid, the cyan grid blocks are active while the blue ones are inactive; left column is the corner-point grid without inactive grid blocks and right column is the corner-point grid with inactive grid blocks) into a finite element grid (Figure 6 shows the finite element grid in which the black dots are the model nodes and cyan blocks are the grid blocks). For each grid, we numbered the material attributes corresponding to the grid and correlated the parameters of the grid attributes so that different grid attributes (Young's modulus, Poisson ratio, and density) are assigned (Figure 6 shows number material attributes, in which each color refers to a unique material attribute, and assign material attribute, in which each grid block has a unique material attribute), and the heterogeneous within a reservoir or rock body is effectively retained, as shown in Figure 6.Figure 6
Case study of the preprocessing procedure: the corner-point grid without inactive grid blocks (left column) and with inactive grid blocks (right column, the cyan grid blocks are active and blue ones are inactive); the finite element grid nodes (black dots) and grid blocks (cyan cells); every color refers to a individual material attribute and every grid block has a unique color (some colors are duplicated due to limited number of colors); every grid block has a unique material attribute.Considering about a simplified grid model, there are two independent grid blocks (I and II) and 16 nodes (1-12, 2′, 5′, 8′, and 11′), as shown in Figure7. If the model does not have inactive grid, as shown in Figure 7(a), for the corner-point grid, node 1-8 belongs to grid block I and nodes 2′, 3, 5′, 6, 8′, 9, 11′, and 12 belong to grid block II. However, for the corresponding finite element grid, grid block I contains nodes 1, 2, 4, 5, 7, 8, 10, and 11 and grid block II contains nodes 2, 5, 8, 11, 3, 6, 9, and 12. Grid block I and II share the same node of 2, 5, 8, and 11. Talking about the conversion algorithm, nodes are first sorted along Nx, Ny, and Nz direction, and some nodes are eliminated through the original mapping relationship, and then a new node sequence and mapping relationship between node and grid block is settled. Thus, a new model has formed based on the finite element grid. Considering about the model with inactive grid blocks, for example, grid block I is inactive, as shown in Figure 7(b), a flag is marked for the grid block I during the process of nodes and grid blocks. When the new mapping relationship is rebuilt, the grid block with the flag is not connected with its related nodes. As a result, only the active grid blocks organize the new model. Some case models are also employed as shown in Figure 8.Figure 7
Illustration of grid conversion from corner-point grid to the finite element grid using a conceptual model without and with inactive grid blocks: (1) model without inactive grid block and (2) model with inactive grid blocks.
(a)
(b)Figure 8
Illustration of grid conversion from the corner-point grid to the finite element grid using two case study models without and with inactive grid blocks. (a) Grid conversion without inactive grid blocks. (b) Grid conversion with inactive grid blocks.
(a)
(b)There are two different case models in Figure8, the model without inactive grid blocks (left column in Figure 8) and with inactive grid blocks (right column in Figure 8). For the model without inactive grid blocks, there are 84 grid blocks (7 × 6 × 2) and 168 nodes (8 × 7 × 3). For the corner-point grid-based model, there are 672 nodes stored to form the 84 grid blocks, while for the finite element grid-based model, there are only 168 nodes. For the model with inactive grid blocks, there are 126 grid blocks (9 × 7 × 2, 80 grid blocks are active and 36 of them are inactive) and 240 nodes (10 × 8 × 3). For the corner-point grid-based model, there are 1008 nodes stored to form the 126 grid blocks (80 grid blocks are active and 36 of them are inactive), while for the finite element grid-based model, there are only 240 nodes (the inactive grid blocks’ nodes are preserved to maintain the new mapping relationship with the former one) and only the 80 active grid blocks are formed.


### 4.3. FEM Simulation and Postprocessing Procedure

In the FEM simulation procedure, first, with stress or displacement constraints at known target points, the proper boundary loads and boundary constraints are automatically defined in the functions GenerateLoad() and GenerateDisplacement() using the built-in parameters in the algorithm (which is integrated in the functions GenerateLoad() and GenerateDisplacement()). Then, StressOperation() is used to establish a 3D finite element in situ stress model using the FEM simulator, ANSYS or ABAQUS, that gives a presentation of the in situ stress distribution. Next, OutputCtrl() automatically stores the simulation result and various information contained in the 3D finite element in situ stress model. Finally, in the postprocessing procedure, PostProcessAlgorithm() is used to further convert the model into a 3D in situ stress model based on the corner-point grid, including the magnitude and direction, as shown in Figure9.Figure 9
FEM simulation and postprocessing procedure.We established a finite element grid model using the preprocessing procedure. Then, according to the stress levels at known points, we inverted the grid boundary constraints and stresses, defined the loads and constraints for the target reservoir, simulated the in situ stress of the target reservoir with FEM, and resampled the resulted 3D finite element grid in situ stress models (including the magnitude and direction of the stress) into the corresponding corner-point grid models in a seamless manner to provide support for the subsequent numerical simulation of the reservoir or other efforts, thereby guaranteeing the procedure continuity, as shown in Figure10, with a case study. The variations and differences of the stress models in Figure 10 are originated from different visualization software, the finite element grid model for ABAQUS and the corner-point grid for Petrel. The values of every grid blocks are the same in both finite element grid-based stress model and corner-point grid-based stress model.Figure 10
Case study of the postprocessing procedure (finite element models are visualized by ABAQUS and corner-point grid model is visualized by Petrel).


## 4.1. Grid Conversion Algorithm

We designed a grid conversion algorithm that makes conversion between the two grids possible and guarantees the procedure continuity and data consistency. First, we defined two structural bodies, GPoint{} and GVolume{}, to store the node information and grid information, respectively. GPoint{}, which is used to store the position information of each point in the corner-point grid, mainly contains theX, Y, and Z coordinate information; GVolume{}, which is used to store cell-related information, mainly contains the numbering of the grids, the corner points contained in the grids, and the attribute information of the grids (including the grid effectiveness and density and attributes of the rock mechanical parameters).The following functions are contained in the bidirectional grid conversion algorithm (Figure4): PreProcessAlgorithm() is used to optimize the input corner-point grid model files and check their validity. InputModelOperation() is used to process the input corner-point grid models and record information relating to them. GenerateModelNode() is used to convert corner-point grid models into finite element grid nodes. GenerateFiniteGrid() is used to establish 3D finite element grid models with the resulting finite element nodes. GenerateMaterialNumber() and GenerateMateiralAttribute() are used to number the 3D finite element grids and define the material attribute parameters of each finite element grid against attribute information such as the grid density and rock mechanical parameters at the corresponding corner points. GenerateLoad() and GenerateDisplacement() are used to define the grid boundary loads and boundary constraints. StressOperation() involves multiple reversions of the model boundary loads and boundary constraints at known target points and the execution of the finite element simulation algorithm to obtain a 3D finite element in situ stress model. OutputCtrl() is used to automatically store the coordinates as well as the stresses and strains of the nodes in each finite element grid and to process and store the magnitude and direction of the stress in each finite element grid. Finally, PostProcessAlgorithm() reads the stored 3D finite element in situ stress models and converts them into 3D in situ stress models based on the corner-point grid.Figure 4
Algorithm functions and their uses.


## 4.2. Preprocessing Procedure

First, PreProcessAlgorithm() is used to check the validity of the input model. Then, InputModelOperation() is used to read the original 3D corner-point grid models, record all node and grid information, and extract the node information from the model. Next, GenerateModelNode() is used to deal with the repetitive grid nodes and convert them into corresponding finite element model nodes. After obtaining a finite element node model, GenerateFiniteGrid() is used to read the grid validity information in the grid model, correlate the nodes corresponding to these valid grids, and obtain a 3D finite element grid model. Finally, GenerateMateiralNumber() and GenerateMateiralAttribute() are used to correlate the attribute information in the grids to the material attributes of the finite element grids, number the grid materials in the individual grids, assign an attribute to each grid material (including Young’s modulus, Poisson ratio, and density), and obtain a complete 3D finite element model, as shown in Figure5.Figure 5
Schematic flowchart showing the preprocessing procedure.We tested the 3D corner-point grids without inactive grids (grids are useful in data organization while are not useful in model description and analysis) and with inactive grids using the grid conversion algorithm, read each of them into the corner-point grid model, and converted the grid nodes into the corresponding finite element node models. According to the correspondences between the nodes and grids and the activity of the grids, we connected the nodes into a finite element grid, thereby successfully converting a corner-point grid (Figure6 shows in the corner-point grid, the cyan grid blocks are active while the blue ones are inactive; left column is the corner-point grid without inactive grid blocks and right column is the corner-point grid with inactive grid blocks) into a finite element grid (Figure 6 shows the finite element grid in which the black dots are the model nodes and cyan blocks are the grid blocks). For each grid, we numbered the material attributes corresponding to the grid and correlated the parameters of the grid attributes so that different grid attributes (Young's modulus, Poisson ratio, and density) are assigned (Figure 6 shows number material attributes, in which each color refers to a unique material attribute, and assign material attribute, in which each grid block has a unique material attribute), and the heterogeneous within a reservoir or rock body is effectively retained, as shown in Figure 6.Figure 6
Case study of the preprocessing procedure: the corner-point grid without inactive grid blocks (left column) and with inactive grid blocks (right column, the cyan grid blocks are active and blue ones are inactive); the finite element grid nodes (black dots) and grid blocks (cyan cells); every color refers to a individual material attribute and every grid block has a unique color (some colors are duplicated due to limited number of colors); every grid block has a unique material attribute.Considering about a simplified grid model, there are two independent grid blocks (I and II) and 16 nodes (1-12, 2′, 5′, 8′, and 11′), as shown in Figure7. If the model does not have inactive grid, as shown in Figure 7(a), for the corner-point grid, node 1-8 belongs to grid block I and nodes 2′, 3, 5′, 6, 8′, 9, 11′, and 12 belong to grid block II. However, for the corresponding finite element grid, grid block I contains nodes 1, 2, 4, 5, 7, 8, 10, and 11 and grid block II contains nodes 2, 5, 8, 11, 3, 6, 9, and 12. Grid block I and II share the same node of 2, 5, 8, and 11. Talking about the conversion algorithm, nodes are first sorted along Nx, Ny, and Nz direction, and some nodes are eliminated through the original mapping relationship, and then a new node sequence and mapping relationship between node and grid block is settled. Thus, a new model has formed based on the finite element grid. Considering about the model with inactive grid blocks, for example, grid block I is inactive, as shown in Figure 7(b), a flag is marked for the grid block I during the process of nodes and grid blocks. When the new mapping relationship is rebuilt, the grid block with the flag is not connected with its related nodes. As a result, only the active grid blocks organize the new model. Some case models are also employed as shown in Figure 8.Figure 7
Illustration of grid conversion from corner-point grid to the finite element grid using a conceptual model without and with inactive grid blocks: (1) model without inactive grid block and (2) model with inactive grid blocks.
(a)
(b)Figure 8
Illustration of grid conversion from the corner-point grid to the finite element grid using two case study models without and with inactive grid blocks. (a) Grid conversion without inactive grid blocks. (b) Grid conversion with inactive grid blocks.
(a)
(b)There are two different case models in Figure8, the model without inactive grid blocks (left column in Figure 8) and with inactive grid blocks (right column in Figure 8). For the model without inactive grid blocks, there are 84 grid blocks (7 × 6 × 2) and 168 nodes (8 × 7 × 3). For the corner-point grid-based model, there are 672 nodes stored to form the 84 grid blocks, while for the finite element grid-based model, there are only 168 nodes. For the model with inactive grid blocks, there are 126 grid blocks (9 × 7 × 2, 80 grid blocks are active and 36 of them are inactive) and 240 nodes (10 × 8 × 3). For the corner-point grid-based model, there are 1008 nodes stored to form the 126 grid blocks (80 grid blocks are active and 36 of them are inactive), while for the finite element grid-based model, there are only 240 nodes (the inactive grid blocks’ nodes are preserved to maintain the new mapping relationship with the former one) and only the 80 active grid blocks are formed.


## 4.3. FEM Simulation and Postprocessing Procedure

In the FEM simulation procedure, first, with stress or displacement constraints at known target points, the proper boundary loads and boundary constraints are automatically defined in the functions GenerateLoad() and GenerateDisplacement() using the built-in parameters in the algorithm (which is integrated in the functions GenerateLoad() and GenerateDisplacement()). Then, StressOperation() is used to establish a 3D finite element in situ stress model using the FEM simulator, ANSYS or ABAQUS, that gives a presentation of the in situ stress distribution. Next, OutputCtrl() automatically stores the simulation result and various information contained in the 3D finite element in situ stress model. Finally, in the postprocessing procedure, PostProcessAlgorithm() is used to further convert the model into a 3D in situ stress model based on the corner-point grid, including the magnitude and direction, as shown in Figure9.Figure 9
FEM simulation and postprocessing procedure.We established a finite element grid model using the preprocessing procedure. Then, according to the stress levels at known points, we inverted the grid boundary constraints and stresses, defined the loads and constraints for the target reservoir, simulated the in situ stress of the target reservoir with FEM, and resampled the resulted 3D finite element grid in situ stress models (including the magnitude and direction of the stress) into the corresponding corner-point grid models in a seamless manner to provide support for the subsequent numerical simulation of the reservoir or other efforts, thereby guaranteeing the procedure continuity, as shown in Figure10, with a case study. The variations and differences of the stress models in Figure 10 are originated from different visualization software, the finite element grid model for ABAQUS and the corner-point grid for Petrel. The values of every grid blocks are the same in both finite element grid-based stress model and corner-point grid-based stress model.Figure 10
Case study of the postprocessing procedure (finite element models are visualized by ABAQUS and corner-point grid model is visualized by Petrel).


## 5. A Case Study of RegionX

In this section, we tested our algorithm on real oil field data of regionX (study area) using a fine 3D rock mechanical parameter model based on corner-point grids. We converted the grid data and obtained a 3D finite element stress model of this area using hexahedral grid block through constraint inversion and finite element stress simulation. Then, we converted the model back and obtained a stress field model of the area based on corner-point grids. Against this model, we validated the model in terms of the magnitude and direction of the simulated stress.The tight sand reservoir in regionX lies in a city in Xinjiang, China, where the oil-bearing area is 32.4 km2. The location is an ultralow-abundance, medium-productivity, intermediate-deep medium-sized oil reservoir with a geological reserve of 100,954,500 barrel oil and recoverable reserve of 10,954,500 barrel oil. The reservoir mainly comprises gray and brown fine conglomerate, with lesser amounts of grayish green pebbly inequigranular lithic sandstone and pebbly intermediate-coarse lithic sandstone of varying sizes. The oil layers stand 100-190 m thick, with the porosity ranging from 1.9% to 16.8% with an average of 8.48% and the permeability spanning from 0.014 to 464 mD with an average of 5.54 mD. The oil layers are widely distributed and highly heterogeneous. The case study model based on corner-point grid consists 3960000 grid blocks and 4139486 nodes. The grid size is 25 m both in x and y direction, and the grid block’s height varies from 4.93 m to 16.61 m. Utilizing an in-house developed software platform, we established a fine corner-point grid rock mechanical parameter model of this area (Figures 11(a) and 11(b)) using collocated sequential Gaussian simulation; the stress is then simulated using FEM simulator, ABAQUS. According to the magnitude and direction of the principal stress in the study area, we established maximum horizontal principal stress and minimum horizontal stress models (Figures 11(c) and 11(d)) and a stress direction model (Figures 11(e) and 11(f)) of the area.Figure 11
3D corner-point grid attribute model for regionX: (a) Poisson ratio, (b) Young’s modulus, (c) minimum horizontal principal stress, (d) maximum horizontal principal stress, (e) overall perspective of the maximum principal stress direction, and (f) side view of the maximum principal stress direction.
(a)
(b)
(c)
(d)
(e)
(f)Using the measured minimum principal stresses of seven samples from regionX and the paleomagnetic-based maximum principal stress direction data at three observation points, we extracted the corresponding grids in the simulated stress fields. By comparing the simulated result with the measured result (Tables 1 and 2), we discover a modest relative error between them, with both the positive and negative deviations within 5%. The results show a high agreement between the simulated result and measured result, validating that our new algorithm is both effective and accurate. From the results, the simulated values well describe the magnitude of 6 real samples and direction of 3 test points from region X and our method for the finite element simulation of in situ stress based on the 3D corner-point grid model can accurately simulate the stress distribution and it well reproduces the values connected to the real samples.Table 1
Minimum horizontal principal stresses: measured versus simulated results.
Sample ID
Minimum horizontal principal stress
Error (%)
Measured (MPa)
Simulated (MPa)
1
56.86
54.84
3.56
2
56.08
53.34
4.88
3
57.26
55.21
3.57
4
58.43
60.32
3.23
5
59.32
62.12
4.73
6
57.94
60.43
4.29Table 2
Maximum horizontal principal stress direction: measured versus simulated results.
Observation point
Maximum horizontal principal stress direction
Measured (°)
Mean (°)
Simulated (°)
Error (%)
1
32.61
34.71
35.62
2.62
36.81
35.91
33.81
34.41
2
33.21
32.93
34.43
4.56
32.81
32.91
32.91
32.91
32.81
3
31.71
32.93
34.73
5.47
32.91
37.21
33.61
34.41
27.71The error between the simulated result and the measured result, however, may be attributable to the following facts: (1) the 3D Poisson ratio model and Young modulus model of the study area are based on logging data and spatial interpolations, and the accuracies of these models directly affect the accuracy of the material parameters of the finite element grid model and consequently the accuracy of the simulated result. In the case study, we applied collocated sequential Gaussian simulation approach with both logging data on the wells and prestack inversion seismic data between wells to reduce the error in modeling; meanwhile, several geologists and engineers are employed to carry out the quality control process to ensure the models’ accuracy. Still, the approach may not exactly represent the attributes distribution due to the uncertainty in geology; (2) the measured in situ stresses are the stress levels for a particular depth, whereas the simulated values are the stress data of the grid corresponding to this depth. As the grid data represent the mean stress in a range around this point, error is practically inevitable; (3) by defining proper boundary loads and boundary constraints, it is possible to obtain a 3D stress model of the area. Different combination of the boundary load data and differences in the boundary grid constraints can result in errors between the simulated result and measured result; (4) the test equipment may also contribute to the errors between the simulated result and measured result.


## 6. Conclusions and Discussion

In this paper, we present a finite element in situ stress simulation method based on the 3D corner-point grid model by designing a simulation procedure, analyzing the grid structure, and implementing related algorithms. Using this method, we were able to both retain 3D fine geological models and provide an accurate presentation of the 3D stress distribution. Using data of theX oil field as an example, we validated the effectiveness and accuracy of the proposed method. From our investigation, the following conclusions can be drawn:(1)
Using a grid conversion algorithm based on the data structures of both corner-point grids and finite element grids, it is possible to convert between these two types of grids while retaining the grid properties; the attribute information of corner-point grids can also be accurately resampled into the corresponding finite element grids to allow the automatic definition of the corresponding material attribute parameters.(2)
The design and implementation of the grid conversion algorithm contributes to establishing a finite element in situ stress simulation method based on corner-point grids. This method not only retains the reservoir structure properties based on the corner-point grid model but also allows the accurate acquisition of the in situ stress distribution in the target area through FEM. The stress distribution in the target area was simulated based on the measured result of the oil fields in regionX, and the simulated result was compared with the measured result to validate the effectiveness and accuracy of the proposed method.(3)
The implementation of the finite element in situ stress simulation method based on corner-point grids offers important references for sweet spot evaluation, fracture morphology, and distribution simulation and well pattern optimization and guarantees the procedure continuity and data consistency in modeling operations for the exploration and development of oil and gas fields.



---

*Source: 7384085-2020-02-25.xml*