Meshing of domains with complex internal geometries - CiteSeerX
NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS Numer. Linear Algebra Appl. 2006; 00:1–14
Prepared using nlaauth.cls [Version: 2002/09/18 v1.02]
Meshing of domains with complex internal geometries Randi Holm∗ , Roland Kaufmann, Bjørn-Ove Heimsund, Erlend Øian and Magne S. Espedal Center for Integrated Petroleum Research Allgt. 41, N-5007 Bergen, Norway
SUMMARY This paper presents a meshing algorithm for domains with internal boundaries. It is an extension of the gridding algorithm presented by Persson and Strang. The resulting triangulation matches all boundaries, and the triangles are all nearly equilateral. Equilateral triangles are beneficial for a finite volume discretization, as fluid flow between elements of very different size is only possible at small timesteps. The mesh generator is compared with the well regarded Triangle program, where both element quality and simulation performance are checked. It is shown that c 2006 John Wiley & Sons, Ltd. our mesh generator consistently delivers better meshes. Copyright KEY WORDS :
physically-based mesh generation; fractured porous media; subsurface flow simulation
1. INTRODUCTION Flexible meshes are necessary in order to conduct detailed numerical studies of fluid flow and transport in geometrically complex geological systems. Our goal is to generate computational meshes of high quality for fracture systems and layered porous media. We focus on two mesh quality criteria, geometry adaptation and element shape. For general problems these criteria are in conflict, i.e. good geometrical fitting can lead to poor element quality. Recently, a mesh generation algorithm was presented by Persson and Strang [14]. This algorithm iteratively improves an initial mesh by assigning repulsive forces between the nodes on either side of an edge, and then moving the nodes such that global force-equilibrium is obtained. At each step of the algorithm, a new triangulation is obtained by using a Delaunay triangulation method [3, 5]. We present an extension to this algorithm such that the mesh will conform to internal boundaries such as material discontinuities and fractures, in addition to the external boundaries of the domain. The idea of the extension is to restrict the repulsive forces acting on boundary nodes such that the nodes may only move along the boundary, and not leave it. This simple modification uses a constrained Delaunay triangulation method [18] to ensure that no edges of the mesh cross an internal boundary. In the next section we give a geology oriented motivation for this work on unstructured mesh generation. The fundamentals of the meshing algorithm described by Persson [14] are then given in
∗ Correspondence to: [email protected]
c 2006 John Wiley & Sons, Ltd. Copyright
2
R. HOLM ET AL.
Figure 1. The figure shows a domain that contains a stochastically generated set of fractures.
Section 3. That section continues with our extension of the formulation that enables adaptation of the mesh to internal boundaries. This gridding method is then applied to some representative and realistic fracture networks in Section 4, and a comparison to a popular meshing algorithm [16] is made. The comparisons are rounded off by performing some simple water/oil fluid flow simulations.
2. GEOLOGICAL MOTIVATION Geological systems often exhibit complex geometrical structures that complicate the mesh generation. Examples of such include faults with associated fractured damage zones, naturally fractured rocks and layer truncations and pinch-outs in porous media. In particular, faults and fractures in the subsurface may act as barriers or conduits for fluid flow. Consequently, such structures are important factors for determining the flow paths of fluids in porous and non-porous rock masses. In order to get reliable flow simulations, it is then vital to develop robust meshing algorithms that accurately represent the geometrical properties of these features. Important steps in a work flow from problem definition to flow simulation are data acquisition, characterization, petrophysical modeling, meshing and simulation. In this paper we only consider the task of generating a well adapted mesh from some digital characterization of the geometrical structures, and consider its impact on simulations. For concreteness, we motivate the mesh generation task from challenges of meshing faulted and fractured porous rocks. In particular, an ongoing work within our research group is concerned with fluid flow near and in fault zones [20]. Note that by fractures we understand both high permeable open fractures and low permeable deformation bands (“anti-fractures”). The major task is to create a high-quality mesh that is well adapted to these fractures. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
3
Figure 2. A photograph of a meter scale outcrop and the corresponding bitmap of (manually) interpreted layers and fractures/deformation bands.
2.1. Characterization Characterizations of fractured porous media can be categorized into two main classes, a stochastic and a deterministic approach. It is also possible to combine the two into a hybrid stochastic-deterministic approach. The latter is common to apply when particular faults or fractures are known. In the stochastic approach, synthetic fracture systems are generated based on statistical parameters, see Figure 1. The generating algorithms can either be based on purely statistical distributions or on more elaborate schemes that try to mimic geometrical structures of nature [13]. In the stochastic approach, all fractures typically have an analytical geometrical representation. Consequently, these can be used directly as input in the mesh generation procedure. The deterministic approach includes using e.g. outcrop photographs or seismics combined with some image processing technique to recognize the fractures [6, 8]. In Figure 2 below, this is illustrated for a real outcrop example, where the “image processing” technique used is simply manual interpretation by geologists within an image manipulation program. 2.2. Meshing challenges In both characterization approaches described above, there are some common features of most fracture realizations. Reichenberger [15] describes four scenarios which are difficult to handle for mesh generators; see Figure 3. The flow connectivity is affected by slightly crossing, Figure 3(a), and almost crossing fractures, Figure 3(b), and can be important to resolve accurately. It is also challenging to generate high quality elements that adapt to close, (almost) parallel fractures, Figure 3(c), and fractures crossing at a small angle, Figure 3(d). They can all be identified from the example characterizations in Figure 1 and Figure 2 above. On the other hand, the fluid flow simulation performance depends on well-shaped elements. In a control volume formulation, fluxes between adjacent cells are calculated from the cell pressures. For a homogeneous medium consisting of equilateral triangular elements, the flux can be calculated using just the two neighboring cell values [1, 2], which gives a two-point flux stencil. Deviations from an equilateral form implies that multiple cells are coupled for the calculation of fluid fluxes, which is a multi-point flux stencil. In general, multi-point flux stencils lead to more work for flow simulators [7]. The focus of this paper is to present additions to a mesh generation algorithm such that domains with complex internal geometries can be treated and be amenable to flow simulations. Specifically, we seek an algorithm that honors the geometries of faulted and fractured rocks to yield well-formed elements. In this case, fracture end-points and fracture intersections must be accurately represented. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
4
R. HOLM ET AL.
(a)
(b)
(c)
(d)
Figure 3. The figure shows different classes of fracture configurations that are difficult to mesh: (a) slightly crossing fractures; (b) almost crossing fractures; (c) close and (almost) parallel fractures; and, (d) fractures crossing at a small angle.
These are particularly important in order to capture the effects of fluids accumulating in regions limited by low-permeable fractures (compartmentalization) and fluid distributions caused by the connectivity of high-permeable fractures.
3. FORCE EQUILIBRIUM MESHING We will now briefly review the meshing algorithm given by Persson and Strang [14], and then extend it to handle internal interfaces, such as layering or fractures. The algorithm is iterative, where in each step the nodes are moved towards a force equilibrium position and a new topology is formed by a Delaunay triangulation method. The initial condition is an initial point distribution connected by a triangular mesh. The discussion starts with the algorithm for moving the nodes. This is followed by methods for adapting the nodes to internal boundaries along with suitable Delaunay triangulation schemes. 3.1. The edge forces To each edge we assign a force function f (`, `0 ), whose strength depends on the current edge length ` and the target length `0 . As argued by Persson [14], we may set ( `0 − ` ` < `0 , f (`, `0 ) = 0, ` ≥ `0 . Positive forces are repulsive, and hence only edges shorter than their target lengths are assigned a force. Attractive forces could be included by just setting f = `0 − `, but this may reduce the mesh mesh quality. The target length `0 is just a relative quantity that can be used for mesh adaptivity. Its absolute value carries little significance, since the mesh density is determined by the number of nodes in the domain. Consequently, a varying `0 serves to cluster nodes (mesh adaptivity), but will not determine edge lengths. 3.2. The node forces It is the nodes of the mesh which will be moved, and not the edges. Nodes are moved using the summed edge forces of all the edges connected to that node. After moving the nodes, a new set of edges are c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
5
found from a Delaunay triangulation algorithm. The repulsive forces acting on node i is X f ` j , ` j,0 pEi − pEi, j , FEi = j
where the summation is over all the edges connected to node i. pEi is the node coordinate, and pEi, j is the opposing node coordinate along edge j. Since f ≥ 0, the summation into FEi is positively proportional to the vectors pointing at node i. Hence the node forces are positive and repulsive, just as the edge forces, and therefore the nodes will tend to move away from each other. 3.3. An explicit algorithm After assembling all the nodal forces, it remains to solve the coupled system FEi ( pE1 , pE2 , . . . , pEm ) = 0,
for i = 1, . . . , m,
where m is the number of nodes in the mesh. The node forces are dependent on the topology, which in turn has been generated by a Delaunay triangulation. It is consequently hard to solve the equilibrium equations directly since FE is a discontinuous function of the points. Instead, a pseudo time stepping approach is adopted, whereby ∂ pEi = FEi , ∂t
for i = 1, . . . , m,
which is discretized explicitly as pEin+1 = pEin + 1t FEin ,
for i = 1, . . . , m.
(1)
After each time step, new edges are generated from pEin+1 by a Delaunay triangulation algorithm. The iteration terminates once k FEi k ≈ 0, ∀i. 3.4. Delaunay triangulation Perhaps the most common way of triangulating the mesh is the Delaunay scheme. For a Delaunay scheme, triangles are chosen such that no nodes are within the circumcircle, i.e. the circle that goes through all three nodes of the triangle [5]. This is done by picking the nearest neighbor to a node as the one with which it will form an element [4]. An illustration of the concept is given in Figure 4. This small example contains a domain of four nodes. Figure 4(a) shows the circumcircles for each of the four different combinations of three nodes, whereas Figure 4(b) displays the final triangles. As point C is within the circumcircle of AB D, this disqualifies the points A, B and D from forming a triangle. Instead, 4ABC is grown. It has been showed that the Delaunay triangulation gives the largest minimum angle in any element, which creates as equilateral angles as possible [11]. Note however, that it does not necessarily give the largest sum of minimum angles overall. 3.5. Adjusting the nodes to boundaries We now turn to the major topic of the paper: aligning the mesh to boundaries. In our case, there are two types of boundaries: external and internal. External boundaries enclose the domain, while the internal c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
6
R. HOLM ET AL.
A
A
A C
C
C
B
B
B D
D
D
(b) Regular Delaunay
(a) Circumcircles
(c) Constrained Delaunay
Figure 4. Triangulation of four points.
boundaries are the fractures which are to be honored. Points which are pushed outside the domain are simply projected back onto the boundary. For a boxed domain, this projection simply clamps the coordinate to lie within the box. Internal boundaries (fractures) are handled differently. Using a target density, points are initially distributed evenly throughout the domain. The same density is applied along each fracture, and points off the fractures but deemed too close to them are removed to obtain the desired density. Crossing fractures are split at the intersections, which yields a set of non-intersecting fractures. The nodes at the endpoints of the fractures are locked to their position, and the nodal forces acting on them are always set to zero. A point initially placed on a fracture will always be projected back onto that fracture after Equation (1) has been solved. Referring to Figure 5, let the vector uE be the fracture and vE be the vector that connects one endpoint of the fracture to the point to be projected back. From standard vector calculus, the desired vector vEu is vEu =
vE · uE uE. uE · uE
This uses the property that vE · uE = kE v kkE u k cos θ . 3.6. Constrained Delaunay methods While the nodes will now properly obey all geometrical constraints, the resulting Delaunay triangulation can generate edges which cross fractures. To avoid this, we resort to a constrained triangulation where edges connecting nodes along fractures are kept, and the triangles which are generated cannot have edges crossing the fractures. A constrained Delaunay scheme relaxes the demands on the triangles. The triangles may now have other points inside their circumcircle if those points are on the other side of a constraining barrier. Figure 4(c) depicts how the triangulation may change after introducing constraints. Let us assume that it is desirable to have the line AD as an edge because this represents a fracture. Now, the point C is occluded behind this subsegment, and 4AB D is considered legal. Care must be taken to start the algorithm from a boundary, so that a situation where a triangle crosses the boundary does not arise. Triangles are otherwise grown similarly, and keep many of the desirable properties, such as for instance having the largest minimum angle [12], even if they are no longer strictly Delaunay. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
7
vE C C : uE C C : vEu θ Figure 5. Projecting a point at vE back onto its associated fracture uE, yielding vEu .
There are also some other triangulation methods which we will briefly discuss: Conforming Delaunay combines refinement with constraints by adding points along internal boundaries until the fractures are part of the triangulation while simultaneously maintaining the Delaunay property [17]. Similar is the almost Delaunay procedure, which add points at the intersections between edges and constraining segments, thus ensuring that the edges will be split during the next iteration. Since we prefer to adjust the node geometry using the force equilibrium approach, the constrained Delaunay method fits well into our approach. The conforming and almost Delaunay schemes add points, and may disturb the equilibrium.
4. EXPERIMENTS Based on the preceding discussion, we have implemented an augmented version of Persson’s algorithm with the following enhancements: • Internal boundaries are divided into subsegments with approximately the same size as the expected lengths of interior elements. • Interior points that are close to internal boundaries are removed to avoid crowding. • End-points of internal boundaries remain fixed. • Points that originated at an internal boundary are projected back to that boundary after force application. • A constrained Delaunay triangulation method ensures that no edges cross internal boundaries. 4.1. Quality measures When it comes to using a triangulation in a finite volume method, size does matter. Fluid flow between elements of very different size is only possible at small timesteps, which increase simulation time or even worse, makes it hard to get an accurate result, if any at all. To measure the size quality of an element, we compare the deviation of its size from the mean size of all elements, and adjust this so that the quality ends up between zero (for infinitely different element sizes) and one (for perfectly equal element sizes). All terms are squared to penalize elements further away from the mean. Thus, we end up with the expression: qsize = where ai is the size of element i and a¯ = c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
1 n
a¯ 2
P
i
a¯ 2 , − (a¯ − ai )2
ai is the mean size of all elements. Numer. Linear Algebra Appl. 2006; 00:1–14
8
R. HOLM ET AL.
Another way of measuring size quality would be to use the ratio between sizes of adjacent triangles only. However, due to time constraints this has not been investigated. Since the Delaunay method is used to maximize the minimum angle, element quality has traditionally been measured by looking at the angles. We use the same geometrically intuitive measurement as Persson [14], namely the ratio between the smallest circle that encompasses the triangle (the circumcircle) and the largest circle that fits inside the triangle (the incircle). This ratio is divided by two so that it is normalized between zero (for infinitely thin triangles) and one (for equilateral triangles): (h 1 + h 2 − h 3 )(h 1 + h 3 − h 2 )(h 2 + h 3 − h 1 ) . h1h2h3 An alternative quality measure for triangles that is commonly used is given by Johnson [9]: √ 4 3a qtriangle = 2 , h 1 + h 22 + h 23 qangle =
where a is the area and h 1 , h 2 and h 3 the side lengths of the triangle. If h 1 = h 2 = h 3 , the triangle quality is optimal and q = 1. Quality is usually deemed acceptable if qtriangle > 0.6. We study the properties of our approach by meshing domains of different origin and compare with grids generated by a conventional method. The methods will be compared using the quality measures described above in conjunction with visual inspection. For each of the examples we provide statistics on the mean over all cells, the standard deviation and the first percentile of the quality. The mean, µ, will give us indication of the overall quality of the generated mesh, whereas the standard deviation, σ , tells us to which degree that quality is homogeneous throughout the elements. As we will see, a mesh can have a very good overall quality but still have some regions of lower quality. It will typically be those elements that pose problems during a simulation. The first percentile indicates in which range the elements with lowest quality lie. We could also have used the minimum quality found, but in cases with a significant density of fractures there will typically be at least one element of very low quality. By using the first percentile instead, we are less prone to be biased by those pathological cases in our analysis. 4.2. Benchmark As a benchmark to which we can compare our own approach we use Jonathan Shewchuk’s wellrenowned Triangle program, [16]† , which may be considered the state of the art of triangulation. It is our experience that this program in general generates very good triangulations. However, to maintain quality it always refines the initial mesh by splitting elements into smaller ones, thus having a tendency to end up with “cobwebs” around fractures that crosses at a steep angle. In contrast, our approach attempts to minimize the global tension between points, which means that tricky areas are avoided rather than successively refined. Some of these problems may be rectified by altering triangulation parameters. For instance, relaxing the angle quality constraint may in some cases give an overall better quality. These tweaks tend to be very mesh-specific, as the program already does an excellent job for the general case out of the box with the default settings.
† Version 1.6 released on July 28, 2005
c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
9
We have also created a patch to the Triangle program that inhibits the refinement cascade for elements below a certain specified size‡ . Preferring size over angle quality, this method will yield long, thin “needles”. The reader should be informed that this method is without any theoretical guarantees. 4.3. Deterministic The algorithms were used to create a triangulation of the outcrop in Figure 2 that was interpreted manually. We strived to get triangulations that contained about the same number of elements as possible to get comparable results. Figure 6(a) is made by Triangle. Note the little cobweb that appears around the middle of the upper right quadrant. Usually, these clusters indicate that the lines do not quite match up or make steep angle. We have manually inspected this location and there is nothing special about it; it simply seems to be an unfortunate angle. When inhibiting refinement in Triangle (Patched), we get Figure 6(c). At first sight there appears to be no irregularities, but upon closer inspection a lot of thin triangles can be found. None of these are detrimental, but by having so many of them, they represent an impediment for an effective flow simulation. The best results were obtained using the revised force-equilibrium algorithm (Persson), which is shown in Figure 6(e). Most of the thin triangles are gone; those that are left are typically in places that would be hard to triangulate anyhow. Table I shows how the various methods compare to each other. Not surprisingly, the revised forceequilibrium and the patched refinement algorithms are able to generate higher size quality on average, but they also have a generally better distribution of the sizes as well, in terms of lower standard deviation and higher percentile mark. A result that we had not expected was that the revised force-equilibrium and the patched refinement algorithms also perform better when it comes to angle quality. Whereas the Triangle algorithm may produce the best result when viewing the possibilities for one triangle isolated (given the triangulation that has occurred up till this point), it apparently does not translate into finding a good mesh overall. Graphically, the same trend may be shown in a sensitivity plot. Figure 7 shows the qualities for meshes of various sizes. All algorithms are able to get better qualities for more elements, which is natural given that the distance between the constraints increases. However, the force-equilibrium and the patched refinement maintain a steady lead, for both measures. We are unsure of why the patched refinement program manages to do so well for small meshes, and we would advise against reading too much into that result. Also note that the graphs are not strictly increasing; in some cases the quality may actually decrease a bit even if more elements are made. This shows that there is still some sensitivity for the actual configuration itself. However, we think that the trend should be clear. 4.4. Stochastic To test out the algorithms on a relatively random example, we used the program Frac3D [19], developed at the University of Stuttgart, to generate the set of fractures in Figure 1. The program works by drawing fracture length, orientation and placement from statistical distributions.
‡ Currently we use refinement based on the centroid as an estimator of the child triangles’ sizes, which admittedly may be very
inaccurate for bad triangles. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
10
R. HOLM ET AL.
(a) Triangle, N = 3164
(b) Triangle, N = 5663
(c) Patched, N = 3108
(d) Patched, N = 5583
(e) Persson, N = 3134
(f) Persson, N = 5614
Figure 6. Meshes for the deterministic case (left column) and the stochastic case (right column). The fractures are drawn with a larger linewidth. About N ≈ 3100 triangles were used for the deterministic case, and about N ≈ 5600 triangles in the stochastic case. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
11
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
Table I. Quality statistics for the meshes generated from the deterministic fracture network, with about N ≈ 3100 triangles. All statistics have been scaled by 100 and given as percentages.
triangle patched persson
µ
qangle σ
1%
µ
86.5 86.8 90.2
14.1 13.8 11.3
32.0 31.8 48.6
88.7 88.9 91.6
qtriangle σ 1% 11.4 11.2 9.5
46.3 45.9 55.4
µ
qsize σ
1%
94.0 94.8 94.9
7.9 5.9 7.5
57.9 76.0 64.4
92 96 91 95
90
89
94
88 93 87 92 86 91
85
84
90 persson triangle patched
83
82
1500
2000
2500
3000
3500
(a) Angle quality, qangle
4000
persson triangle patched
89 1500
2000
2500
3000
3500
4000
(b) Size quality, qsize
Figure 7. Mean quality vs. number of triangles for a deterministic outcrop.
In this case Triangle performed poorly. As shown in Figure 6(b), many small elements were required to match the fracture set. Most of these are removed and replaced by single triangles in our modification, as depicted in Figure 6(d). In contrast, this is a case where the force-equilibrium variant really shines. There are a lot fewer bad cases, as can be seen by inspecting some of the intersections in Figure 6(f). By studying the numbers in Table II, the first conclusion is that this case is quantifiably worse. Not only is the expected quality lower, but the statistical distribution is also flatter. Our proposed algorithms are able to keep their lead, nonetheless. In fact, the lead is even more obvious in this case as visualized in Figure 8; the patched refinement program is able to keep up with the force-equilibrium algorithm on size quality, but looks more like the classical one when it comes to angle quality. 4.5. Simulation performance As the primary use of a mesh is in simulation, we will now check how the different methods compare when used in a reservoir simulator. In the literature, both node-centered finite volume-finite element c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
12
R. HOLM ET AL.
Table II. Quality statistics for the meshes generated from the stochastic fracture network, with about N ≈ 5600 triangles. All statistics have been scaled by 100 and given as percentages.
triangle patched persson
µ
qangle σ
1%
µ
80.7 83.5 88.5
16.5 18.2 14.8
30.5 14.9 22.3
83.3 85.6 89.9
90
qtriangle σ 1% 13.8 16.5 13.5
44.3 17.1 28.4
µ
qsize σ
1%
73.4 91.9 93.4
19.1 10.6 9.7
36.2 51.1 54.9
95
88
90
86
85
84
80
82
75
80
70
78
persson triangle patched 4500
5000
5500
6000
(a) Angle quality, qangle
6500
7000
persson triangle patched
65
4500
5000
5500
6000
6500
7000
(b) Size quality, qsize
Figure 8. Mean quality vs. number of triangles for stochastically generated fractures.
schemes [15] and cell-centered finite volume schemes [10] have been used. Here, we use a cell-centered scheme where the inter-cell fluxes are calculated with a multi-point stencil [1, 2]. For equilateral triangles, this reduces to the conventional two-point method. The test setup is the classical quarter-of-a-five spot: the medium is initially oil filled, and there is one water injector in the lower right corner and correspondingly one producer in the upper right corner. Flux across fractures is reduced by setting a transmissibility multiplier equal to 0.1 over them. And finally, both the oil and water have pressure dependent densities and viscosities, which have been taken from empirical data. We only consider the deterministic fracture network, as the stochastic network was not amenable to practical simulations due to excessively small triangles. The water saturation for the deterministic network is given in Figure 9. As one might expect, the grids yield somewhat different contours, and this may be due to grid orientation effects and accuracy. But as shown in Table III, the water breakthrough occurs at about the same time. Also, the simulation running on the mesh produced by Persson’s approach was somewhat quicker than Triangle or the patched version. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
(a) Triangle
(b) Patched
13
(c) Persson
Figure 9. Water saturation at the end of the simulation.
Table III. Run statistics from the simulation. Time usage is relative to the computer simulation time spent on the Persson grid. The number of failures is the number of times that negative masses are found. Water breakthrough is the earliest time when water is recorded in the producer, given as a percentage of the total runtime.
Relative time
# Failures
Water breakthrough
1.130 1.085 1.000
10 12 8
70.66% 72.52% 73.23%
triangle patched persson
Our simulator uses the implicit pressure, explicit mass formulation. Hence it may have stability problems at larger timesteps. We have tried to use the same timestep in all cases, but this resulted in some failed timesteps due to negative mass numbers. In that case, the timestep had to be retried with a smaller timestep size. Notice that the simulation on the force-equilibrium mesh has fewer failures than the other two methods.
5. CONCLUDING REMARKS Based on the preceding section, we may conclude: • The use of the force-equilibrium algorithm in conjunction with a constrained Delaunay triangulation consistently yielded meshes of high quality. • In comparisons with the state-of-the-art program Triangle, our approach almost always produced higher quality meshes by all conventional measures. • When used as a basis for porous media flow, the force-equilibrium mesh yielded faster solution time and less numerical instability. In summary, the force-equilibrium based meshing excels at moving nodes which results in a high mesh quality. The combination with a constrained Delaunay triangulation ensures perfect fit to internal boundaries, without severely reducing the triangle qualities. Implementing this combined algorithm is also relatively straight-forward. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
14
R. HOLM ET AL.
The logical next step is to extend this work to three dimensions. This primarily requires the development of a constrained Delaunay triangulation in three dimensions [17]. Also, the forceequilibrium algorithm would have to project points from the 3D space onto a 2D plane, but this should be a minor undertaking.
REFERENCES 1. I. Aavatsmark, T. Barkve, Ø. Bøe, and T. Mannseth. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods. SIAM Journal on Scientific Computing, 19:1700–1716, 1998. 2. I. Aavatsmark, T. Barkve, Ø. Bøe, and T. Mannseth. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results. SIAM Journal on Scientific Computing, 19:1717–1736, 1998. 3. F. Aurenhammer. Voronoi diagrams - A survey of a fundamental geometric data structure. ACM Computing Surveys, 23(3):345–405, 1991. 4. A. K. Cline and R. J. Renka. A constrained two dimensional triangulation and the solution of closest node problems in the presence of barriers. SIAM Journal of Numerical Analysis, 27:1305–1321, 1990. 5. H. Edelsbrunner. Geometry and Topology for Mesh Generation. Cambridge University Press, 2001. 6. E. A. Flodin, L. J. Durlofsky, and A. Aydin. Upscaled models of flow and transport in faulted sandstone: Boundary condition effects and explicit fracture modelling. Petroleum Geoscience, 10(2):173–181, 2004. 7. D. Gunasekera, P. Childs, J. Herring, and J. Cox. A multi-point flux discretization scheme for general polyhedral grids. In International Oil and Gas Conference and Exhibition in China. SPE, November 1998. SPE 48855. 8. D. Hale. Atomic meshes: from seismic imaging to reservoir simulation. In Proceedings of the 8th European conference on the mathematics of oil recovery. European Association of Geoscientists & Engineers, 2002. ISBN 90-73781-24-8. 9. C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. Studentlitteratur, 1987. 10. M. Karimi-Fard, L. J. Durlofsky, and K. Aziz. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE Journal, 9(2):227–236, June 2004. SPE 88812. 11. C. L. Lawson. Software for C 1 surface interpolation. In John R. Rice, editor, Mathematical Software III, pages 161–194. Academic Press, 1977. 12. D.-T. Lee and A.K. Lin. Generalized delaunay triangulations for planar graphs. Discrete & Computational Geometry, 1:201–217, 1986. 13. N. E. Odling, S. D. Harris, and R. J. Knipe. Permeability scaling properties of fault damage zones in siliclastic rocks. Journal of Structural Geology, 26(9):1727–1747, 2004. 14. P.-O. Persson and G. Strang. A simple mesh generator in MATLAB. SIAM Review, 46(2):329–345, June 2004. 15. V. Reichenberger, H. Jakobs, P. Bastian, R. Helmig, and J. Niessner. Complex gas-water processes in discrete fracturematrix systems. Upscaling, mass-conservative discretization and efficient multilevel solution. Technischer Bericht Heft 130, IWS, 2004. ISBN: 3-9337 61-22-6. 16. J. R. Shewchuk. Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In Ming C. Lin and Dinesh Manocha, editors, Applied Computational Geometry: Towards Geometric Engineering, volume 1148 of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag, May 1996. From the First ACM Workshop on Applied Computational Geometry. 17. J. R. Shewchuk. Constrained Delaunay tetrahedralizations and provably good boundary recovery. In Proceedings of the Eleventh International Meshing Roundtable (Ithaca, New York), pages 193–204. Sandia National Laboratories, September 2002. 18. J. R. Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry, 22(1-3):21– 74, May 2002. 19. A. Silberhorn-Hemminger. Modeling of fracture aquifier systems; geostatistical analysis and deterministic-stochastic fracture generation. Phd thesis, Universit¨at Stuttgart, January 2003. 20. J. Tveranger, A. Braathen, T. Skar, and A. Skauge. Centre for Integrated Petroleum Research - Research activities with emphasis on fluid flow in fault zones. Norwegian Journal of Geology, 85(1-2):63 – 71, 2005.
c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
Prepared using nlaauth.cls [Version: 2002/09/18 v1.02]
Meshing of domains with complex internal geometries Randi Holm∗ , Roland Kaufmann, Bjørn-Ove Heimsund, Erlend Øian and Magne S. Espedal Center for Integrated Petroleum Research Allgt. 41, N-5007 Bergen, Norway
SUMMARY This paper presents a meshing algorithm for domains with internal boundaries. It is an extension of the gridding algorithm presented by Persson and Strang. The resulting triangulation matches all boundaries, and the triangles are all nearly equilateral. Equilateral triangles are beneficial for a finite volume discretization, as fluid flow between elements of very different size is only possible at small timesteps. The mesh generator is compared with the well regarded Triangle program, where both element quality and simulation performance are checked. It is shown that c 2006 John Wiley & Sons, Ltd. our mesh generator consistently delivers better meshes. Copyright KEY WORDS :
physically-based mesh generation; fractured porous media; subsurface flow simulation
1. INTRODUCTION Flexible meshes are necessary in order to conduct detailed numerical studies of fluid flow and transport in geometrically complex geological systems. Our goal is to generate computational meshes of high quality for fracture systems and layered porous media. We focus on two mesh quality criteria, geometry adaptation and element shape. For general problems these criteria are in conflict, i.e. good geometrical fitting can lead to poor element quality. Recently, a mesh generation algorithm was presented by Persson and Strang [14]. This algorithm iteratively improves an initial mesh by assigning repulsive forces between the nodes on either side of an edge, and then moving the nodes such that global force-equilibrium is obtained. At each step of the algorithm, a new triangulation is obtained by using a Delaunay triangulation method [3, 5]. We present an extension to this algorithm such that the mesh will conform to internal boundaries such as material discontinuities and fractures, in addition to the external boundaries of the domain. The idea of the extension is to restrict the repulsive forces acting on boundary nodes such that the nodes may only move along the boundary, and not leave it. This simple modification uses a constrained Delaunay triangulation method [18] to ensure that no edges of the mesh cross an internal boundary. In the next section we give a geology oriented motivation for this work on unstructured mesh generation. The fundamentals of the meshing algorithm described by Persson [14] are then given in
∗ Correspondence to: [email protected]
c 2006 John Wiley & Sons, Ltd. Copyright
2
R. HOLM ET AL.
Figure 1. The figure shows a domain that contains a stochastically generated set of fractures.
Section 3. That section continues with our extension of the formulation that enables adaptation of the mesh to internal boundaries. This gridding method is then applied to some representative and realistic fracture networks in Section 4, and a comparison to a popular meshing algorithm [16] is made. The comparisons are rounded off by performing some simple water/oil fluid flow simulations.
2. GEOLOGICAL MOTIVATION Geological systems often exhibit complex geometrical structures that complicate the mesh generation. Examples of such include faults with associated fractured damage zones, naturally fractured rocks and layer truncations and pinch-outs in porous media. In particular, faults and fractures in the subsurface may act as barriers or conduits for fluid flow. Consequently, such structures are important factors for determining the flow paths of fluids in porous and non-porous rock masses. In order to get reliable flow simulations, it is then vital to develop robust meshing algorithms that accurately represent the geometrical properties of these features. Important steps in a work flow from problem definition to flow simulation are data acquisition, characterization, petrophysical modeling, meshing and simulation. In this paper we only consider the task of generating a well adapted mesh from some digital characterization of the geometrical structures, and consider its impact on simulations. For concreteness, we motivate the mesh generation task from challenges of meshing faulted and fractured porous rocks. In particular, an ongoing work within our research group is concerned with fluid flow near and in fault zones [20]. Note that by fractures we understand both high permeable open fractures and low permeable deformation bands (“anti-fractures”). The major task is to create a high-quality mesh that is well adapted to these fractures. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
3
Figure 2. A photograph of a meter scale outcrop and the corresponding bitmap of (manually) interpreted layers and fractures/deformation bands.
2.1. Characterization Characterizations of fractured porous media can be categorized into two main classes, a stochastic and a deterministic approach. It is also possible to combine the two into a hybrid stochastic-deterministic approach. The latter is common to apply when particular faults or fractures are known. In the stochastic approach, synthetic fracture systems are generated based on statistical parameters, see Figure 1. The generating algorithms can either be based on purely statistical distributions or on more elaborate schemes that try to mimic geometrical structures of nature [13]. In the stochastic approach, all fractures typically have an analytical geometrical representation. Consequently, these can be used directly as input in the mesh generation procedure. The deterministic approach includes using e.g. outcrop photographs or seismics combined with some image processing technique to recognize the fractures [6, 8]. In Figure 2 below, this is illustrated for a real outcrop example, where the “image processing” technique used is simply manual interpretation by geologists within an image manipulation program. 2.2. Meshing challenges In both characterization approaches described above, there are some common features of most fracture realizations. Reichenberger [15] describes four scenarios which are difficult to handle for mesh generators; see Figure 3. The flow connectivity is affected by slightly crossing, Figure 3(a), and almost crossing fractures, Figure 3(b), and can be important to resolve accurately. It is also challenging to generate high quality elements that adapt to close, (almost) parallel fractures, Figure 3(c), and fractures crossing at a small angle, Figure 3(d). They can all be identified from the example characterizations in Figure 1 and Figure 2 above. On the other hand, the fluid flow simulation performance depends on well-shaped elements. In a control volume formulation, fluxes between adjacent cells are calculated from the cell pressures. For a homogeneous medium consisting of equilateral triangular elements, the flux can be calculated using just the two neighboring cell values [1, 2], which gives a two-point flux stencil. Deviations from an equilateral form implies that multiple cells are coupled for the calculation of fluid fluxes, which is a multi-point flux stencil. In general, multi-point flux stencils lead to more work for flow simulators [7]. The focus of this paper is to present additions to a mesh generation algorithm such that domains with complex internal geometries can be treated and be amenable to flow simulations. Specifically, we seek an algorithm that honors the geometries of faulted and fractured rocks to yield well-formed elements. In this case, fracture end-points and fracture intersections must be accurately represented. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
4
R. HOLM ET AL.
(a)
(b)
(c)
(d)
Figure 3. The figure shows different classes of fracture configurations that are difficult to mesh: (a) slightly crossing fractures; (b) almost crossing fractures; (c) close and (almost) parallel fractures; and, (d) fractures crossing at a small angle.
These are particularly important in order to capture the effects of fluids accumulating in regions limited by low-permeable fractures (compartmentalization) and fluid distributions caused by the connectivity of high-permeable fractures.
3. FORCE EQUILIBRIUM MESHING We will now briefly review the meshing algorithm given by Persson and Strang [14], and then extend it to handle internal interfaces, such as layering or fractures. The algorithm is iterative, where in each step the nodes are moved towards a force equilibrium position and a new topology is formed by a Delaunay triangulation method. The initial condition is an initial point distribution connected by a triangular mesh. The discussion starts with the algorithm for moving the nodes. This is followed by methods for adapting the nodes to internal boundaries along with suitable Delaunay triangulation schemes. 3.1. The edge forces To each edge we assign a force function f (`, `0 ), whose strength depends on the current edge length ` and the target length `0 . As argued by Persson [14], we may set ( `0 − ` ` < `0 , f (`, `0 ) = 0, ` ≥ `0 . Positive forces are repulsive, and hence only edges shorter than their target lengths are assigned a force. Attractive forces could be included by just setting f = `0 − `, but this may reduce the mesh mesh quality. The target length `0 is just a relative quantity that can be used for mesh adaptivity. Its absolute value carries little significance, since the mesh density is determined by the number of nodes in the domain. Consequently, a varying `0 serves to cluster nodes (mesh adaptivity), but will not determine edge lengths. 3.2. The node forces It is the nodes of the mesh which will be moved, and not the edges. Nodes are moved using the summed edge forces of all the edges connected to that node. After moving the nodes, a new set of edges are c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
5
found from a Delaunay triangulation algorithm. The repulsive forces acting on node i is X f ` j , ` j,0 pEi − pEi, j , FEi = j
where the summation is over all the edges connected to node i. pEi is the node coordinate, and pEi, j is the opposing node coordinate along edge j. Since f ≥ 0, the summation into FEi is positively proportional to the vectors pointing at node i. Hence the node forces are positive and repulsive, just as the edge forces, and therefore the nodes will tend to move away from each other. 3.3. An explicit algorithm After assembling all the nodal forces, it remains to solve the coupled system FEi ( pE1 , pE2 , . . . , pEm ) = 0,
for i = 1, . . . , m,
where m is the number of nodes in the mesh. The node forces are dependent on the topology, which in turn has been generated by a Delaunay triangulation. It is consequently hard to solve the equilibrium equations directly since FE is a discontinuous function of the points. Instead, a pseudo time stepping approach is adopted, whereby ∂ pEi = FEi , ∂t
for i = 1, . . . , m,
which is discretized explicitly as pEin+1 = pEin + 1t FEin ,
for i = 1, . . . , m.
(1)
After each time step, new edges are generated from pEin+1 by a Delaunay triangulation algorithm. The iteration terminates once k FEi k ≈ 0, ∀i. 3.4. Delaunay triangulation Perhaps the most common way of triangulating the mesh is the Delaunay scheme. For a Delaunay scheme, triangles are chosen such that no nodes are within the circumcircle, i.e. the circle that goes through all three nodes of the triangle [5]. This is done by picking the nearest neighbor to a node as the one with which it will form an element [4]. An illustration of the concept is given in Figure 4. This small example contains a domain of four nodes. Figure 4(a) shows the circumcircles for each of the four different combinations of three nodes, whereas Figure 4(b) displays the final triangles. As point C is within the circumcircle of AB D, this disqualifies the points A, B and D from forming a triangle. Instead, 4ABC is grown. It has been showed that the Delaunay triangulation gives the largest minimum angle in any element, which creates as equilateral angles as possible [11]. Note however, that it does not necessarily give the largest sum of minimum angles overall. 3.5. Adjusting the nodes to boundaries We now turn to the major topic of the paper: aligning the mesh to boundaries. In our case, there are two types of boundaries: external and internal. External boundaries enclose the domain, while the internal c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
6
R. HOLM ET AL.
A
A
A C
C
C
B
B
B D
D
D
(b) Regular Delaunay
(a) Circumcircles
(c) Constrained Delaunay
Figure 4. Triangulation of four points.
boundaries are the fractures which are to be honored. Points which are pushed outside the domain are simply projected back onto the boundary. For a boxed domain, this projection simply clamps the coordinate to lie within the box. Internal boundaries (fractures) are handled differently. Using a target density, points are initially distributed evenly throughout the domain. The same density is applied along each fracture, and points off the fractures but deemed too close to them are removed to obtain the desired density. Crossing fractures are split at the intersections, which yields a set of non-intersecting fractures. The nodes at the endpoints of the fractures are locked to their position, and the nodal forces acting on them are always set to zero. A point initially placed on a fracture will always be projected back onto that fracture after Equation (1) has been solved. Referring to Figure 5, let the vector uE be the fracture and vE be the vector that connects one endpoint of the fracture to the point to be projected back. From standard vector calculus, the desired vector vEu is vEu =
vE · uE uE. uE · uE
This uses the property that vE · uE = kE v kkE u k cos θ . 3.6. Constrained Delaunay methods While the nodes will now properly obey all geometrical constraints, the resulting Delaunay triangulation can generate edges which cross fractures. To avoid this, we resort to a constrained triangulation where edges connecting nodes along fractures are kept, and the triangles which are generated cannot have edges crossing the fractures. A constrained Delaunay scheme relaxes the demands on the triangles. The triangles may now have other points inside their circumcircle if those points are on the other side of a constraining barrier. Figure 4(c) depicts how the triangulation may change after introducing constraints. Let us assume that it is desirable to have the line AD as an edge because this represents a fracture. Now, the point C is occluded behind this subsegment, and 4AB D is considered legal. Care must be taken to start the algorithm from a boundary, so that a situation where a triangle crosses the boundary does not arise. Triangles are otherwise grown similarly, and keep many of the desirable properties, such as for instance having the largest minimum angle [12], even if they are no longer strictly Delaunay. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
7
vE C C : uE C C : vEu θ Figure 5. Projecting a point at vE back onto its associated fracture uE, yielding vEu .
There are also some other triangulation methods which we will briefly discuss: Conforming Delaunay combines refinement with constraints by adding points along internal boundaries until the fractures are part of the triangulation while simultaneously maintaining the Delaunay property [17]. Similar is the almost Delaunay procedure, which add points at the intersections between edges and constraining segments, thus ensuring that the edges will be split during the next iteration. Since we prefer to adjust the node geometry using the force equilibrium approach, the constrained Delaunay method fits well into our approach. The conforming and almost Delaunay schemes add points, and may disturb the equilibrium.
4. EXPERIMENTS Based on the preceding discussion, we have implemented an augmented version of Persson’s algorithm with the following enhancements: • Internal boundaries are divided into subsegments with approximately the same size as the expected lengths of interior elements. • Interior points that are close to internal boundaries are removed to avoid crowding. • End-points of internal boundaries remain fixed. • Points that originated at an internal boundary are projected back to that boundary after force application. • A constrained Delaunay triangulation method ensures that no edges cross internal boundaries. 4.1. Quality measures When it comes to using a triangulation in a finite volume method, size does matter. Fluid flow between elements of very different size is only possible at small timesteps, which increase simulation time or even worse, makes it hard to get an accurate result, if any at all. To measure the size quality of an element, we compare the deviation of its size from the mean size of all elements, and adjust this so that the quality ends up between zero (for infinitely different element sizes) and one (for perfectly equal element sizes). All terms are squared to penalize elements further away from the mean. Thus, we end up with the expression: qsize = where ai is the size of element i and a¯ = c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
1 n
a¯ 2
P
i
a¯ 2 , − (a¯ − ai )2
ai is the mean size of all elements. Numer. Linear Algebra Appl. 2006; 00:1–14
8
R. HOLM ET AL.
Another way of measuring size quality would be to use the ratio between sizes of adjacent triangles only. However, due to time constraints this has not been investigated. Since the Delaunay method is used to maximize the minimum angle, element quality has traditionally been measured by looking at the angles. We use the same geometrically intuitive measurement as Persson [14], namely the ratio between the smallest circle that encompasses the triangle (the circumcircle) and the largest circle that fits inside the triangle (the incircle). This ratio is divided by two so that it is normalized between zero (for infinitely thin triangles) and one (for equilateral triangles): (h 1 + h 2 − h 3 )(h 1 + h 3 − h 2 )(h 2 + h 3 − h 1 ) . h1h2h3 An alternative quality measure for triangles that is commonly used is given by Johnson [9]: √ 4 3a qtriangle = 2 , h 1 + h 22 + h 23 qangle =
where a is the area and h 1 , h 2 and h 3 the side lengths of the triangle. If h 1 = h 2 = h 3 , the triangle quality is optimal and q = 1. Quality is usually deemed acceptable if qtriangle > 0.6. We study the properties of our approach by meshing domains of different origin and compare with grids generated by a conventional method. The methods will be compared using the quality measures described above in conjunction with visual inspection. For each of the examples we provide statistics on the mean over all cells, the standard deviation and the first percentile of the quality. The mean, µ, will give us indication of the overall quality of the generated mesh, whereas the standard deviation, σ , tells us to which degree that quality is homogeneous throughout the elements. As we will see, a mesh can have a very good overall quality but still have some regions of lower quality. It will typically be those elements that pose problems during a simulation. The first percentile indicates in which range the elements with lowest quality lie. We could also have used the minimum quality found, but in cases with a significant density of fractures there will typically be at least one element of very low quality. By using the first percentile instead, we are less prone to be biased by those pathological cases in our analysis. 4.2. Benchmark As a benchmark to which we can compare our own approach we use Jonathan Shewchuk’s wellrenowned Triangle program, [16]† , which may be considered the state of the art of triangulation. It is our experience that this program in general generates very good triangulations. However, to maintain quality it always refines the initial mesh by splitting elements into smaller ones, thus having a tendency to end up with “cobwebs” around fractures that crosses at a steep angle. In contrast, our approach attempts to minimize the global tension between points, which means that tricky areas are avoided rather than successively refined. Some of these problems may be rectified by altering triangulation parameters. For instance, relaxing the angle quality constraint may in some cases give an overall better quality. These tweaks tend to be very mesh-specific, as the program already does an excellent job for the general case out of the box with the default settings.
† Version 1.6 released on July 28, 2005
c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
9
We have also created a patch to the Triangle program that inhibits the refinement cascade for elements below a certain specified size‡ . Preferring size over angle quality, this method will yield long, thin “needles”. The reader should be informed that this method is without any theoretical guarantees. 4.3. Deterministic The algorithms were used to create a triangulation of the outcrop in Figure 2 that was interpreted manually. We strived to get triangulations that contained about the same number of elements as possible to get comparable results. Figure 6(a) is made by Triangle. Note the little cobweb that appears around the middle of the upper right quadrant. Usually, these clusters indicate that the lines do not quite match up or make steep angle. We have manually inspected this location and there is nothing special about it; it simply seems to be an unfortunate angle. When inhibiting refinement in Triangle (Patched), we get Figure 6(c). At first sight there appears to be no irregularities, but upon closer inspection a lot of thin triangles can be found. None of these are detrimental, but by having so many of them, they represent an impediment for an effective flow simulation. The best results were obtained using the revised force-equilibrium algorithm (Persson), which is shown in Figure 6(e). Most of the thin triangles are gone; those that are left are typically in places that would be hard to triangulate anyhow. Table I shows how the various methods compare to each other. Not surprisingly, the revised forceequilibrium and the patched refinement algorithms are able to generate higher size quality on average, but they also have a generally better distribution of the sizes as well, in terms of lower standard deviation and higher percentile mark. A result that we had not expected was that the revised force-equilibrium and the patched refinement algorithms also perform better when it comes to angle quality. Whereas the Triangle algorithm may produce the best result when viewing the possibilities for one triangle isolated (given the triangulation that has occurred up till this point), it apparently does not translate into finding a good mesh overall. Graphically, the same trend may be shown in a sensitivity plot. Figure 7 shows the qualities for meshes of various sizes. All algorithms are able to get better qualities for more elements, which is natural given that the distance between the constraints increases. However, the force-equilibrium and the patched refinement maintain a steady lead, for both measures. We are unsure of why the patched refinement program manages to do so well for small meshes, and we would advise against reading too much into that result. Also note that the graphs are not strictly increasing; in some cases the quality may actually decrease a bit even if more elements are made. This shows that there is still some sensitivity for the actual configuration itself. However, we think that the trend should be clear. 4.4. Stochastic To test out the algorithms on a relatively random example, we used the program Frac3D [19], developed at the University of Stuttgart, to generate the set of fractures in Figure 1. The program works by drawing fracture length, orientation and placement from statistical distributions.
‡ Currently we use refinement based on the centroid as an estimator of the child triangles’ sizes, which admittedly may be very
inaccurate for bad triangles. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
10
R. HOLM ET AL.
(a) Triangle, N = 3164
(b) Triangle, N = 5663
(c) Patched, N = 3108
(d) Patched, N = 5583
(e) Persson, N = 3134
(f) Persson, N = 5614
Figure 6. Meshes for the deterministic case (left column) and the stochastic case (right column). The fractures are drawn with a larger linewidth. About N ≈ 3100 triangles were used for the deterministic case, and about N ≈ 5600 triangles in the stochastic case. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
11
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
Table I. Quality statistics for the meshes generated from the deterministic fracture network, with about N ≈ 3100 triangles. All statistics have been scaled by 100 and given as percentages.
triangle patched persson
µ
qangle σ
1%
µ
86.5 86.8 90.2
14.1 13.8 11.3
32.0 31.8 48.6
88.7 88.9 91.6
qtriangle σ 1% 11.4 11.2 9.5
46.3 45.9 55.4
µ
qsize σ
1%
94.0 94.8 94.9
7.9 5.9 7.5
57.9 76.0 64.4
92 96 91 95
90
89
94
88 93 87 92 86 91
85
84
90 persson triangle patched
83
82
1500
2000
2500
3000
3500
(a) Angle quality, qangle
4000
persson triangle patched
89 1500
2000
2500
3000
3500
4000
(b) Size quality, qsize
Figure 7. Mean quality vs. number of triangles for a deterministic outcrop.
In this case Triangle performed poorly. As shown in Figure 6(b), many small elements were required to match the fracture set. Most of these are removed and replaced by single triangles in our modification, as depicted in Figure 6(d). In contrast, this is a case where the force-equilibrium variant really shines. There are a lot fewer bad cases, as can be seen by inspecting some of the intersections in Figure 6(f). By studying the numbers in Table II, the first conclusion is that this case is quantifiably worse. Not only is the expected quality lower, but the statistical distribution is also flatter. Our proposed algorithms are able to keep their lead, nonetheless. In fact, the lead is even more obvious in this case as visualized in Figure 8; the patched refinement program is able to keep up with the force-equilibrium algorithm on size quality, but looks more like the classical one when it comes to angle quality. 4.5. Simulation performance As the primary use of a mesh is in simulation, we will now check how the different methods compare when used in a reservoir simulator. In the literature, both node-centered finite volume-finite element c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
12
R. HOLM ET AL.
Table II. Quality statistics for the meshes generated from the stochastic fracture network, with about N ≈ 5600 triangles. All statistics have been scaled by 100 and given as percentages.
triangle patched persson
µ
qangle σ
1%
µ
80.7 83.5 88.5
16.5 18.2 14.8
30.5 14.9 22.3
83.3 85.6 89.9
90
qtriangle σ 1% 13.8 16.5 13.5
44.3 17.1 28.4
µ
qsize σ
1%
73.4 91.9 93.4
19.1 10.6 9.7
36.2 51.1 54.9
95
88
90
86
85
84
80
82
75
80
70
78
persson triangle patched 4500
5000
5500
6000
(a) Angle quality, qangle
6500
7000
persson triangle patched
65
4500
5000
5500
6000
6500
7000
(b) Size quality, qsize
Figure 8. Mean quality vs. number of triangles for stochastically generated fractures.
schemes [15] and cell-centered finite volume schemes [10] have been used. Here, we use a cell-centered scheme where the inter-cell fluxes are calculated with a multi-point stencil [1, 2]. For equilateral triangles, this reduces to the conventional two-point method. The test setup is the classical quarter-of-a-five spot: the medium is initially oil filled, and there is one water injector in the lower right corner and correspondingly one producer in the upper right corner. Flux across fractures is reduced by setting a transmissibility multiplier equal to 0.1 over them. And finally, both the oil and water have pressure dependent densities and viscosities, which have been taken from empirical data. We only consider the deterministic fracture network, as the stochastic network was not amenable to practical simulations due to excessively small triangles. The water saturation for the deterministic network is given in Figure 9. As one might expect, the grids yield somewhat different contours, and this may be due to grid orientation effects and accuracy. But as shown in Table III, the water breakthrough occurs at about the same time. Also, the simulation running on the mesh produced by Persson’s approach was somewhat quicker than Triangle or the patched version. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
MESHING OF DOMAINS WITH COMPLEX INTERNAL GEOMETRIES
(a) Triangle
(b) Patched
13
(c) Persson
Figure 9. Water saturation at the end of the simulation.
Table III. Run statistics from the simulation. Time usage is relative to the computer simulation time spent on the Persson grid. The number of failures is the number of times that negative masses are found. Water breakthrough is the earliest time when water is recorded in the producer, given as a percentage of the total runtime.
Relative time
# Failures
Water breakthrough
1.130 1.085 1.000
10 12 8
70.66% 72.52% 73.23%
triangle patched persson
Our simulator uses the implicit pressure, explicit mass formulation. Hence it may have stability problems at larger timesteps. We have tried to use the same timestep in all cases, but this resulted in some failed timesteps due to negative mass numbers. In that case, the timestep had to be retried with a smaller timestep size. Notice that the simulation on the force-equilibrium mesh has fewer failures than the other two methods.
5. CONCLUDING REMARKS Based on the preceding section, we may conclude: • The use of the force-equilibrium algorithm in conjunction with a constrained Delaunay triangulation consistently yielded meshes of high quality. • In comparisons with the state-of-the-art program Triangle, our approach almost always produced higher quality meshes by all conventional measures. • When used as a basis for porous media flow, the force-equilibrium mesh yielded faster solution time and less numerical instability. In summary, the force-equilibrium based meshing excels at moving nodes which results in a high mesh quality. The combination with a constrained Delaunay triangulation ensures perfect fit to internal boundaries, without severely reducing the triangle qualities. Implementing this combined algorithm is also relatively straight-forward. c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
14
R. HOLM ET AL.
The logical next step is to extend this work to three dimensions. This primarily requires the development of a constrained Delaunay triangulation in three dimensions [17]. Also, the forceequilibrium algorithm would have to project points from the 3D space onto a 2D plane, but this should be a minor undertaking.
REFERENCES 1. I. Aavatsmark, T. Barkve, Ø. Bøe, and T. Mannseth. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods. SIAM Journal on Scientific Computing, 19:1700–1716, 1998. 2. I. Aavatsmark, T. Barkve, Ø. Bøe, and T. Mannseth. Discretization on unstructured grids for inhomogeneous, anisotropic media. Part II: Discussion and numerical results. SIAM Journal on Scientific Computing, 19:1717–1736, 1998. 3. F. Aurenhammer. Voronoi diagrams - A survey of a fundamental geometric data structure. ACM Computing Surveys, 23(3):345–405, 1991. 4. A. K. Cline and R. J. Renka. A constrained two dimensional triangulation and the solution of closest node problems in the presence of barriers. SIAM Journal of Numerical Analysis, 27:1305–1321, 1990. 5. H. Edelsbrunner. Geometry and Topology for Mesh Generation. Cambridge University Press, 2001. 6. E. A. Flodin, L. J. Durlofsky, and A. Aydin. Upscaled models of flow and transport in faulted sandstone: Boundary condition effects and explicit fracture modelling. Petroleum Geoscience, 10(2):173–181, 2004. 7. D. Gunasekera, P. Childs, J. Herring, and J. Cox. A multi-point flux discretization scheme for general polyhedral grids. In International Oil and Gas Conference and Exhibition in China. SPE, November 1998. SPE 48855. 8. D. Hale. Atomic meshes: from seismic imaging to reservoir simulation. In Proceedings of the 8th European conference on the mathematics of oil recovery. European Association of Geoscientists & Engineers, 2002. ISBN 90-73781-24-8. 9. C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. Studentlitteratur, 1987. 10. M. Karimi-Fard, L. J. Durlofsky, and K. Aziz. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. SPE Journal, 9(2):227–236, June 2004. SPE 88812. 11. C. L. Lawson. Software for C 1 surface interpolation. In John R. Rice, editor, Mathematical Software III, pages 161–194. Academic Press, 1977. 12. D.-T. Lee and A.K. Lin. Generalized delaunay triangulations for planar graphs. Discrete & Computational Geometry, 1:201–217, 1986. 13. N. E. Odling, S. D. Harris, and R. J. Knipe. Permeability scaling properties of fault damage zones in siliclastic rocks. Journal of Structural Geology, 26(9):1727–1747, 2004. 14. P.-O. Persson and G. Strang. A simple mesh generator in MATLAB. SIAM Review, 46(2):329–345, June 2004. 15. V. Reichenberger, H. Jakobs, P. Bastian, R. Helmig, and J. Niessner. Complex gas-water processes in discrete fracturematrix systems. Upscaling, mass-conservative discretization and efficient multilevel solution. Technischer Bericht Heft 130, IWS, 2004. ISBN: 3-9337 61-22-6. 16. J. R. Shewchuk. Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In Ming C. Lin and Dinesh Manocha, editors, Applied Computational Geometry: Towards Geometric Engineering, volume 1148 of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag, May 1996. From the First ACM Workshop on Applied Computational Geometry. 17. J. R. Shewchuk. Constrained Delaunay tetrahedralizations and provably good boundary recovery. In Proceedings of the Eleventh International Meshing Roundtable (Ithaca, New York), pages 193–204. Sandia National Laboratories, September 2002. 18. J. R. Shewchuk. Delaunay refinement algorithms for triangular mesh generation. Computational Geometry, 22(1-3):21– 74, May 2002. 19. A. Silberhorn-Hemminger. Modeling of fracture aquifier systems; geostatistical analysis and deterministic-stochastic fracture generation. Phd thesis, Universit¨at Stuttgart, January 2003. 20. J. Tveranger, A. Braathen, T. Skar, and A. Skauge. Centre for Integrated Petroleum Research - Research activities with emphasis on fluid flow in fault zones. Norwegian Journal of Geology, 85(1-2):63 – 71, 2005.
c 2006 John Wiley & Sons, Ltd. Copyright Prepared using nlaauth.cls
Numer. Linear Algebra Appl. 2006; 00:1–14
