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

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. .e 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. .e stress model is then resampled to corresponding corner-point grid geological models using the reverse algorithm. .e 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, the X region. Comparing the simulated result with the measured result, the high agreement validated the effectiveness and accuracy of the proposed method.


Introduction
e 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]. e 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][3][4]. e 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]. e 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][9][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][14][15]. e 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][21][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.
is 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. us, the finite element simulation of in situ stress based on a 3D cornerpoint 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. en, 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 region X and compared our simulated result with the measured result. e result validates the effectiveness and accuracy of the proposed method.

Overall Procedural Design
e 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. en, 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. e conversion algorithm extracts a series of information of all nodes of the cornerpoint model, including node number, node coordinate, and node attribute. e mapping relationships of nodes and grid blocks are also extracted and saved. e 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 cornerpoint grid block attributes and finite element grid block material properties has been established. e 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.
us, the corner-point grid is converted into a finite element grid. e 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. is conversion algorithm does not introduce any instability and reduces accuracy result in the model. e detailed reverse algorithm is almost similar to that of the conversion from corner-point grid to finite grid while it is an opposite process. e in situ stress attributes including magnitude and direction are also linked to each grid blocks. e overall procedure for the design is shown in Figure 1.
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. e 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 gridbased 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]. e current existing FEM simulator cannot directly set up a fine heterogeneous attribute model-based geological statistics methods to represent heterogeneous attribute distribution [30,32]. is 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][9][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. e 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 cornerpoint grid geological models. e corner-point gridbased 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.

Corner-Point Grid.
e 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-andnonstructured 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]. e 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].
e 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 Logically, there are (N x + 1), (N y + 1), and (N z + 1) coordinate lines in the X, Y, and Z (depth) directions, respectively, and the grid model is divided into N x , N y , and N z grid cells:

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. e 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. e 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).
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. Figure 3(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.

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 cornerpoint grid, mainly contains the X, 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). e following functions are contained in the bidirectional grid conversion algorithm ( Figure 4): 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 Mathematical Problems in Engineering 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.

Preprocessing Procedure. First, PreProcessAlgorithm() is
used to check the validity of the input model. en, 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 Figure 5.
We tested the 3D corner-point grids without inactive grids (grids are useful in data organization while are not   Material-61 -material-80 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. 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.
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 Figure 7. 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 N x , N y , and N z 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. us, 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.
ere are two different case models in Figure 8, 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 (

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 Gen-erateDisplacement()). en, 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, Post-ProcessAlgorithm() 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 Figure 9.
We established a finite element grid model using the preprocessing procedure. en, 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 Figure 10, with a case study. e 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. e values of every grid blocks are the same in both finite element grid-based stress model and corner-point grid-based stress model.

A Case Study of Region X
In this section, we tested our algorithm on real oil field data of region X (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   (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.
Using the measured minimum principal stresses of seven samples from region X 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%. e 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. e 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.

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 the X 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) e design and implementation of the grid conversion algorithm contributes to establishing a finite element in situ stress simulation method based on corner-point grids. is method not only retains the reservoir structure properties based on the cornerpoint grid model but also allows the accurate acquisition of the in situ stress distribution in the target area through FEM. e stress distribution in the target area was simulated based on the measured result of the oil fields in region X, and the simulated result was compared with the measured result to validate the effectiveness and accuracy of the proposed method.
(3) e 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.

Data Availability
All the data is within this paper and supporting files.

Conflicts of Interest
e authors declare that they have no conflicts of interest.