- Email: [email protected]
Development of a flow based dynamic gridding approach for fluid flow modeling in heterogeneous reservoirs
Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Development of a flow based dynamic gridding approach for fluid flow modeling in heterogeneous reservoirs Mahdi Amirsardari a, Bahram Dabir a, b, c, *, Abbas Naderifar b a
Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Iran Department of Chemical Engineering, Amirkabir University of Technology, Tehran, Iran c Energy Research Center, Amirkabir University of Technology, Tehran, Iran b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 11 December 2015 Received in revised form 23 March 2016 Accepted 24 March 2016 Available online 28 March 2016
The non-uniform flow based gridding is a well-known approach to create coarse-scale dynamic model to obtain results with high accuracy and low computation cost. In this approach, the grid refinement is performed only in the flow feature regions. Since the flow features will change during the reservoir life, the dynamic gridding method makes it possible to capture the accuracy. In this study, the combination of the grid generation and the reservoir simulation sciences has been used to create the coarse-scale dynamic mesh model with high accuracy and low simulation running time. In this model, the mesh will be updated based on the dynamic interest regions during the simulation time. The simulation results of the presented methodology have been compared with the conventional static gridding models and fine-scale models for oilewater and gasewater flows in the heterogeneous oil and gas reservoirs. © 2016 Elsevier B.V. All rights reserved.
Keywords: Dynamic mesh Heterogeneous reservoir Reservoir simulation Non-uniform coarsening Flow based gridding
1. Introduction Dynamic reservoir simulation models are used to compare the different development scenarios and to predict the hydrocarbon recovery on the basis of a given production strategy. The numerical dynamic simulation modeling is the final stage in the field development planning process. The reservoir modeling section in the field development planning process includes the three important stages: creating a fine-grid geological model (static model), upscaling and coarsening of static model to dynamic model and history matching. The upscaling technique can be defined as assigning the effective parameters from a fine-grid model with a large number of elements to a coarse-grid model such that the simulation results are similar in both models (Portella and Hewett, 1998). It is clear that the coarse-scale grids must be constructed before the upscaling process. Usually coarse-grid model is formed from joining several adjacent elements in fine-grid model (King et al., 2005). This process is called Coarsening. Clearly, the upscaling process is closely related to the generation of the coarse
* Corresponding author. Chemical Engineering Department, Amirkabir University of Technology, No.424, Hafez Ave., P. O. Box 15875-4413, Tehran, Iran. Tel.: þ98 9121304379. E-mail address: [email protected] (B. Dabir). http://dx.doi.org/10.1016/j.jngse.2016.03.077 1875-5100/© 2016 Elsevier B.V. All rights reserved.
elements. Dynamic model is a fine-grid geological model under the coarsening process in which the effective parameters have been assigned from a fine-grid model to a coarse-grid model (Chen, 2009). In the dynamic model, the running time of the computer simulation for the solution of the flow equations considerably decreases in comparison to the fine-grid model. After the history matching process, this model is used for the simulation of the different field development scenarios. The generation of the coarse elements from fine-grid model is a key factor in the coarsening process which is affecting the upscaling and the solution accuracy. Obviously, an ideal coarsening should be implemented in such a way that the simulation results of a coarse-grid model with minimum number of the elements be similar to a fine-grid model. In an ideal coarsening, the coarse elements are not uniform in the computational domain. In other words the elements have different sizes and their sizes depend on some dynamic parameters. The non-uniform coarsening process generates an irregular mesh. Defined uniform coarsening processes in the commercial software are not optimized. Many researchers have interested in proposing automatic upscaling and coarsening methods to generate a coarsegrid model with the minimum number of elements and high accuracy in the simulation results. To achieve this purpose, the combination of the grid generation and the reservoir simulation sciences will be helpful. Various methods for the discretization of
716
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 1. The algorithm of dynamic reservoir simulation using dynamic mesh.
Fig. 3. The function for computation of element size from ESI.
Fig. 2. ESI based on well effect.
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
717
Fig. 5. Iterative process in Bossen method.
Fig. 4. Initial mesh based on well effect.
the computational domain have been presented by many researchers. Among them, Flow Based Gridding method has the particular popularity (Aarnes et al., 2009). Some researchers developed non-uniform structured coarsegrid models by solving the fluid flow equations on the fine-scale model and the use of the streamline method (Verma and Aziz, mez-Herna ndez, 1997; Portella and Hewett, 1998; Wen and Go 1998). Due to using structured gridding concept in these studies, orthogonality will not be preserved. Garcia et al. (Garcia et al., 1992) proposed a method for the construction of the structured mesh which has been called elastic mesh. In this approach an element in the coarse-scale model was constructed from the joining of elements with the similar permeability in the fine-scale model. This idea has been used by other researchers (Tran, 1995; Farmer, 2002; Wen et al., 2003; Ebrahimi and Sahimi, 2004). For the first time, Heinemann et al. (Heinemann et al., 1991; Heinemann and Brand, 1989) developed Voronoi or PEBI (perpendicular bisection) grids to simulate fluid flow in porous media. Mlacnik et al. (Mlacnik et al., 2003) used PEBI window method for flow based unstructured gridding. In this method the pressure solution on the fine-scale model and the streamlines have been used to generate the coarse-scale model. Then, a hybrid procedure with the gradient and Laplacian smoothing steps was applied on the coarse-scale model. Palagi (Palagi, 1992) presented a twodimensional modular gridding approach to generate PEBI grids. In this approach, the geological features such as faults were considered in the mesh construction process. Cirpka et al. (Cirpka et al., 1999) and Cao and Kitanidis (Cao and Kitanidis, 1999) used flow based gridding approach to solve the contaminate transport and the pressure equations in the hydrology context. Cao and Kitanidis (Cao and Kitanidis, 1999) considered a posteriori error estimator as a criterion for the mesh refinement. The construction of a background grid based on the velocity and the permeability has been proposed for unstructured coarse-grid generation (Evazi and Mahani, 2010). Amirsardari et al. (Amirsardari et al., 2015) developed a non-uniform coarsening method based on the anisotropic elements generation in the high permeability regions. They used static gridding approach to build the dynamic coarse-scale model.
In this study, in order to capture high gradients of the dynamic parameters in a reservoir, non-uniform unstructured mesh is used to minimize simulation running time and to obtain accurate results. Dynamic parameters such as saturation and pressure change during the reservoir life especially in water and gas injection cases. Adaptive Dynamic Unstructured meshes make it possible to obtain the accurate results with the minimum computation time. In this approach the structure of the mesh is updated on the basis of the steep gradients of the dynamic parameters. The non-uniform mesh is used to capture the phase interfaces with the small elements. In this method, the fluid flow solver and the mesh generator are coupled together. The structure of this paper will be discussed in the followings. In the next section, the algorithm of the proposed method will be described. The main steps of the algorithm are the construction of the initial mesh based on the well effects; solving the fluid flow equations and the grid generator. These main components of the method will be explained in detail in the different sections. Finally, in order to investigate the accuracy of the proposed method, twophase flow simulation results for the oilewater and gasewater flow cases in the heterogeneous oil and gas reservoirs are compared to the fine-grid model and uniform coarse-grid static mesh model.
2. Methodology The algorithm of the proposed method is shown in Fig. 1. The input data are the reservoir boundary and the wells location. In the first step the initial mesh will be constructed by the mesh generator which will be discussed in the following sections. The elements are small near wells in the initial mesh. The initial conditions (pressure and saturation) are assigned to each node of the mesh. The second step is the solution of the fluid flow equations to determine the pressure and the saturation of each node in the first time step. The important point in this algorithm is the interpolation error of the saturation parameter. The curvature of this flow variable controls the element size in the mesh. If the saturation curvature increases, a higher resolution or smaller elements will be needed to improve the accuracy of the solution. Based on the saturation parameter, the interpolation error is computed. In case of this error be greater than a certain threshold value, the initial mesh will be adapted based on the saturation parameter. The adaptation process includes computation of the metric tensor (element size) and the generation of the mesh. The metric tensor is the combination of Hessian matrix of the saturation variable and the acceptable interpolation error. The metric tensor controls the size of the elements in the mesh. The powerful mesh generator is developed to adapt the mesh based on the computed metric tensor. In the mesh adaptation process, some operations such as the node insertion and deletion, smoothing and edge swapping will be performed. The adaptation process will be repeated until the threshold value is
718
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 6. The algorithm of Bossen Method.
met. Then the pressure and saturation of the previous time step should be interpolated for the new mesh as initial conditions of the new time step. After interpolation the fluid flow equation should be solved for the new time step according to the new adapted mesh. This process is repeated for all the time steps during the reservoir simulation. 3. Initial mesh generation based on well effect The near-well velocity depends on the inertia forces and has a high gradient in this region. Therefore smaller elements are required near the wells to capture the high gradient velocity. For this purpose, we define Element Size Indicator (ESI) concept which is maximum in production and injection wells and then smoothly
vanishes with increasing the distance from the wells. The mathematical function that provides ESI can be shown as follows.
ESI ¼ expðbdÞ
(1)
where, ESI is Element Size Indicator, d is the distance from the well and b is the constant that controls damping intensity from the well. Fig. 2 shows ESI based on the well effects for a reservoir with two wells; injection well in bottom-left and production well in upperright of the reservoir. The ESI is converted into the element size function by assigning the desired minimum element size to maximum ESI and the maximum element size to minimum ESI (Fig. 3). In order to build the initial mesh, the resulted element size function and the reservoir boundary are imported to the mesh
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
719
Fig. 8. A connection in data structure.
Fig. 7. Bossen's smoothing function.
generator which will be described in the section 5. Fig. 4 shows the resulted initial mesh for ESI map in Fig. 2. 4. Fluid flow modeling Control Volume Finite Approximation (CVFA) method is used to discretize the governing equations in porous media (Li et al., 2003). The fully implicit approach is considered to discretize equations in time. The combination of Darcy's law and the individual phase conservation law are employed to describe the flow of wetting (w) and non-wetting (nw) phases in the porous media:
vðfraS Sa =Ba Þ r ¼ V: aS ua þ qa vt Ba kra ua ¼ kVpa
ma
a ¼ w; nw
a ¼ w; nw
Fig. 9. A node in data structure.
(2)
(3)
Sw þ Snw ¼ 1
(4)
Pc ¼ pnw pw
(5)
where, f is the porosity of the porous media, Pc is the capillary pressure, Sa is the saturation of each phase, Ba is the formation volume factor of each phase and raS is the density of each phase at the standard condition, k is the absolute permeability tensor; kra, pa and ma are the relative permeability, pressure and viscosity of each phase. For each node in the discretized domain, the initial pressure and the water saturation are defined. For the outer boundary (reservoir boundary), no-flow boundary condition and for the inner boundary, constant bottom-hole pressure and constant injection flow rate are considered.
all the element sizes will be adapted to the element size function (Fig. 5). The input functions to the grid generator are the boundary of the computational domain and an element size function which is defined on the fine-scale cartesian mesh. First Delaunay Triangulation is implemented in the computational domain. Then, in each step of the iterative process, a node is selected randomly and its position is smoothed according to its neighbors. Finally, Delaunay Triangulation condition is applied and then by considering the nodes density, a node will be inserted or deleted. This process will be repeated until final mesh achieved. Fig. 6 shows the mentioned algorithm details that will be elaborated in the followings.
5. Mesh generation algorithm In this study, Bossen method (Bossen, 1996) in accompany with some modifications which are applied to this method is used to create the dynamic mesh. In this section after reviewing Bossen method, the proposed modifications to decrease the gridding time will be described. 5.1. Review of Bossen method In this method, first an initial triangulation is formed then retriangulation, fining, coarsening and smoothing processes in an iterative loop lead the initial triangulation to the final case so that
Fig. 10. Initial polygon and relation between connections.
720
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
function has been suggested as below.
4 f ðdÞ ¼ 1 d4 :ed
(9)
This smoothing function exerts the attractive and the repulsive forces on the nodes same as the intermolecular forces. If the normalized length is less than one then the nodes will attract each other whereas if the normalized length is greater than one then the nodes will repel each other. In addition, unlike the intermolecular forces, the forces values among the nodes in which have the distances close to zero will have limited reduction values. Fig. 7 shows the applied forces between the two nodes versus the normalized length. Also the neighboring points of node i, were defined as the points had connections with point i and the smoothing constant is 0.2. Due to the triangular structure maintenance during the mesh generation process, the new position of point i will only be acceptable when the triangular 0 structure is n't violated. If the segment between pi and pi does not 0 cross any existing edge then pi node position will change into pi . Otherwise, if the segment crosses a border connection then the node will be transferred to the junction location of the border connection. If none of the cases happen, then the node is n't moved from its location.
Fig. 11. Gradual triangulation and connection updating.
5.1.2. Inserting or deleting node As shown in Fig. 6, if in the computational mesh generation, the mesh is sparse or dense in a location then some nodes will be inserted in or deleted from the mesh. In order to express a quantitative measurement of the sparseness for each node or connection a quantitative parameter which is called Extent is defined. The extent has a normalized area unit. The normalized area parameter is similar to the normalized length in the surface and it can be calculated from the area divided by the desired area (r2) which is computed from the element size function. The following definitions have been presented for the extent of the nodes and the connections.
Fig. 12. Completion of polygon triangulation and connection relation statues.
5.1.1. Smoothing In order to generate an isotropic mesh by considering a specified element size function, the nodes distances can be determined using the following equation. In this equation, k k is a norm operator, and r is the desired distance which can be accessible using an element size function.
dðp1 ; p2 Þ ¼
kp1 p2 k r
(a) Border Connection Extent: half the square of its normalized length (b) Border Node Extent: half the square of the two neighboring connections normalized lengths sum (c) Inner Connection Extent: 0.8 times two adjacent triangular normalized areas pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (d) Inner Node Extent: 2:3= nðn 2Þ times the n surrounding triangular normalized areas sum In the algorithm which is shown in Fig. 6, if i is the selected node
(6)
The proposed equation can be called Normalized Length. Based on this definition, the normalized length of the desired mesh is 1 and the smoothing can be determined using the following first order motion equation. 0
pi ¼ pi þ ai
X f d pi ; pj :uij ;
(7)
j2Ni
pi pj uij ¼ d pi pj
(8)
In the Eq. (7) and Eq. (8), pi is the spatial position of node i, uij is a unit normal vector, Ni is the set of node i adjacent nodes, f is a smoothing function, ai is a smoothing constant. Bossen smoothing
Fig. 13. Initial guess in one-dimensional gridding.
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
721
assumed active. After smoothing pi node and by considering the 0 displacement value of pi pi , the node will be considered as active or inactive. If this displacement value exceeds a threshold value, then the node will be remained active and adjacent nodes will be assumed active. Otherwise it is considered as inactive. In case of the node insertion or deletion, its surrounding nodes will be active. If all nodes be inactive, the algorithm process will be stopped. 5.2. Implementation of Bossen method
Fig. 14. An Element size function and a line segment.
in the algorithm, its extent will be calculated and if its extent is less than 1 then this node will be deleted and the mesh will be locally re-triangulated. If node i is not deleted then, its surrounding connection extent will be calculated. If the extents largest value is greater than 1, its relevant connection will be split through the middle and the node will be added to the mesh and others remained re-triangulated using Delaunay method. The initial border points are never deleted or repositioned. It should be noticed; re-triangulation processes in case of the node insertion will be performed using Delaunay gradual triangulation method. In case of node deletion, the simple achieved polygon is triangulated and by swapping the connections in the polygon, Delaunay triangulation condition will be applied.
5.1.3. Convergence The entire process of smoothing, triangulation, deletion and insertion node, lead all the normalized lengths of the connections to 1. In order to assess the mesh convergence statues, entire node classified as active or inactive nodes. At first all the nodes are
The object-oriented programming and Cþþ language program are used for mesh generation process. A new data structure is used to implement this method. In this data structure, connections are considered as the main components. These components control the data relationship. In each connection two pointers to two nodes at both ends of the connection and two pointers to two adjacent connections exist with assuming counter-clockwise rotation around the mentioned nodes. In order to complete the data structure and retrieving the nodes location, at each node only a pointer to a connection is considered. Accordingly, there is no single node in data structure and each node is connected to a connection (Figs. 8 and 9). Based on the mentioned data structure, in the first step, connections in the polygon which are defining the boundaries of computational domain are established. Then an initial triangulation within the original polygon is formed by simple triangulation method. The important item at this stage is the updating of the data structure gradually along with the completion of the initial triangulation. After the initial triangulation, the data structure between nodes and connections will be completed. The data structure will be updated by adding or removing nodes. Fig. 10 shows the input polygon and the relation between connections at the initial state. Fig. 11 shows the gradual triangulation and the updating data structure and Fig. 12 shows the completion of the initial triangulation and the relationship between connections. It should be noticed when removing a node, a simple polygon is formed and so, the data structure updating is a similar process to the initial start status. Advancing Front method and Delaunay condition are also used to Delaunay triangulation in a simple polygon. 5.3. Modification on Bossen method In order to improve the efficiency of Bossen method, two
Fig. 15. Initial one-dimensional gridding and final gridding after smoothing.
722
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 16. First stage algorithm in two-stage gridding.
modifications are applied to the mentioned method. The first is the placement of the nodes on the edges of the polygon (boundary nodes) based on the element size function and separation of this stage as one-dimensional gridding. The second modification is the implementation of the original algorithm (Fig. 6) in two stages to speed up the grid generation process. In the first stage an estimation of final mesh will be obtained quickly by limited smoothing and in the second stage final smoothing will be performed. In the two following sections these variations are described.
5.3.1. Boundary gridding separation of inner gridding After importing computational domain boundaries as a polygon, each edge is considered as one-dimensional space and the position of the boundary nodes will be fixed based on the element size function on the edges. As aforementioned, in the proposed algorithm in Figs. 10e12, a triangle structure will be preserved during mesh generation. By applying mentioned modification, we will only check the condition of the existence of the node position insides it's surrounded polygon in the smoothing process and it is not
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
723
Fig. 17. Second stage algorithm in two-stage gridding.
0
required to find a possible collision line segment between pi and pi with polygon boundaries surrounding it (pi is node location before 0 smoothing and pi is node location after smoothing). Also, it doesn't have the necessity to check the crossed edge whether is a boundary edge or not. In order to perform one-dimensional gridding which leads to the determination of the boundary nodes, the node position of the two ends of the segment will be fixed and then an initial guess of the intermediate nodes position estimated. In the next step, in an iterative process including node removing and smoothing, the initial guess will be optimized based on the element
size function. To estimate the initial guess, a simple algorithm is used. First the desired lengths based on the element size function in the boundary nodes are computed. If the summation of these lengths half is less than the segment length of the two boundary nodes, then a node will be added in the middle of the free zone. Free zone is defined as the region between the segment of the two boundary nodes and out of the desired lengths half. Fig. 13 shows a schematic of a step in estimating the initial guess. In this figure d1 and d2 are desired lengths based on the element size function for both ends of the
724
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 18. Boundary of computational space and initial guess of mesh.
Fig. 21. Permeability map of a heterogeneous reservoir. Fig. 19. Final mesh based on the element size function.
with the largest absolute deviation from the normal lengths will be identified and removed. Then smoothing is performed again on the nodes and the total absolute deviation from the normal lengths of the nodes will be calculated. If the absolute deviation is decreased, the last position of the nodes will be stored. This process continues to minimize the absolute deviation from the normal length. Fig. 14 shows an element size function and a line segment in the
Table 1 Water relative permeability in oilewater system.
Fig. 20. An ellipse for metric tensor definition.
segment. Above method is repeated in an iterative process for the new segments until no free zone between the nodes exists. After performing the initial guess, smoothing conducted on the intermediate nodes simultaneously. The position of the nodes will be stored after the initial smoothing. The total absolute deviation from the normal lengths of the nodes is calculated and the nodes
Sw
krw
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.8 0.9 1
0 0.0028 0.025 0.07 0.14 0.22 0.33 0.62 0.80 1
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729 Table 2 Oil relative permeability.
725
Table 3 Simulation running times for three models in oilewater flow example.
So
kro
Model
Run time (second)
0 0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9 0.95
0 0.011 0.044 0.177 0.277 0.399 0.543 0.709 0.897 1
Fine-grid Uniform isotropic static grid Non-uniform dynamic grid
21,282 380 378
computational space. Fig. 15 shows the initial guess grid and final grid using the mentioned algorithm. 5.3.2. Two-stage gridding Boundary gridding separation of the inner gridding in the computational domain provides a more coherent smoothing process. The smoothing in terms of the computation time is the most costly step in the algorithm shown in Fig. 6. Therefore, we decided to conduct the smoothing process in the two stages. First it is attempted to provide a primary estimation of the nodes position in the mesh and then primary mesh will be lead to the final mesh. Smoothing of each node position based on the first order motion equation is performed in several steps. The smoothing process is as follows. f
pi /smoothing/p1i /smoothing/p2i …/pi
(10)
where, pi is the first position of node i. This position is transferred to the next position during the smoothing process and finally converges to pfi position. If the distance between the two nodes in the two consecutive steps is defined as Dn pi ¼ pnþ1 pni , then by i using a suitable smoothing function, it is expected that the following relationship will be established.
D0 > D1 > … > Df 1
after the first step in the smoothing process is considered as the final position (pfi ¼ p1i ). Using the method in the algorithm which is shown in Fig. 6, an estimation of the final mesh is achieved quickly. After determining the initial position of the nodes, the final smoothing process will be conducted. Now, all the nodes except the boundary nodes are considered as active nodes and smoothing is conducted on all the active nodes simultaneously. After smoothing operation, Delaunay Triangulation condition is applied to the grid connections. Then the extent of all the active nodes will be examined. According to the computed extents, nodes may become inactive, be removed or one node be inserted. These steps are continued until all nodes in the mesh will become inactive. The steps to remove, deactivate, activate and insertion of the new node in the mesh are conducted according to the algorithm which are shown in Figs. 16 and 17. These figures show the new algorithm based on the two-step gridding. An example of the element size map, boundary of the computational domain and the initial estimation of final mesh are depicted in Fig. 18. Fig. 19 shows the final mesh in accordance with the element size function.
5.4. Element size function computation The Metric Tensor which controls size, shape and orientation of the elements is required for the gridding. The metric tensor will be defined as an ellipse with the two radiuses r1, r2 and the angle q (Fig. 20). The corresponding metric tensor is (Bossen, 1996):
(11)
According to a threshold value the smoothing process is stopped and pfi is determined. In the initial estimation of the node position, a limited smoothing is performed. It means that the node position
M¼ ¼Q
cos q sin q
g1 0
sin q cos q 0 QT
g2
g1 0
0
g2
cos q sin q
sin q cos q
(12)
where, g1 ¼ 1=r12 and g2 ¼ 1=r22 are eigenvalues of metric tensor and Q is the rotation matrix of angle q. The following Metric Tensor in two-dimension is computed as (Amirsardari et al., 2015; Hecht, 1998):
M¼
1 jHðuÞj ε0 supðuÞ infðuÞ
(13)
where, ε0, is user defined interpolation error and H(u) is Hessian matrix of water saturation function:
2
v2 Sw
6 6 vx2 HðuÞ ¼ 6 6 2 4 v Sw vyvx
Fig. 22. Water Saturation distribution and dynamic mesh at 40th day.
3 v2 Sw 7 vxvy 7 7 7 v2 Sw 5
(14)
vy2
In order to control the minimum and maximum size of the elements in the mesh, these parameters should be defined (CastroDiaz et al., 1997). For this purpose, gi(i ¼ 1,2) of the metric in Eq. (13) is limited as follows:
726
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 23. Water-cut for fine-scale, uniform coarse-scale and proposed dynamic mesh models.
Fig. 24. Water saturation distribution at 57th day for fine-scale model (right) and uniform coarse-scale static mesh model (left).
6. Results and discussions
1
g~ 1;2 ¼ min max
g1;2
; 2
hmax
! ;
1 h2min
! (15)
where, hmin and hmax are the minimum and maximum edge lengths. Due to the usage of isotropic elements in the study, the largest eigenvalue is considered for both diagonal eigenvalues. So, the isotropic metric is defined as below:
M¼
~2Þ ~1; g maxðg 0
0 ~2 Þ ~1; g maxðg
(16)
In order to validate the proposed method in the study, the twophase flow simulation in the oil and gas reservoirs was performed for the three models: the fine-grid geological model, the uniform coarse-grid model with static gridding and the non-uniform adaptive coarse-grid model with dynamic gridding. The simulation results of the models are used to determine the efficiency and the accuracy of the proposed method in the study. For this purpose, the fine-grid model was considered as a reference model and the simulation running time of the uniform static gridding model is in the same level of the dynamic gridding approach. We present two examples of oilewater and gasewater flows in the oil and gas reservoirs, respectively. The permeability map is shown in Fig. 21. In
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
727
Fig. 25. Water saturation distribution at 57th day for fine-scale model (right) and non-uniform coarse-scale dynamic mesh model (left).
Table 4 Water relative permeability in gasewater system. Sw
krw
0.22 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.01 0.03 0.07 0.15 0.24 0.33 0.65 1
Table 5 Gas relative permeability. Sg
krg
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.78 1
0.04 0.1 0.2 0.3 0.399 0.42 0.62 0.81 1
both examples the reservoir porosity is 0.1 and at the initial state the water saturation is 0. 6.1. Example 1: oilewater flow in an oil reservoir We consider oilewater flow simulation in this example. The initial reservoir pressure is 9014.7 psi. The water and oil relative permeability are shown in Table 1 and Table 2, respectively. A water injection well (Inj0) and a production well (Prod0) considered in the model in the bottom-left and upper-right locations, respectively (Fig. 21). The production well produces with the constant bottom-hole pressure of 5000 psia and in the injection well water is injected with the constant water rate of 450 BBL/D.
Fig. 26. Water saturation distribution and dynamic mesh at 28th day in gasewater flow.
Table 6 Simulation running times for three models in gasewater flow example. Model
Run time (second)
Fine-grid Uniform isotropic static grid Non-uniform dynamic grid
20,878 350 338
The proposed dynamic model approach is used to simulate the oilewater flow in the reservoir. The maximum number of the elements in the model is 2528. The generated adapted mesh at 40th day is shown in Fig. 22. As it can be observed in this figure, smaller elements in the regions with high saturation gradients and near wells were generated using the proposed algorithm. To validate the proposed method, the simulation results were compared to the fine-grid model and the uniform coarse-grid model (static mesh for both models) with 93,690 and 5894 elements, respectively. The number of the elements in the uniform
728
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 27. Water production rate for fine-scale, uniform coarse-scale and proposed dynamic mesh models at Prod 1.
Fig. 28. Water production rate for fine-scale, uniform coarse-scale and proposed dynamic mesh models at Prod 2.
coarse-grid model is determined in such a way that the simulation running times in both non-uniform coarse-grid dynamic mesh model and uniform coarse-grid static mesh model be at the same level (Table 3). Table 3 shows that the simulation running times for both coarse-scale models (static and dynamic mesh models) are in portion of 1/55 in comparison to the fine-grid geological model. Water cut (water rate to liquid rate) at the production well during the simulation time in three mentioned models has been shown in Fig. 23. Fig. 24 shows water distribution at 57th day (water breakthrough time in the reference model) for the uniform coarse scale static gridding model and fine-scale static mesh model (reference model). Also, water distribution at 57th day for the nonuniform coarse scale dynamic mesh model and fine-scale static
mesh model (reference model) have been shown in Fig. 25. As can be observed in Figs. 23e25, using dynamic mesh model the water cut and water breakthrough time can be obtained with the high accuracy in comparison to the uniform static mesh model (about 20% improvement in the water breakthrough time). 6.2. Example 2: gasewater flow in a gas reservoir In this example we simulate gasewater flow in a gas reservoir. The initial reservoir pressure is 5300 psi. The water and gas relative permeability are shown in Table 4 and Table 5, respectively. A water injection well (Inj1) and two production wells (Prod1 and Prod2) considered in the model (Fig. 21). The production wells
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
produce with the constant bottom-hole pressure of 3000 psia and in the injection well water is injected with the constant water rate of 800 BBL/D. The proposed dynamic model approach is used to simulate gasewater flow in the reservoir. The maximum number of the elements in the model is 2379. The generated adapted mesh at 28th day is shown in Fig. 26. As it can be observed in this figure, the smaller elements in the regions with the high saturation gradients and near wells were generated using the proposed algorithm. To validate the proposed method, the simulation results were compared to the fine-grid model and the uniform coarse-grid model (static mesh for both models) with 93,690 and 4848 elements, respectively. The number of the elements in the uniform coarse-grid model is determined in such a way that the simulation running times in both non-uniform coarse-grid dynamic mesh model and uniform coarse-grid static mesh model be at the same level (Table 6). The water rates for two production wells in three mentioned models have been shown in Fig. 27 and Fig. 28. As can be observed in these figures, using the dynamic mesh approach the water breakthrough time and the water production rate can be computed with the high accuracy. 7. Conclusions and recommendations A new automatic reservoir simulation approach based on the dynamic mesh concept has been introduced. The advantages of this methodology and the recommendations for the future works have been summarized in the following items: 1. One of the main issues in the use of the dynamic mesh is the simulation running time. In this study, the non-uniform adaptive coarse gridding is used to have the smaller elements only in the dynamic interest regions in order to optimize the simulation running time. 2. A flexible two-dimensional computation mesh generator is used for the dynamic gridding. In order to improve the efficiency and decrease the gridding time, some modifications have been applied to this mesh generator algorithm. 3. The creation of dynamic mesh is based on the saturation variations and the well effects. Any other parameters (static or dynamic) can be used to build the dynamic mesh for the future works. 4. The numerical tests of the presented methodology in the oil and gas reservoirs show that the use of the mentioned approach makes it possible to obtain the high accurate results with the low computation cost. 5. As all the mentioned steps in the proposed algorithm can be extended to 3D models, this method will be applicable in case of 3D modeling. Nomenclature M f pa S
a k kra
raS Ba
Metric tensor Smoothing function Pressure of phase a Saturation Phase (oil/water/gas) Permeability tensor Relative permeability of phase a Surface density of phase a Formation volume factor of phase a
Pc ua t
ma p H ESI ε0 d
f g hmin hmax sup inf
729
Capillary pressure Velocity of phase a Time Viscosity of phase a Spatial point Hessian matrix Element Size Indicator Interpolation error Distance Reservoir porosity Eigenvalue of metric tensor Minimum edge length Maximum edge length Supremum Infimum
References Aarnes, J.E., Lie, K.A., Kippe, V., Krogstad, S., 2009. Multiscale Modeling and Simulation in Science, vol. 66. Springer Verlag, p. 27. Amirsardari, M., Dabir, B., Naderifar, A., 2015. Anisotropic grid generation for dynamic simulation in heterogeneous hydrocarbon reservoirs. J. Nat. Gas. Sci. Eng. 27, 1523e1535. Bossen, F.J., 1996. Anisotropic Mesh Generation with Particles. Msc. Thesis. School of Computer Science Carnegie Mellon University. Cao, J., Kitanidis, P.K., 1999. Adaptive-grid simulation of groundwater flow in heterogeneous aquifers. Adv. Water Res. 22, 681e696. Castro-Diaz, M.J., Hecht, F., Mohammadi, B., Pironneau, O., 1997. Anisotropic unstructured mesh adaption for flow simulations. Int. J. Numer. Methods Fluids 25, 475e491. Chen, T., 2009. New Methods for Accurate Upscaling with Full-tensor Effects. Ph.D. Thesis. Stanford University, Stanford, CA. Cirpka, O.A., Frind, E.O., Helming, R., 1999. Streamline-oriented grid generation for transport modelling in two-dimensional domains including wells. Adv. Water Res. 22, 697e710. Ebrahimi, F., Sahimi, M., 2004. Multiresolution wavelet scale up of unstable miscible displacements in flow through heterogeneous porous media. Transp. Porous Media 57, 75e102. Evazi, M., Mahani, H., 2010. Unstructured-coarse-grid generation using background-grid approach. SPE J. 15, 326e340. Farmer, C., 2002. Upscaling: a review. Int. J. Numer. Methods Eng. 40, 63e78. Garcia, M., Journel, A.G., Aziz, K., 1992. Automatic grid generation for modeling reservoir heterogeneities. SPE Reserv. Eng. 7, 278e284. Heinemann, Z.E., Brand, C.V., 1989. Gridding techniques in reservoir simulation. In: Proceedings of the First and Second International Forum on Reservoir Simulation, Alpbach, Austria. Heinemann, Z.E., Brand, C.W., Munka, M., 1991. Modeling reservoir geometry with irregular grids. SPE Reserv. Eng. 6, 225e232. Hecht, F., 1998. Bidimensional Anisotropic Mesh Generator. Technical Report, INRIA, Rocquencourt. King, M.J., Burn, K.S., Wang, P., Muralidharan, V., Alvarado, F.E., Ma, X., DattaGupta, A., 2005. Optimal coarsening of 3D reservoir models for flow simulation. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Li, B., Chen, Z., Huan, G., 2003. Control volume function approximation methods and their applications to modeling porous media flow. Adv. Water Res. 26, 435e444. Mlacnik, M.J., Harrer, A.W., Heinemann, Z.E., 2003. Locally streamline-pressurepotential-based PEBI grids. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Palagi, C.L., 1992. Generation and Application of Voronoi Grids to Model Flow in Heterogeneous Reservoirs. Ph.D. Thesis. Stanford University, Stanford, CA. Portella, R., Hewett, T.A., 1998. Upscaling, gridding, and simulation using streamtubes. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Tran, T.T., 1995. Stochastic Simulation of Permeability Fields and Their Scale-up for Flow Modeling. Ph.D. thesis. Stanford University, Stanford, CA. Verma, S., Aziz, K., 1997. A control volume scheme for flexible grids in reservoir simulation. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. mez-Herna ndez, J.J., 1998. Upscaling hydraulic conductivities in Wen, X.-H., Go cross-bedded formations. Math. Geol. 30, 181e211. Wen, X., Durlofsky, L., Edwards, M., 2003. Upscaling of channel systems in two dimensions using flow-based grids. Transp. Porous Media 51, 343e366.
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Development of a flow based dynamic gridding approach for fluid flow modeling in heterogeneous reservoirs Mahdi Amirsardari a, Bahram Dabir a, b, c, *, Abbas Naderifar b a
Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Iran Department of Chemical Engineering, Amirkabir University of Technology, Tehran, Iran c Energy Research Center, Amirkabir University of Technology, Tehran, Iran b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 11 December 2015 Received in revised form 23 March 2016 Accepted 24 March 2016 Available online 28 March 2016
The non-uniform flow based gridding is a well-known approach to create coarse-scale dynamic model to obtain results with high accuracy and low computation cost. In this approach, the grid refinement is performed only in the flow feature regions. Since the flow features will change during the reservoir life, the dynamic gridding method makes it possible to capture the accuracy. In this study, the combination of the grid generation and the reservoir simulation sciences has been used to create the coarse-scale dynamic mesh model with high accuracy and low simulation running time. In this model, the mesh will be updated based on the dynamic interest regions during the simulation time. The simulation results of the presented methodology have been compared with the conventional static gridding models and fine-scale models for oilewater and gasewater flows in the heterogeneous oil and gas reservoirs. © 2016 Elsevier B.V. All rights reserved.
Keywords: Dynamic mesh Heterogeneous reservoir Reservoir simulation Non-uniform coarsening Flow based gridding
1. Introduction Dynamic reservoir simulation models are used to compare the different development scenarios and to predict the hydrocarbon recovery on the basis of a given production strategy. The numerical dynamic simulation modeling is the final stage in the field development planning process. The reservoir modeling section in the field development planning process includes the three important stages: creating a fine-grid geological model (static model), upscaling and coarsening of static model to dynamic model and history matching. The upscaling technique can be defined as assigning the effective parameters from a fine-grid model with a large number of elements to a coarse-grid model such that the simulation results are similar in both models (Portella and Hewett, 1998). It is clear that the coarse-scale grids must be constructed before the upscaling process. Usually coarse-grid model is formed from joining several adjacent elements in fine-grid model (King et al., 2005). This process is called Coarsening. Clearly, the upscaling process is closely related to the generation of the coarse
* Corresponding author. Chemical Engineering Department, Amirkabir University of Technology, No.424, Hafez Ave., P. O. Box 15875-4413, Tehran, Iran. Tel.: þ98 9121304379. E-mail address: [email protected] (B. Dabir). http://dx.doi.org/10.1016/j.jngse.2016.03.077 1875-5100/© 2016 Elsevier B.V. All rights reserved.
elements. Dynamic model is a fine-grid geological model under the coarsening process in which the effective parameters have been assigned from a fine-grid model to a coarse-grid model (Chen, 2009). In the dynamic model, the running time of the computer simulation for the solution of the flow equations considerably decreases in comparison to the fine-grid model. After the history matching process, this model is used for the simulation of the different field development scenarios. The generation of the coarse elements from fine-grid model is a key factor in the coarsening process which is affecting the upscaling and the solution accuracy. Obviously, an ideal coarsening should be implemented in such a way that the simulation results of a coarse-grid model with minimum number of the elements be similar to a fine-grid model. In an ideal coarsening, the coarse elements are not uniform in the computational domain. In other words the elements have different sizes and their sizes depend on some dynamic parameters. The non-uniform coarsening process generates an irregular mesh. Defined uniform coarsening processes in the commercial software are not optimized. Many researchers have interested in proposing automatic upscaling and coarsening methods to generate a coarsegrid model with the minimum number of elements and high accuracy in the simulation results. To achieve this purpose, the combination of the grid generation and the reservoir simulation sciences will be helpful. Various methods for the discretization of
716
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 1. The algorithm of dynamic reservoir simulation using dynamic mesh.
Fig. 3. The function for computation of element size from ESI.
Fig. 2. ESI based on well effect.
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
717
Fig. 5. Iterative process in Bossen method.
Fig. 4. Initial mesh based on well effect.
the computational domain have been presented by many researchers. Among them, Flow Based Gridding method has the particular popularity (Aarnes et al., 2009). Some researchers developed non-uniform structured coarsegrid models by solving the fluid flow equations on the fine-scale model and the use of the streamline method (Verma and Aziz, mez-Herna ndez, 1997; Portella and Hewett, 1998; Wen and Go 1998). Due to using structured gridding concept in these studies, orthogonality will not be preserved. Garcia et al. (Garcia et al., 1992) proposed a method for the construction of the structured mesh which has been called elastic mesh. In this approach an element in the coarse-scale model was constructed from the joining of elements with the similar permeability in the fine-scale model. This idea has been used by other researchers (Tran, 1995; Farmer, 2002; Wen et al., 2003; Ebrahimi and Sahimi, 2004). For the first time, Heinemann et al. (Heinemann et al., 1991; Heinemann and Brand, 1989) developed Voronoi or PEBI (perpendicular bisection) grids to simulate fluid flow in porous media. Mlacnik et al. (Mlacnik et al., 2003) used PEBI window method for flow based unstructured gridding. In this method the pressure solution on the fine-scale model and the streamlines have been used to generate the coarse-scale model. Then, a hybrid procedure with the gradient and Laplacian smoothing steps was applied on the coarse-scale model. Palagi (Palagi, 1992) presented a twodimensional modular gridding approach to generate PEBI grids. In this approach, the geological features such as faults were considered in the mesh construction process. Cirpka et al. (Cirpka et al., 1999) and Cao and Kitanidis (Cao and Kitanidis, 1999) used flow based gridding approach to solve the contaminate transport and the pressure equations in the hydrology context. Cao and Kitanidis (Cao and Kitanidis, 1999) considered a posteriori error estimator as a criterion for the mesh refinement. The construction of a background grid based on the velocity and the permeability has been proposed for unstructured coarse-grid generation (Evazi and Mahani, 2010). Amirsardari et al. (Amirsardari et al., 2015) developed a non-uniform coarsening method based on the anisotropic elements generation in the high permeability regions. They used static gridding approach to build the dynamic coarse-scale model.
In this study, in order to capture high gradients of the dynamic parameters in a reservoir, non-uniform unstructured mesh is used to minimize simulation running time and to obtain accurate results. Dynamic parameters such as saturation and pressure change during the reservoir life especially in water and gas injection cases. Adaptive Dynamic Unstructured meshes make it possible to obtain the accurate results with the minimum computation time. In this approach the structure of the mesh is updated on the basis of the steep gradients of the dynamic parameters. The non-uniform mesh is used to capture the phase interfaces with the small elements. In this method, the fluid flow solver and the mesh generator are coupled together. The structure of this paper will be discussed in the followings. In the next section, the algorithm of the proposed method will be described. The main steps of the algorithm are the construction of the initial mesh based on the well effects; solving the fluid flow equations and the grid generator. These main components of the method will be explained in detail in the different sections. Finally, in order to investigate the accuracy of the proposed method, twophase flow simulation results for the oilewater and gasewater flow cases in the heterogeneous oil and gas reservoirs are compared to the fine-grid model and uniform coarse-grid static mesh model.
2. Methodology The algorithm of the proposed method is shown in Fig. 1. The input data are the reservoir boundary and the wells location. In the first step the initial mesh will be constructed by the mesh generator which will be discussed in the following sections. The elements are small near wells in the initial mesh. The initial conditions (pressure and saturation) are assigned to each node of the mesh. The second step is the solution of the fluid flow equations to determine the pressure and the saturation of each node in the first time step. The important point in this algorithm is the interpolation error of the saturation parameter. The curvature of this flow variable controls the element size in the mesh. If the saturation curvature increases, a higher resolution or smaller elements will be needed to improve the accuracy of the solution. Based on the saturation parameter, the interpolation error is computed. In case of this error be greater than a certain threshold value, the initial mesh will be adapted based on the saturation parameter. The adaptation process includes computation of the metric tensor (element size) and the generation of the mesh. The metric tensor is the combination of Hessian matrix of the saturation variable and the acceptable interpolation error. The metric tensor controls the size of the elements in the mesh. The powerful mesh generator is developed to adapt the mesh based on the computed metric tensor. In the mesh adaptation process, some operations such as the node insertion and deletion, smoothing and edge swapping will be performed. The adaptation process will be repeated until the threshold value is
718
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 6. The algorithm of Bossen Method.
met. Then the pressure and saturation of the previous time step should be interpolated for the new mesh as initial conditions of the new time step. After interpolation the fluid flow equation should be solved for the new time step according to the new adapted mesh. This process is repeated for all the time steps during the reservoir simulation. 3. Initial mesh generation based on well effect The near-well velocity depends on the inertia forces and has a high gradient in this region. Therefore smaller elements are required near the wells to capture the high gradient velocity. For this purpose, we define Element Size Indicator (ESI) concept which is maximum in production and injection wells and then smoothly
vanishes with increasing the distance from the wells. The mathematical function that provides ESI can be shown as follows.
ESI ¼ expðbdÞ
(1)
where, ESI is Element Size Indicator, d is the distance from the well and b is the constant that controls damping intensity from the well. Fig. 2 shows ESI based on the well effects for a reservoir with two wells; injection well in bottom-left and production well in upperright of the reservoir. The ESI is converted into the element size function by assigning the desired minimum element size to maximum ESI and the maximum element size to minimum ESI (Fig. 3). In order to build the initial mesh, the resulted element size function and the reservoir boundary are imported to the mesh
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
719
Fig. 8. A connection in data structure.
Fig. 7. Bossen's smoothing function.
generator which will be described in the section 5. Fig. 4 shows the resulted initial mesh for ESI map in Fig. 2. 4. Fluid flow modeling Control Volume Finite Approximation (CVFA) method is used to discretize the governing equations in porous media (Li et al., 2003). The fully implicit approach is considered to discretize equations in time. The combination of Darcy's law and the individual phase conservation law are employed to describe the flow of wetting (w) and non-wetting (nw) phases in the porous media:
vðfraS Sa =Ba Þ r ¼ V: aS ua þ qa vt Ba kra ua ¼ kVpa
ma
a ¼ w; nw
a ¼ w; nw
Fig. 9. A node in data structure.
(2)
(3)
Sw þ Snw ¼ 1
(4)
Pc ¼ pnw pw
(5)
where, f is the porosity of the porous media, Pc is the capillary pressure, Sa is the saturation of each phase, Ba is the formation volume factor of each phase and raS is the density of each phase at the standard condition, k is the absolute permeability tensor; kra, pa and ma are the relative permeability, pressure and viscosity of each phase. For each node in the discretized domain, the initial pressure and the water saturation are defined. For the outer boundary (reservoir boundary), no-flow boundary condition and for the inner boundary, constant bottom-hole pressure and constant injection flow rate are considered.
all the element sizes will be adapted to the element size function (Fig. 5). The input functions to the grid generator are the boundary of the computational domain and an element size function which is defined on the fine-scale cartesian mesh. First Delaunay Triangulation is implemented in the computational domain. Then, in each step of the iterative process, a node is selected randomly and its position is smoothed according to its neighbors. Finally, Delaunay Triangulation condition is applied and then by considering the nodes density, a node will be inserted or deleted. This process will be repeated until final mesh achieved. Fig. 6 shows the mentioned algorithm details that will be elaborated in the followings.
5. Mesh generation algorithm In this study, Bossen method (Bossen, 1996) in accompany with some modifications which are applied to this method is used to create the dynamic mesh. In this section after reviewing Bossen method, the proposed modifications to decrease the gridding time will be described. 5.1. Review of Bossen method In this method, first an initial triangulation is formed then retriangulation, fining, coarsening and smoothing processes in an iterative loop lead the initial triangulation to the final case so that
Fig. 10. Initial polygon and relation between connections.
720
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
function has been suggested as below.
4 f ðdÞ ¼ 1 d4 :ed
(9)
This smoothing function exerts the attractive and the repulsive forces on the nodes same as the intermolecular forces. If the normalized length is less than one then the nodes will attract each other whereas if the normalized length is greater than one then the nodes will repel each other. In addition, unlike the intermolecular forces, the forces values among the nodes in which have the distances close to zero will have limited reduction values. Fig. 7 shows the applied forces between the two nodes versus the normalized length. Also the neighboring points of node i, were defined as the points had connections with point i and the smoothing constant is 0.2. Due to the triangular structure maintenance during the mesh generation process, the new position of point i will only be acceptable when the triangular 0 structure is n't violated. If the segment between pi and pi does not 0 cross any existing edge then pi node position will change into pi . Otherwise, if the segment crosses a border connection then the node will be transferred to the junction location of the border connection. If none of the cases happen, then the node is n't moved from its location.
Fig. 11. Gradual triangulation and connection updating.
5.1.2. Inserting or deleting node As shown in Fig. 6, if in the computational mesh generation, the mesh is sparse or dense in a location then some nodes will be inserted in or deleted from the mesh. In order to express a quantitative measurement of the sparseness for each node or connection a quantitative parameter which is called Extent is defined. The extent has a normalized area unit. The normalized area parameter is similar to the normalized length in the surface and it can be calculated from the area divided by the desired area (r2) which is computed from the element size function. The following definitions have been presented for the extent of the nodes and the connections.
Fig. 12. Completion of polygon triangulation and connection relation statues.
5.1.1. Smoothing In order to generate an isotropic mesh by considering a specified element size function, the nodes distances can be determined using the following equation. In this equation, k k is a norm operator, and r is the desired distance which can be accessible using an element size function.
dðp1 ; p2 Þ ¼
kp1 p2 k r
(a) Border Connection Extent: half the square of its normalized length (b) Border Node Extent: half the square of the two neighboring connections normalized lengths sum (c) Inner Connection Extent: 0.8 times two adjacent triangular normalized areas pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (d) Inner Node Extent: 2:3= nðn 2Þ times the n surrounding triangular normalized areas sum In the algorithm which is shown in Fig. 6, if i is the selected node
(6)
The proposed equation can be called Normalized Length. Based on this definition, the normalized length of the desired mesh is 1 and the smoothing can be determined using the following first order motion equation. 0
pi ¼ pi þ ai
X f d pi ; pj :uij ;
(7)
j2Ni
pi pj uij ¼ d pi pj
(8)
In the Eq. (7) and Eq. (8), pi is the spatial position of node i, uij is a unit normal vector, Ni is the set of node i adjacent nodes, f is a smoothing function, ai is a smoothing constant. Bossen smoothing
Fig. 13. Initial guess in one-dimensional gridding.
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
721
assumed active. After smoothing pi node and by considering the 0 displacement value of pi pi , the node will be considered as active or inactive. If this displacement value exceeds a threshold value, then the node will be remained active and adjacent nodes will be assumed active. Otherwise it is considered as inactive. In case of the node insertion or deletion, its surrounding nodes will be active. If all nodes be inactive, the algorithm process will be stopped. 5.2. Implementation of Bossen method
Fig. 14. An Element size function and a line segment.
in the algorithm, its extent will be calculated and if its extent is less than 1 then this node will be deleted and the mesh will be locally re-triangulated. If node i is not deleted then, its surrounding connection extent will be calculated. If the extents largest value is greater than 1, its relevant connection will be split through the middle and the node will be added to the mesh and others remained re-triangulated using Delaunay method. The initial border points are never deleted or repositioned. It should be noticed; re-triangulation processes in case of the node insertion will be performed using Delaunay gradual triangulation method. In case of node deletion, the simple achieved polygon is triangulated and by swapping the connections in the polygon, Delaunay triangulation condition will be applied.
5.1.3. Convergence The entire process of smoothing, triangulation, deletion and insertion node, lead all the normalized lengths of the connections to 1. In order to assess the mesh convergence statues, entire node classified as active or inactive nodes. At first all the nodes are
The object-oriented programming and Cþþ language program are used for mesh generation process. A new data structure is used to implement this method. In this data structure, connections are considered as the main components. These components control the data relationship. In each connection two pointers to two nodes at both ends of the connection and two pointers to two adjacent connections exist with assuming counter-clockwise rotation around the mentioned nodes. In order to complete the data structure and retrieving the nodes location, at each node only a pointer to a connection is considered. Accordingly, there is no single node in data structure and each node is connected to a connection (Figs. 8 and 9). Based on the mentioned data structure, in the first step, connections in the polygon which are defining the boundaries of computational domain are established. Then an initial triangulation within the original polygon is formed by simple triangulation method. The important item at this stage is the updating of the data structure gradually along with the completion of the initial triangulation. After the initial triangulation, the data structure between nodes and connections will be completed. The data structure will be updated by adding or removing nodes. Fig. 10 shows the input polygon and the relation between connections at the initial state. Fig. 11 shows the gradual triangulation and the updating data structure and Fig. 12 shows the completion of the initial triangulation and the relationship between connections. It should be noticed when removing a node, a simple polygon is formed and so, the data structure updating is a similar process to the initial start status. Advancing Front method and Delaunay condition are also used to Delaunay triangulation in a simple polygon. 5.3. Modification on Bossen method In order to improve the efficiency of Bossen method, two
Fig. 15. Initial one-dimensional gridding and final gridding after smoothing.
722
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 16. First stage algorithm in two-stage gridding.
modifications are applied to the mentioned method. The first is the placement of the nodes on the edges of the polygon (boundary nodes) based on the element size function and separation of this stage as one-dimensional gridding. The second modification is the implementation of the original algorithm (Fig. 6) in two stages to speed up the grid generation process. In the first stage an estimation of final mesh will be obtained quickly by limited smoothing and in the second stage final smoothing will be performed. In the two following sections these variations are described.
5.3.1. Boundary gridding separation of inner gridding After importing computational domain boundaries as a polygon, each edge is considered as one-dimensional space and the position of the boundary nodes will be fixed based on the element size function on the edges. As aforementioned, in the proposed algorithm in Figs. 10e12, a triangle structure will be preserved during mesh generation. By applying mentioned modification, we will only check the condition of the existence of the node position insides it's surrounded polygon in the smoothing process and it is not
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
723
Fig. 17. Second stage algorithm in two-stage gridding.
0
required to find a possible collision line segment between pi and pi with polygon boundaries surrounding it (pi is node location before 0 smoothing and pi is node location after smoothing). Also, it doesn't have the necessity to check the crossed edge whether is a boundary edge or not. In order to perform one-dimensional gridding which leads to the determination of the boundary nodes, the node position of the two ends of the segment will be fixed and then an initial guess of the intermediate nodes position estimated. In the next step, in an iterative process including node removing and smoothing, the initial guess will be optimized based on the element
size function. To estimate the initial guess, a simple algorithm is used. First the desired lengths based on the element size function in the boundary nodes are computed. If the summation of these lengths half is less than the segment length of the two boundary nodes, then a node will be added in the middle of the free zone. Free zone is defined as the region between the segment of the two boundary nodes and out of the desired lengths half. Fig. 13 shows a schematic of a step in estimating the initial guess. In this figure d1 and d2 are desired lengths based on the element size function for both ends of the
724
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 18. Boundary of computational space and initial guess of mesh.
Fig. 21. Permeability map of a heterogeneous reservoir. Fig. 19. Final mesh based on the element size function.
with the largest absolute deviation from the normal lengths will be identified and removed. Then smoothing is performed again on the nodes and the total absolute deviation from the normal lengths of the nodes will be calculated. If the absolute deviation is decreased, the last position of the nodes will be stored. This process continues to minimize the absolute deviation from the normal length. Fig. 14 shows an element size function and a line segment in the
Table 1 Water relative permeability in oilewater system.
Fig. 20. An ellipse for metric tensor definition.
segment. Above method is repeated in an iterative process for the new segments until no free zone between the nodes exists. After performing the initial guess, smoothing conducted on the intermediate nodes simultaneously. The position of the nodes will be stored after the initial smoothing. The total absolute deviation from the normal lengths of the nodes is calculated and the nodes
Sw
krw
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.8 0.9 1
0 0.0028 0.025 0.07 0.14 0.22 0.33 0.62 0.80 1
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729 Table 2 Oil relative permeability.
725
Table 3 Simulation running times for three models in oilewater flow example.
So
kro
Model
Run time (second)
0 0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9 0.95
0 0.011 0.044 0.177 0.277 0.399 0.543 0.709 0.897 1
Fine-grid Uniform isotropic static grid Non-uniform dynamic grid
21,282 380 378
computational space. Fig. 15 shows the initial guess grid and final grid using the mentioned algorithm. 5.3.2. Two-stage gridding Boundary gridding separation of the inner gridding in the computational domain provides a more coherent smoothing process. The smoothing in terms of the computation time is the most costly step in the algorithm shown in Fig. 6. Therefore, we decided to conduct the smoothing process in the two stages. First it is attempted to provide a primary estimation of the nodes position in the mesh and then primary mesh will be lead to the final mesh. Smoothing of each node position based on the first order motion equation is performed in several steps. The smoothing process is as follows. f
pi /smoothing/p1i /smoothing/p2i …/pi
(10)
where, pi is the first position of node i. This position is transferred to the next position during the smoothing process and finally converges to pfi position. If the distance between the two nodes in the two consecutive steps is defined as Dn pi ¼ pnþ1 pni , then by i using a suitable smoothing function, it is expected that the following relationship will be established.
D0 > D1 > … > Df 1
after the first step in the smoothing process is considered as the final position (pfi ¼ p1i ). Using the method in the algorithm which is shown in Fig. 6, an estimation of the final mesh is achieved quickly. After determining the initial position of the nodes, the final smoothing process will be conducted. Now, all the nodes except the boundary nodes are considered as active nodes and smoothing is conducted on all the active nodes simultaneously. After smoothing operation, Delaunay Triangulation condition is applied to the grid connections. Then the extent of all the active nodes will be examined. According to the computed extents, nodes may become inactive, be removed or one node be inserted. These steps are continued until all nodes in the mesh will become inactive. The steps to remove, deactivate, activate and insertion of the new node in the mesh are conducted according to the algorithm which are shown in Figs. 16 and 17. These figures show the new algorithm based on the two-step gridding. An example of the element size map, boundary of the computational domain and the initial estimation of final mesh are depicted in Fig. 18. Fig. 19 shows the final mesh in accordance with the element size function.
5.4. Element size function computation The Metric Tensor which controls size, shape and orientation of the elements is required for the gridding. The metric tensor will be defined as an ellipse with the two radiuses r1, r2 and the angle q (Fig. 20). The corresponding metric tensor is (Bossen, 1996):
(11)
According to a threshold value the smoothing process is stopped and pfi is determined. In the initial estimation of the node position, a limited smoothing is performed. It means that the node position
M¼ ¼Q
cos q sin q
g1 0
sin q cos q 0 QT
g2
g1 0
0
g2
cos q sin q
sin q cos q
(12)
where, g1 ¼ 1=r12 and g2 ¼ 1=r22 are eigenvalues of metric tensor and Q is the rotation matrix of angle q. The following Metric Tensor in two-dimension is computed as (Amirsardari et al., 2015; Hecht, 1998):
M¼
1 jHðuÞj ε0 supðuÞ infðuÞ
(13)
where, ε0, is user defined interpolation error and H(u) is Hessian matrix of water saturation function:
2
v2 Sw
6 6 vx2 HðuÞ ¼ 6 6 2 4 v Sw vyvx
Fig. 22. Water Saturation distribution and dynamic mesh at 40th day.
3 v2 Sw 7 vxvy 7 7 7 v2 Sw 5
(14)
vy2
In order to control the minimum and maximum size of the elements in the mesh, these parameters should be defined (CastroDiaz et al., 1997). For this purpose, gi(i ¼ 1,2) of the metric in Eq. (13) is limited as follows:
726
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 23. Water-cut for fine-scale, uniform coarse-scale and proposed dynamic mesh models.
Fig. 24. Water saturation distribution at 57th day for fine-scale model (right) and uniform coarse-scale static mesh model (left).
6. Results and discussions
1
g~ 1;2 ¼ min max
g1;2
; 2
hmax
! ;
1 h2min
! (15)
where, hmin and hmax are the minimum and maximum edge lengths. Due to the usage of isotropic elements in the study, the largest eigenvalue is considered for both diagonal eigenvalues. So, the isotropic metric is defined as below:
M¼
~2Þ ~1; g maxðg 0
0 ~2 Þ ~1; g maxðg
(16)
In order to validate the proposed method in the study, the twophase flow simulation in the oil and gas reservoirs was performed for the three models: the fine-grid geological model, the uniform coarse-grid model with static gridding and the non-uniform adaptive coarse-grid model with dynamic gridding. The simulation results of the models are used to determine the efficiency and the accuracy of the proposed method in the study. For this purpose, the fine-grid model was considered as a reference model and the simulation running time of the uniform static gridding model is in the same level of the dynamic gridding approach. We present two examples of oilewater and gasewater flows in the oil and gas reservoirs, respectively. The permeability map is shown in Fig. 21. In
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
727
Fig. 25. Water saturation distribution at 57th day for fine-scale model (right) and non-uniform coarse-scale dynamic mesh model (left).
Table 4 Water relative permeability in gasewater system. Sw
krw
0.22 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.01 0.03 0.07 0.15 0.24 0.33 0.65 1
Table 5 Gas relative permeability. Sg
krg
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.78 1
0.04 0.1 0.2 0.3 0.399 0.42 0.62 0.81 1
both examples the reservoir porosity is 0.1 and at the initial state the water saturation is 0. 6.1. Example 1: oilewater flow in an oil reservoir We consider oilewater flow simulation in this example. The initial reservoir pressure is 9014.7 psi. The water and oil relative permeability are shown in Table 1 and Table 2, respectively. A water injection well (Inj0) and a production well (Prod0) considered in the model in the bottom-left and upper-right locations, respectively (Fig. 21). The production well produces with the constant bottom-hole pressure of 5000 psia and in the injection well water is injected with the constant water rate of 450 BBL/D.
Fig. 26. Water saturation distribution and dynamic mesh at 28th day in gasewater flow.
Table 6 Simulation running times for three models in gasewater flow example. Model
Run time (second)
Fine-grid Uniform isotropic static grid Non-uniform dynamic grid
20,878 350 338
The proposed dynamic model approach is used to simulate the oilewater flow in the reservoir. The maximum number of the elements in the model is 2528. The generated adapted mesh at 40th day is shown in Fig. 22. As it can be observed in this figure, smaller elements in the regions with high saturation gradients and near wells were generated using the proposed algorithm. To validate the proposed method, the simulation results were compared to the fine-grid model and the uniform coarse-grid model (static mesh for both models) with 93,690 and 5894 elements, respectively. The number of the elements in the uniform
728
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
Fig. 27. Water production rate for fine-scale, uniform coarse-scale and proposed dynamic mesh models at Prod 1.
Fig. 28. Water production rate for fine-scale, uniform coarse-scale and proposed dynamic mesh models at Prod 2.
coarse-grid model is determined in such a way that the simulation running times in both non-uniform coarse-grid dynamic mesh model and uniform coarse-grid static mesh model be at the same level (Table 3). Table 3 shows that the simulation running times for both coarse-scale models (static and dynamic mesh models) are in portion of 1/55 in comparison to the fine-grid geological model. Water cut (water rate to liquid rate) at the production well during the simulation time in three mentioned models has been shown in Fig. 23. Fig. 24 shows water distribution at 57th day (water breakthrough time in the reference model) for the uniform coarse scale static gridding model and fine-scale static mesh model (reference model). Also, water distribution at 57th day for the nonuniform coarse scale dynamic mesh model and fine-scale static
mesh model (reference model) have been shown in Fig. 25. As can be observed in Figs. 23e25, using dynamic mesh model the water cut and water breakthrough time can be obtained with the high accuracy in comparison to the uniform static mesh model (about 20% improvement in the water breakthrough time). 6.2. Example 2: gasewater flow in a gas reservoir In this example we simulate gasewater flow in a gas reservoir. The initial reservoir pressure is 5300 psi. The water and gas relative permeability are shown in Table 4 and Table 5, respectively. A water injection well (Inj1) and two production wells (Prod1 and Prod2) considered in the model (Fig. 21). The production wells
M. Amirsardari et al. / Journal of Natural Gas Science and Engineering 31 (2016) 715e729
produce with the constant bottom-hole pressure of 3000 psia and in the injection well water is injected with the constant water rate of 800 BBL/D. The proposed dynamic model approach is used to simulate gasewater flow in the reservoir. The maximum number of the elements in the model is 2379. The generated adapted mesh at 28th day is shown in Fig. 26. As it can be observed in this figure, the smaller elements in the regions with the high saturation gradients and near wells were generated using the proposed algorithm. To validate the proposed method, the simulation results were compared to the fine-grid model and the uniform coarse-grid model (static mesh for both models) with 93,690 and 4848 elements, respectively. The number of the elements in the uniform coarse-grid model is determined in such a way that the simulation running times in both non-uniform coarse-grid dynamic mesh model and uniform coarse-grid static mesh model be at the same level (Table 6). The water rates for two production wells in three mentioned models have been shown in Fig. 27 and Fig. 28. As can be observed in these figures, using the dynamic mesh approach the water breakthrough time and the water production rate can be computed with the high accuracy. 7. Conclusions and recommendations A new automatic reservoir simulation approach based on the dynamic mesh concept has been introduced. The advantages of this methodology and the recommendations for the future works have been summarized in the following items: 1. One of the main issues in the use of the dynamic mesh is the simulation running time. In this study, the non-uniform adaptive coarse gridding is used to have the smaller elements only in the dynamic interest regions in order to optimize the simulation running time. 2. A flexible two-dimensional computation mesh generator is used for the dynamic gridding. In order to improve the efficiency and decrease the gridding time, some modifications have been applied to this mesh generator algorithm. 3. The creation of dynamic mesh is based on the saturation variations and the well effects. Any other parameters (static or dynamic) can be used to build the dynamic mesh for the future works. 4. The numerical tests of the presented methodology in the oil and gas reservoirs show that the use of the mentioned approach makes it possible to obtain the high accurate results with the low computation cost. 5. As all the mentioned steps in the proposed algorithm can be extended to 3D models, this method will be applicable in case of 3D modeling. Nomenclature M f pa S
a k kra
raS Ba
Metric tensor Smoothing function Pressure of phase a Saturation Phase (oil/water/gas) Permeability tensor Relative permeability of phase a Surface density of phase a Formation volume factor of phase a
Pc ua t
ma p H ESI ε0 d
f g hmin hmax sup inf
729
Capillary pressure Velocity of phase a Time Viscosity of phase a Spatial point Hessian matrix Element Size Indicator Interpolation error Distance Reservoir porosity Eigenvalue of metric tensor Minimum edge length Maximum edge length Supremum Infimum
References Aarnes, J.E., Lie, K.A., Kippe, V., Krogstad, S., 2009. Multiscale Modeling and Simulation in Science, vol. 66. Springer Verlag, p. 27. Amirsardari, M., Dabir, B., Naderifar, A., 2015. Anisotropic grid generation for dynamic simulation in heterogeneous hydrocarbon reservoirs. J. Nat. Gas. Sci. Eng. 27, 1523e1535. Bossen, F.J., 1996. Anisotropic Mesh Generation with Particles. Msc. Thesis. School of Computer Science Carnegie Mellon University. Cao, J., Kitanidis, P.K., 1999. Adaptive-grid simulation of groundwater flow in heterogeneous aquifers. Adv. Water Res. 22, 681e696. Castro-Diaz, M.J., Hecht, F., Mohammadi, B., Pironneau, O., 1997. Anisotropic unstructured mesh adaption for flow simulations. Int. J. Numer. Methods Fluids 25, 475e491. Chen, T., 2009. New Methods for Accurate Upscaling with Full-tensor Effects. Ph.D. Thesis. Stanford University, Stanford, CA. Cirpka, O.A., Frind, E.O., Helming, R., 1999. Streamline-oriented grid generation for transport modelling in two-dimensional domains including wells. Adv. Water Res. 22, 697e710. Ebrahimi, F., Sahimi, M., 2004. Multiresolution wavelet scale up of unstable miscible displacements in flow through heterogeneous porous media. Transp. Porous Media 57, 75e102. Evazi, M., Mahani, H., 2010. Unstructured-coarse-grid generation using background-grid approach. SPE J. 15, 326e340. Farmer, C., 2002. Upscaling: a review. Int. J. Numer. Methods Eng. 40, 63e78. Garcia, M., Journel, A.G., Aziz, K., 1992. Automatic grid generation for modeling reservoir heterogeneities. SPE Reserv. Eng. 7, 278e284. Heinemann, Z.E., Brand, C.V., 1989. Gridding techniques in reservoir simulation. In: Proceedings of the First and Second International Forum on Reservoir Simulation, Alpbach, Austria. Heinemann, Z.E., Brand, C.W., Munka, M., 1991. Modeling reservoir geometry with irregular grids. SPE Reserv. Eng. 6, 225e232. Hecht, F., 1998. Bidimensional Anisotropic Mesh Generator. Technical Report, INRIA, Rocquencourt. King, M.J., Burn, K.S., Wang, P., Muralidharan, V., Alvarado, F.E., Ma, X., DattaGupta, A., 2005. Optimal coarsening of 3D reservoir models for flow simulation. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Li, B., Chen, Z., Huan, G., 2003. Control volume function approximation methods and their applications to modeling porous media flow. Adv. Water Res. 26, 435e444. Mlacnik, M.J., Harrer, A.W., Heinemann, Z.E., 2003. Locally streamline-pressurepotential-based PEBI grids. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Palagi, C.L., 1992. Generation and Application of Voronoi Grids to Model Flow in Heterogeneous Reservoirs. Ph.D. Thesis. Stanford University, Stanford, CA. Portella, R., Hewett, T.A., 1998. Upscaling, gridding, and simulation using streamtubes. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Tran, T.T., 1995. Stochastic Simulation of Permeability Fields and Their Scale-up for Flow Modeling. Ph.D. thesis. Stanford University, Stanford, CA. Verma, S., Aziz, K., 1997. A control volume scheme for flexible grids in reservoir simulation. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. mez-Herna ndez, J.J., 1998. Upscaling hydraulic conductivities in Wen, X.-H., Go cross-bedded formations. Math. Geol. 30, 181e211. Wen, X., Durlofsky, L., Edwards, M., 2003. Upscaling of channel systems in two dimensions using flow-based grids. Transp. Porous Media 51, 343e366.