ACCURATE STREAMLINE TRACING AND COVERAGE
A REPORT SUBMITTED TO THE DEPARTMENT OF PETROLEUM ENGINEERING OF STANFORD UNIVERSITY IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE
By Sébastien François Matringe June 2004
I certify that I have read this report and that in my opinion it is fully adequate, in scope and in quality, as partial fulfillment of the degree of Master of Science in Petroleum Engineering.
__________________________________ Prof. Margot G. Gerritsen (Principal Advisor)
iii
Abstract Because of their efficiency, streamline methods are increasingly used for, amongst others, flow based upscaling, history matching, optimal well placement, and compositional simulation, for which they show great potential. The accuracy of the streamline solutions depends on several factors including the quality of the pressure solve, the numerical scheme used to propagate solution components along streamlines, the mappings between streamlines and pressure grid, and naturally, the quality of the streamline tracing algorithm. Using the homogeneous quarter five spot pattern, for which the exact streamline pattern is known analytically, we show that the errors in streamline location, computed arc length and time of flight can be significant, and can negatively affect the simulation results. The tracing errors are caused by the inaccuracy of the fluxes at the edges of the gridblocks, the choice of streamline launching pattern, and the order of flux interpolation used to compute the velocity field. We have investigated the use of the Mixed Hybrid Finite Element method for solving the pressure field that is known to provide more accurate fluxes, and discuss when this method leads to superior tracing results over the traditional finite difference method. We also show that a judicious choice of launch locations leads to lower tracing errors, and present two new methods to improve streamline tracing within grid cells. The first is an Adaptive-Mesh-Refinement-inspired tracing method that provides more accurate tracing where needed. The second method is a higher order variation on Pollock’s tracing that uses a second order interpolation of the velocity field and, like Pollock’s tracing, is efficient and flux conservative. Next, the problem of streamline coverage will be studied. To provide the optimal streamline density in the reservoir, an error indicator is provided that measures the accuracy of the mapping algorithm. Using such an error indicator is necessary in the design of an optimal streamline coverage methodology.
v
Acknowledgments My deepest gratitude is addressed to Prof. Margot G. Gerritsen, for her constant time, guidance and support. I also wish to thank Prof. Ruben Juanes for his help and support on the theory of Mixed Finite Element Methods, Bradley T. Mallison for his constant interest and help and Dr. Mathieu Prévost for his early help and guidance on streamline tracing. I wish to acknowledge funding from the Stanford University Petroleum Research Institute – Gas Injection (SUPRI-C) program. I feel extremely lucky to belong to such a talented and dynamic group of individuals. At a more personal level, I wish to thank Nitin Srvivastava for frequent animated scientific discussions late at night and Marie Lasnier for her every day support. I would like to dedicate this work, to my family, Ghislaine, François, Laurence and Hadrien for their support and love.
vii
Contents Abstract.............................................................................................................................. v Acknowledgments ........................................................................................................... vii Contents ............................................................................................................................ ix List of Tables .................................................................................................................. xiii List of Figures.................................................................................................................. xv
1.
Introduction............................................................................................................... 1 1.1. Motivation............................................................................................................... 1 1.2. Streamline simulation ............................................................................................. 1 1.3. Objectives ............................................................................................................... 3
2.
Tracing errors in the streamline method ................................................................ 5 2.1. The current tracing algorithm: Pollock’s method ................................................... 5 2.1.1. A first order interpolation scheme. ................................................................. 5 2.1.2. Finding the streamline location ...................................................................... 6 2.2. Analytical solution .................................................................................................. 7 2.2.1. P and v fields................................................................................................... 8 2.2.2. Finding the streamline location .................................................................... 10 2.3. Comparison of Pollock’s tracing to the analytical solution .................................. 10 2.3.1. Errors in location, arc length and time of flight........................................... 10 2.3.2. Influence of the interpolation scheme ........................................................... 12 2.3.3. Influence of the launching point ................................................................... 14 2.3.4. Influence of the accuracy of the fluxes.......................................................... 16 2.4. Conclusions........................................................................................................... 17
3.
Improving the accuracy of the fluxes .................................................................... 19 3.1. Governing equations ............................................................................................. 19 3.2. Finite difference method ....................................................................................... 20 3.3. Mixed finite element methods .............................................................................. 22 3.3.1. The finite element approach.......................................................................... 22 3.3.2. The mixed FEM............................................................................................. 24 3.3.3. The mixed hybrid FEM ................................................................................. 24 3.3.4. Advantages of MFEM and MHFEM............................................................. 25 ix
3.4. Comparative results .............................................................................................. 26 3.4.1. Examples of high-permeability streaks ......................................................... 27 3.4.2. Examples of flow barriers............................................................................. 28 3.4.3. Examples of real field ................................................................................... 29 3.4.4. Observations ................................................................................................. 29 3.5. When should we use Mixed Hybrid FEM over Finite Difference? ...................... 29
4.
Improving the interpolation method ..................................................................... 31 4.1. Adaptive tracing refinement.................................................................................. 31 4.1.1. Basic concept ................................................................................................ 32 4.1.2. Defining boundary conditions for the fine scale problem ............................ 32 4.1.3. Solving for the fine scale pressure and velocity field ................................... 37 4.1.4. Upscaled permeability field .......................................................................... 38 4.2. Higher order tracing algorithm.............................................................................. 40 4.2.1. Defining the velocity profile on a gridblock face.......................................... 42 4.2.2. Velocity spaces.............................................................................................. 43 A/ Property of the global velocity field ............................................................................43 B/ Defining the global velocity field from local velocity fields on finite elements .........43 C/ Conforming method to approximate H(div, Ω)............................................................44
4.2.3. 4.2.4. 4.2.5. 4.2.6.
Constructing the higher order velocity field................................................. 44 Computation of the stream function.............................................................. 45 Tracing the streamlines................................................................................. 46 Extension of the method to two dimensional unstructured grids.................. 48
A/ Nonorthogonal quadrilateral grids ...............................................................................48 B/ Triangular grids ............................................................................................................49
4.3. Results................................................................................................................... 49 4.4. Conclusions........................................................................................................... 52
5.
Optimizing the streamline launching location to minimize the tracing errors. 53 5.1. Launching from the injectors/producers ............................................................... 53 5.2. Optimized launching location ............................................................................... 54 5.2.1. Launch on a circle a distance away from the singularities .......................... 54 5.2.2. Launch on the Voronoi grid built around the singularities of the domain ... 54 5.3. Comparison on the quarter five-spot..................................................................... 55
6.
Streamline Coverage............................................................................................... 57 6.1. Error indicator....................................................................................................... 58 6.1.1. The Kriging variance - Limitations .............................................................. 58 6.1.2. Global variance - Definition......................................................................... 60 6.2. Coverage control algorithm using the error indicator ........................................... 60 6.3. Partial streamline tracing ...................................................................................... 61
x
7.
Conclusions.............................................................................................................. 63 7.1. Conclusions........................................................................................................... 63 7.2. Future work........................................................................................................... 64
Nomenclature .................................................................................................................. 65 References ........................................................................................................................ 67
xi
List of Tables
Table 4-1:
Time of flight (TOF) error for the three methods and gain in accuracy when using a higher order method over Pollock’s for unfavorable launching locations ....................................................................................50
Table 4-2:
Time of flight (TOF) error for the three methods and gain in accuracy when using a higher order method over Pollock’s for favorable launching locations ....................................................................................................................52
Table 5-1:
Average errors in time of flight for the three launching methods shown in figure 5-3 for a 10x10 grid.........................................................................56
Table 5-2:
Average errors in arc length for the three launching methods shown in figure 5-3 for a 100x100 grid.....................................................................56
xiii
List of Figures
Figure 2-1:
Representation of a streamline through a two-dimensional cell (from Batycky (1997))............................................................................................6
Figure 2-2:
Normalized pressure field for the homogeneous quarter five spot. .............9
Figure 2-3:
Error in time of flight as a function of the number of gridblocks used......11
Figure 2-4:
Error in arc length as a function of the number of gridblocks used...........11
Figure 2-5:
Pollock’s streamlines (blue) vs. analytical streamlines (black) for a 10x10 grid. ............................................................................................................12
Figure 2-6:
Errors in time of flight as a function of the streamline location, for different levels of grid refinement. The numbering of the streamlines is shown in the next figure.............................................................................13
Figure 2-7:
Numbering of the streamlines for figure 2-6..............................................14
Figure 2-8:
Pollock’s streamlines on 100x100 grid (red), Pollock’s streamlines on 10x10 grid and analytical streamlines (black, hidden by the red streamlines) ................................................................................................15
Figure 2-9:
Streamlines launched near the injector, traced on a 100x100 grid (red) vs. analytical streamlines (black).....................................................................16
Figure 2-10: Pollock’s streamlines from analytical pressure field (red); Pollock’s streamlines from analytical velocity field (blue) and Analytical streamlines (black). .......................................................................................................17 Figure 3-1:
Reservoir model used for the comparison..................................................20
Figure 3-2:
Ghost cells added to the domain to account for the Neumann boundary conditions...................................................................................................21
Figure 3-3:
Homogeneous domain of 1md with a permeability streak of 100 mD ......27
Figure 3-4:
Homogeneous domain of 1mD with three permeability streaks of 500 mD ....................................................................................................................28
xv
Figure 3-5:
Homogeneous domain of 200md with two flow barriers of 1 mD ...........28
Figure 3-6:
Homogeneous domain of 200md with two flow barriers of 1 mD ...........28
Figure 3-7:
40x40 field from slice 37 of SPE10 test case ............................................29
Figure 3-8:
40x40 field from slice 59 of SPE10 test case ............................................29
Figure 4-1:
Boundary conditions for the subgrid are determined from the large scale fluxes..........................................................................................................33
Figure 4-2:
Streamline paths obtained with Pollock’s method (blue) and the proposed method (red) with constant subgrid boundary conditions..........................33
Figure 4-3:
Reconstruction of the subgrid fluxes (black arrow) from the neighboring coarse fluxes (blue arrows) ........................................................................34
Figure 4-4:
Defining the slope of the profile from the neighboring fluxes...................35
Figure 4-5:
Ranges (shown in red) of fine scale velocities allowed by the slope limiter for four cases..............................................................................................36
Figure 4-6:
Reconstructed velocity profile at the face (i-½,j) (red profile) ..................37
Figure 4-7:
Presentation of the downscaling problem ..................................................39
Figure 4-8:
Velocity profiles at gridblock faces in Pollock’s method (left) and the higher order method (right)........................................................................41
Figure 4-9:
Linear velocity profile (red) and Pollock’s constant profile (blue)............42
Figure 4-10:
Stream function (variations of red) that yields two possible exit points....47
Figure 4-11:
The slope limiter reduces the efficiency of the higher order methods on the homogeneous quarter five spot. .................................................................50
Figure 4-12:
Streamlines obtained for launching locations optimizing the influence of the slope limiter..........................................................................................51
Figure 4-13:
A close-up on analytical streamlines (black doted), Pollock’s (red), the subgrid streamlines (blue) and the higher order streamlines (green) .........51
Figure 5-1:
The launching points (red dots) are taken on circles around the wells (blue stars)...........................................................................................................54
xvi
Figure 5-2:
Voronoi grid around three different well configurations. ..........................55
Figure 5-3:
Pollock’s streamlines (blue) and analytical streamlines (red) launched from the well block (left), from a circle away from the well (middle) and from the Voronoi grid (right). ....................................................................55
Figure 6-1: Limitations of the Kriging variance as an error estimator...............................59
xvii
Chapter 1
1. Introduction 1.1. Motivation Because of their efficiency, streamline simulation methods are increasingly used in Reservoir Engineering. Current applications range from flow based upscaling (Durlofsky et al. (1996), Prévost (2003)) to history matching (Emanuel and Milliken (2000), Gross (2003)), from injection optimization (Thiele et al.(2003)) to compositional simulation (Jessen and Orr (2003), Thiele et al. (1997)). As streamline methods become more popular, it becomes increasingly important to carefully study their accuracy and efficiency. This is especially true in compositional simulation of gas injection processes, which is the main focus of research in the Stanford University Petroleum Research Institute – Gas Injection (SUPRI-C). The accuracy of the solutions provided by a streamline method depends on many factors including the quality of the pressure solve, the numerical scheme used to propagate solution components along streamlines, the mappings between streamlines and pressure grid, and naturally, the quality of the streamline tracing algorithm. The efficiency of the compositional streamline method is found to be largely impacted by the computational time spent in the one-dimensional solver used to move the compositions forward in time along the streamlines. It is hence crucial for the method to use an optimal number of streamlines. Too few streamlines will obviously impact the accuracy of the solution. Too many streamlines will penalize the efficiency of the method without gain of accuracy in the solution. This observation raises an interesting question: What is the minimum streamline coverage necessary to provide an accurate enough solution? In this work, we focus on improving the accuracy of the streamline tracing algorithm and propose a novel approach to the streamline coverage problem by introducing a coverage control algorithm.
1.2. Streamline simulation We will here recall the main steps of the streamline method. For simplicity, we will describe the streamline method for immiscible two-phase flow in a heterogeneous porous 1
medium. The compressibility, gravity, capillarity, and dispersive forces are here neglected but can be included in the method. The velocity field in a porous media is governed by Darcy’s law through
1
u=−
µ
k .∇( p ) ,
(1-1)
where u is the Darcy flow velocity, p is the pressure, µ the dynamic viscosity of the fluid and k the permeability tensor of the medium. As the system considered is incompressible, the continuity equation, expressing the mass conservation condition, will be
∇.u = q ,
(1-2)
where q is a source or sink term representing the wells of the domain. Combining equations (1-1) and (1-2), we obtain the second order elliptic pressure equation
∇(
1
µ
k .∇( p )) = q .
(1-3)
Writing the mass conservation of the injected fluid yields a hyperbolic equation in terms of the volumetric saturation S of the injected fluid ∂S + u.∇( S ) = 0 . ∂t
φ
(1-4)
Using the operator identity
u.∇ = u
∂ , ∂ξ
(1-5)
where ξ is the arc length along a streamline, we can reformulate the transport equation as
φ
∂S ∂S +u =0 . ∂t ∂ξ
(1-6)
One can then use the time of flight τ , given by
u
∂ ∂ =φ , ∂ξ ∂τ
to obtain a new formulation of the transport equation:
2
(1-7)
∂S ∂S + =0 . ∂t ∂τ
(1-8)
In streamline simulation, equations (1-3) and (1-8) will be solved sequentially. First, a finite difference or finite volume approach will be used to solve (1-3) for the pressure. The velocity field will then be recovered from the pressure nodes. Using this velocity field, streamlines are traced on the domain. The saturations, known initially on the finite difference grid, will then be mapped to the streamlines. Once the saturations are known on the streamlines, a one-dimensional solver is used to move saturations forward in time along the streamlines. At some point, the pressure field and streamlines location need to be updated. The saturations are then mapped back from the streamlines to the finite difference grid. A new pressure solve is now possible. By decoupling in such a way the pressure and transport equations, two timescales can be defined. The first corresponds to the finite difference solve, the second to the transport problem on the streamline. This sequential approach allows large time steps for the finite difference solve. This provides the streamline method its efficiency.
1.3. Objectives The first focus of this work will be on studying the accuracy of the streamline tracing method employed in streamline methods. The time of flight, arc length and streamline location all depend on the streamline tracing algorithm employed. We hence understand the importance of the accuracy of the streamline tracing method for all applications of streamline simulation. In Chapter 2, the tracing errors introduced by Pollock’s method will be studied. Methods will then be provided to reduce the three sources of tracing errors identified in Chapter 2. Thus, Chapter 3 will study the use of Mixed Hybrid Finite Element method to improve the accuracy of the Darcy fluxes at the faces of the gridblocks. Then, improvements on the traditional tracing method developed by Pollock (1988) will be proposed in Chapter 4. To finish our study on streamline tracing, we will see in Chapter 5 how to choose optimal launching location to reduce as much as possible the tracing errors. We then focus on an efficiency problem in the streamline method: the streamline coverage. In compositional streamline simulation, the streamline coverage significantly impacts the computational efficiency of the method. We will, in Chapter 6, provide a coverage control algorithm that uses an error indicator to provide an optimal coverage for a user-defined level of accuracy. We also, in this chapter, present the concept of partial streamline tracing for improved coverage results.
3
Chapter 2
2. Tracing errors in the streamline method The most commonly used method for streamline tracing is Pollock’s method. Pollock (1988) initially developed a tracing method for Cartesian grids. He used a first order interpolation of the velocity field within a gridblock to trace streamlines. His method is hence consistent with the order of accuracy of the pressure field. Most research in streamline tracing is focused on extending Pollock’s tracing method to different types of grids using the same low order accuracy as Pollock’s tracing. Cordes and Kinzelbach (1992) and Prévost (2000 and 2003) thus proposed extensions of Pollock’s method to trace streamlines on non-orthogonal and unstructured grids. The focus of this chapter will be to study the numerical errors present in the streamline tracing method proposed by Pollock. The objective is to quantify the importance of the tracing errors for streamline simulation applications. To our knowledge, all commercial streamline simulators are today relying on Pollock’s method to trace streamlines. The attraction of this method relies on the fact that it provides flux conservative streamlines and that it is formulated in term of time of flight, which is, as we saw in Chapter 1 a key parameter in streamline simulation. We will here recall Pollock’s method. Then, the analytical solution of a quarter five-spot problem will be presented. Finally, we will compare the results of Pollock’s tracing method to the exact streamline locations given by the analytical solution and analyze the results. We identify the important criteria for tracing accuracy and isolate three possible sources of tracing errors.
2.1. The current tracing algorithm: Pollock’s method 2.1.1. A first order interpolation scheme. A streamline is the path in a flow field that is at all times tangent to the velocity vector. The starting point of all streamline-tracing algorithms is hence to define a velocity field. In Pollock’s method, the velocity field is linearly interpolated from the fluxes at the faces of the gridblock, according to
5
v x = v x ,o + a x .( x − xo )
, ax =
v y = v y ,o + a y .( y − yo ) , a y = v z = v z ,o + a z .( z − z o )
,
v x ,∆x − v x ,o ∆x v y , ∆y − v y , o
∆y v − v z ,o a z = z ,∆z ∆z
, ,
(2-1)
,
where vx, vy and vz are the x-, y- and z- components of the velocity field, respectively, and vx,o, vx,∆x, vy,o, vy,∆y, vz,o and vz,∆z are the fluxes through the faces of the gridblock. Figure 21 is a schematic representing these variables in a two-dimensional case on a Cartesian grid.
Figure 2-1: Representation of a streamline through a two-dimensional cell (from Batycky (1997))
2.1.2. Finding the streamline location The velocity field being defined, the path of a particle entering the gridblock can now be found analytically. As the velocity in each direction is only a function of the coordinate in that direction, we can write v x = v x ,o + a x .( x − xo ) =
∂x dx = , ∂t dt
(2-2)
with similar equations for the y- and z- directions. Rewriting (2-2) gives
dt =
6
v x ,o
dx . + a x .( x − xo )
(2-3)
Integrating the above equation between the abscissa of the entry point xi of the streamline and the abscissa of the possible exit face xe, we have ∆t e , x =
v 1 . ln( x ,e ) , ax v x ,i
(2-4)
where vx,i is the velocity of the particle in the x- direction when entering the gridblock and vx,e is the velocity in the x- direction at the possible exit face. Similarly, we can compute the time required to reach the other faces of the gridblock. The face that will be effectively reached will then be the one corresponding to the minimum positive time (a negative time corresponding to a streamline traced upwind). This time will be the time of flight of the particle in the block considered. This time of flight is then used to compute the coordinates of the point where the streamline exits the gridblock, which are 1 (v x ,i exp(a x .∆te )) + xo ax 1 xe = (v x ,i exp(a x .∆te )) + xo ax 1 z e = (vz ,i exp(a z .∆t e )) + z o az xe =
, ,
(2-5)
.
Some specific cases are to be taken care of when coding this method such as the cases where both fluxes in one direction are pointing inwards or outwards. Details about these special cases can be found in Pollock (1988).
2.2. Analytical solution
To be able to quantify the accuracy of Pollock’s tracing, an analytical model was chosen. Morel-Seytoux (1965 and 1966) found the analytical solution of a variety of injection problems, including the quarter five spot pattern we are interested in. The assumptions necessary to derive the analytical solution are that of tracer flow in a two-dimensional homogeneous domain. The displacement is considered piston-like and the capillary and gravity effects are neglected.
7
2.2.1. P and v fields
To obtain an analytical solution, Morel-Seytoux (1965) makes use of the complex potential w(z), where z is the complex argument defined by z = x + i. y .
(2-6)
The complex potential is found by locating point sources representing the wells, in the complex space. A variety of complex potentials are provided in Morel-Seytoux (1966). Having the complex potential, both the flow potential Φ and the stream function Ψ can be recovered easily through Φ ( x, y ) = Re( w( z )) ,
(2-7)
Ψ ( x, y ) = Im(w( z )) .
In our case, the complex potential was found to be w( z ) =
1 ln(cn( z + K, ko )) , 2π
(2-8)
where cn(z,ko) is the elliptic cosine of modulus ko, and K=K(ko) is the complete elliptic integral of the first kind. A description of these functions can be found in Milne-Thomson (1965) or Byrd and Friedman (1971). The expressions of the flow potential and the stream function are now found to be respectively, 2 2 1 1 − cn x cn y Φ ( x, y ) = ln 4π ko2 cn x2 + ko '2 cn 2y
, (2-9)
and
Ψ ( x, y ) =
sn dn cn 1 arctan y y x 2π sn x dn x cn y
.
Here, cnx=cn(x,ko), cny=cn(y,ko) and ko ' is the complementary modulus. In our case, k o = ko ' =
1 . 2
(2-10)
Morel-Seytoux shows that the pressure p is recovered from the flow potential through p = −β
8
µΦ kh
q + constant .
(2-11)
So, the pressure is a linear function of the flow potential. In equation 2-11, µ represents the dynamic viscosity, k the permeability of the medium, h the thickness of the reservoir, q the flow rate and β is a conversion factor dependent on the unit used. To obtain the velocity field (vx,vy) from the stream function, one can use
∂Ψ ∂y ∂Ψ vy = − ∂x vx =
, (2-12) .
Differentiation of the stream function yields
vx =
2 2 2 2 2 2 cnx dn x sn x sn y dn y + cn y (dn y − k sn y ) . 2π (cn y dn x sn x ) 2 + (cnx dn y sn y ) 2
vy =
sn dn + cn (dn − k sn ) (cn y dn x sn x ) 2 + (cn x dn y sn y ) 2
,
(2-13) cn y dn y sn y 2π
.
2 x
2 x
2 x
2 x
2
2 x
.
The normalized pressure field obtained for the analytical homogeneous quarter five-spot is shown in figure 2-2.
Figure 2-2: Normalized pressure field for the homogeneous quarter five spot.
9
2.2.2. Finding the streamline location
To obtain the analytical streamlines from this velocity field, we will use the stream function, which is constant along a given streamline. So the equation of a streamline passing through a point (xo,yo) is Ψ ( x, y ) = Ψ ( xo , y o ) .
(2-14)
The inversion of this stream function allows for an analytical expression of the streamline path, given by y = F (arccos( − β + β 2 + 1 ) | α ) ,
(2-15)
where F designates the Elliptic integral of the first kind (Tee (1997), Milne-Thomson (1965), Byrd and Friedman (1971)), and with
α=
β=
π 4
,
sn x dnx tan(2π .Ψ ( xo , yo )) . cnx
(2-16)
(2-17)
The time of flight and the arc length of a streamline are calculated using the analytical expression of the velocity field. The integration involved is done numerically using a Matlab function based on a high order method that relies on an adaptive Gauss-Lobatto quadrature rule described by Gander and Gautschi (2000).
2.3. Comparison of Pollock’s tracing to the analytical solution
We here compare the results of Pollock’s tracing algorithm to Morel-Seytoux’s analytical solution. The aim of this part will be to compute typical tracing errors and to isolate each possible source of tracing error in the streamline method. 2.3.1. Errors in location, arc length and time of flight
For use in a streamline simulation framework, the three most important variables in the tracing method are the streamline location, arc length and time of flight. They define the streamline grid used to solve the transport problem. Therefore, errors in these parameters are expected to affect the accuracy of the simulation results. For example, inaccuracy in the location of a streamline can result in that streamline crossing the wrong gridcells in the domain. The consequence will be that the simulation will map saturations away from their real location. The accuracy of the time of flight is of great importance for this is the parameter used to compute the flux along a streamline. The time of flight variable will also have a direct impact on the breakthrough time of a specific streamline. For
10
applications such as streamline-based history matching, the accuracy of the breakthrough time is crucial. It is used to identify the regions of the reservoir in which the geology will be updated in the history matching process. Figures 2-3 and 2-4 show the errors in, respectively, time of flight and arc length. We can see that for a 10x10 grid, errors of the order of 10% are observed in time of flight and arc length calculations. About 900 gridblocks (30x30 grid) are needed to decrease the error in the time of flight below 1%. To obtain the same accuracy in the arc length calculation, 1600 gridblocks (80x80 grid) are needed.
Figure 2-3: Error in time of flight as a function of the number of gridblocks used.
Figure 2-4: Error in arc length as a function of the number of gridblocks used.
11
2.3.2. Influence of the interpolation scheme
In this subsection, we isolate errors due to the first order interpolation algorithm. To do so, we provide Pollock’s algorithm with the analytical Darcy fluxes at the faces of the gridblocks. Therefore, the tracing errors observed here can only be attributed to the velocity field being approximated by a bilinear function in x- and y. Figure 2-5 shows a plot of both Pollock’s and the analytical streamlines for a 10x10 grid. The blue streamlines are Pollock’s and the black ones are the analytical streamlines.
Figure 2-5: Pollock’s streamlines (blue) vs. analytical streamlines (black) for a 10x10 grid.
In figure 2-5, the launching points of the streamlines are indicated by the circles around the first gridblock of the domain. The blue dots along Pollock’s streamlines mark the points where the streamlines are entering and leaving a gridblock computed by Pollock’s streamline tracing algorithm. Figure 2-5 clearly shows the tracing errors due to the first order interpolation assumption. We can even notice that many streamlines do not pass through the same gridblocks. As explained earlier, this is expected to impact the accuracy of the simulation. This particular result motivates a search for a more accurate interpolation algorithm than Pollock’s. In Chapter 4 we propose two alternative tracing methods.
12
Figure 2-6 is a plot of the errors in time of flight for the ten streamlines shown in Figure 2-7. Each line in Figure 2-6 represents a level of refinement of the grid. We can observe that for coarse grids, where the tracing errors are the largest, there exists a trend in the errors: the streamlines located further away from the diagonal of the quarter five spot are showing the largest errors. This shows the limitations of a first order interpolation scheme. As we move away from the diagonal, streamlines bend more. Therefore, a first order interpolation becomes less accurate. On the contrary, close to the diagonal, the streamlines are almost straight lines from the injector to the producer. In this case, a linear interpolation is reasonable.
Figure 2-6: Errors in time of flight as a function of the streamline location, for different levels of grid refinement. The numbering of the streamlines is shown in the next figure.
13
Figure 2-7: Numbering of the streamlines for figure 2-6.
2.3.3. Influence of the launching point
In all previous cases, the launching points of the streamlines were kept constant when refining the grid and were evenly distributed around a square of length 0.1. In this case, we can show convergence of Pollock’s method as the grid is refined. Figure 2-8 represents this classical refinement study. The launching points of the streamlines are represented by the black circles, the blue lines are the streamlines traced on a 10x10 grid. The red and black streamlines, which are respectively the numerical streamlines obtained with a 100x100 grid and the analytical streamlines, are almost identical. This shows the convergence of Pollock’s method.
14
Figure 2-8: Pollock’s streamlines on 100x100 grid (red), Pollock’s streamlines on 10x10 grid and analytical streamlines (black, hidden by the red streamlines)
However, if we now proceed to a refinement without keeping the launching point constant, a very interesting phenomenon arises. Figure 2-9 shows seven streamlines traced using 100x100 grid, but here, the launching points were taken around the closest gridblock to the injector. That is, around a square of length 0.01. Comparing the streamlines obtained from both refinement methods, we observe that for a similar level of refinement, the accuracy of the streamlines is drastically influenced by the launching point: when launching around the first gridblock next to the injector, the tracing errors remain significant regardless of the fineness of the grid. On the other hand, if the launching points are kept at the same location, a refinement of the grid will eliminate the tracing errors. Chapter 5 will focus on the influence of the launching point on the tracing errors. A method will there be proposed to minimize the tracing errors due to inappropriate streamline launching location.
15
Figure 2-9: Streamlines launched near the injector, traced on a 100x100 grid (red) vs. analytical streamlines (black).
2.3.4. Influence of the accuracy of the fluxes
We will here present the last source of tracing errors we were able to isolate. Figure 2-9 illustrates the problem. This figure is a plot of three different kinds of streamlines. As usual, in black, the analytical streamlines are shown. The blue curves represent Pollock’s streamlines traced using the analytical fluxes of the faces of the gridblocks. The streamlines in red are computed using a Pollock’s interpolation algorithm as well, but instead of using the analytical fluxes, the fluxes are computed directly from the analytical pressure field. The differences observed in figure 2-9 between the blue and the red streamlines are solely due to the difference in accuracy of the fluxes at the cell faces. As Pollock’s algorithm uses these fluxes to reconstruct the velocity field within a gridblock, inaccuracy in these fluxes will result in inaccuracy in the velocity field and hence, inaccuracy in the streamline tracing. In Chapter 3, we will use the Mixed Finite Element methods instead of a classical Finite Difference approach to see if we can recover more accurate Darcy fluxes and hence improve the tracing of the streamlines.
16
Figure 2-10: Pollock’s streamlines from analytical pressure field (red); Pollock’s streamlines from analytical velocity field (blue) and Analytical streamlines (black).
2.4. Conclusions
Using a simple reservoir model, analytical streamlines were drawn and used as a reference to test the accuracy of the streamlines traced with the first order tracing algorithm developed by Pollock (1988). Tracing errors were observed in terms of streamline location, time of flight and arc length. Comparing Pollock’s tracing algorithm to the analytical solution, three sources of errors were identified: inaccuracy of the Darcy fluxes, low order of accuracy of the interpolation of the velocity field, and location of the launching point. The three next chapters will focus on finding ways to reduce the tracing errors due to these three sources.
17
Chapter 3
3. Improving the accuracy of the fluxes In the current tracing algorithm, the velocity field used to trace streamlines is reconstructed from the fluxes through the faces of the gridblock. We will in this chapter investigate whether Mixed Finite Element Methods can improve the accuracy of these fluxes. For unstructured grids, Mosé et al. (1994) and Durlofsky (1994) have discussed the benefits of using a Mixed Finite Element Method over a Control Volume Finite Element method (CVFE), a method similar to an integrated FD method. Here, we will first highlight the differences between the Finite Difference, Finite Element (FEM), Mixed Finite Element (MFEM), and Mixed Hybrid Finite Element (MHFEM) methods. Then, comparative results of Finite Difference and Mixed Hybrid FEM will be shown and analyzed.
3.1. Governing equations
We will be using the Finite Difference method to solve a flow problem in a heterogeneous porous domain. If no sources or sinks are present in the domain, and if the system is considered incompressible, the continuity equation, expressing the mass conservation condition, will be
∇.u = 0 ,
(3-1)
where u is the Darcy velocity. We can relate the velocity field to the pressure field using Darcy’s law
u=−
1
µ
k .∇( p ) ,
(3-2)
where µ is the dynamic viscosity of the fluid and k the permeability tensor of the medium. Combining equations (3-1) and (3-2), we obtain the second order elliptic pressure equation
19
1 ∇( k .∇( p )) = 0 .
µ
(3-3)
To compare both methods, we will use a square reservoir in two dimensions. We will use a diagonal permeability tensor and the boundary conditions shown in Figure 3-1. No flow BC
Neumann BC
Neumann BC
No flow BC
Figure 3-1: Reservoir model used for the comparison
We take a unit mobility ratio and test different permeability fields to visualize the accuracy of both algorithms.
3.2. Finite difference method
The Finite Difference method (FD) uses a discretization of the strong form of the governing equation. For a two-dimensional problem and assuming a diagonal permeability tensor the governing equation (3-3) will be transformed to yield
∂2 p ∂2 p kx 2 + k y 2 = 0 . ∂x ∂y
(3-4)
Now, approximating the continuous pressure p by its discrete values ( pi , j ) at the centers of a Cartesian grid, the previous equation can be rewritten using a central discretization scheme Τxi −1/ 2, j .( p i −1, j − pi , j ) + Τxi+1/ 2, j .( p i +1, j − pi , j ) + Τyi , j −1/ 2 .( p i , j −1 − pi , j ) + Τyi , j +1/ 2 .( p i , j +1 − pi , j ) = 0 , (3-5)
where Τxi−1 / 2, j is the transmissibility in the x- direction at face (i-½,j) which forms the interface between block (i-1,j) and block (i,j). This transmissibility can be computed using a weighted harmonic average of the permeability values of the two neighbors in the form 20
Τxi −1/ 2, j =
2.ki , j .ki−1, j ∆y . . ki , j + ki−1, j ∆x
(3-6)
This averaging of permeabilities is of crucial importance for us as it is the way the continuity between the velocities at the two gridblocks is enforced. Once the problem is discretized, we solve for the pressure by writing the system of equation in its matrix form A. p = b ,
(3-7)
where A is the transmissibility matrix, containing the coefficients multiplying the pressures, p is the vector containing the discrete values of the pressure at the centers of the gridblocks, and b is the right hand side vector accounting for the boundary conditions. As explained earlier, we will use in this particular case Neumann type boundary conditions on the whole boundary of the domain. Neumann boundary conditions render the system singular. An extra condition must be enforced to find a unique solution. Here, we use a mass conservation condition on the domain, which is, for an in incompressible fluid, the flux conservation equation Flow rate in = Flow rate out .
(3-8)
When this condition is enforced, a solution (defined up to a constant) will exist. To implement the Neumann boundary conditions, we will use the method of ghost cells. In this method, the domain is extended by adding a layer of blocks around the boundaries of the domain, as shown in figure 3-2. Previous domain Ghost cells added to domain
Figure 3-2: Ghost cells added to the domain to account for the Neumann boundary conditions
21
We will then use the pressure values at these ghost cells to express the Neumann boundary conditions as
VNeumannBC =
( p ghostcell − pboundary ) ∆x
.
(3-9)
Using (3-9) the ghost cell pressures are eliminated from the set of unknowns. Finally, the system can be solved. In our case, the system is symmetric, positive definite and well posed. We solve it using Matlab’s direct solver, which employs a sparse LU factorization algorithm.
3.3. Mixed finite element methods
We have implemented a Mixed Hybrid Finite Element Method (MHFEM) to investigate the importance of the accuracy of the fluxes for the purpose of streamline tracing. Here, we will describe the principles of the Finite Element Method (FEM), omitting the details, which can be found in Hughes (1987) or Brenner and Scott (1994). We will then present the approach of the Mixed FEM. Finally, we will present the Mixed Hybrid FEM as a reformulation of the Mixed FEM that yields a positive definite system of equations. 3.3.1. The finite element approach
The Finite Element Method (FEM) is fundamentally different from the Finite Difference method. In FD, a discretization of the strong form of the governing partial differential equation is used. In FEM, we start by transforming the governing partial differential equation. The main steps of the FEM are succinctly described below. •
We start with the strong form of the PDE
∇(−k .∇( p )) = 0 ,
p∈S ,
(3-10)
where p, the solution of the problem on the domain Ω, belongs to a space of admissible functions that we will call S. •
We then multiply by a test function v belonging to S and integrate over the volume of the domain, giving
∫ ∇(−k .∇( p)).v.dω = 0
Ω
22
, ( p, v ) ∈ S 2 .
(3-11)
•
An integration by parts yields
∫ v.(−k .∇( p)).n.dΓ + ∫ ∇(v).k .∇( p).dω = 0
∂Ω Neumann
.
(3-12)
Ω
• Equation (3-12) is the weak form of the problem. Using the boundary conditions we can write
∫ v.(−k.∇( p)).n.dΓ = ∫ v.Q
out
∂Ω Neumann
.dΓ ,
(3-13)
.dΓ .
(3-14)
∂Ω Neumann
where Qout is the flux out through the boundary of the domain. •
Inserting (3-13) in (3-12) yields
∫ ∇(v).k .∇( p).dω = − ∫ v.Q
out
Ω
•
∂Ω Neumann
We now define the bilinear form a(.,.)as a(v, p) = ∫ ∇(v).k .∇( p).dω ,
(3-15)
∫ v.Q
(3-16)
Ω
and the linear form l(.) as:
l (v ) = −
out
.dΓ ,
∂Ω Neumann
to re-write the weak form of the problem as a ( v, p ) = l ( v )
.
(3-17)
• We now restrict the space S to a finite-dimensional subspace Sm of admissible functions pm, that are solution to the problem on the element Ωm. Using a test function vm defined on Sm, we can write a ( p m , v m ) = l (v m ) . •
(3-18)
If we now define a basis w1, w2,… wn of Sm, we can write pm as n
pm = ∑ α i wi ,
(3-19)
i =1
23
where α i ∈ ℜ, i = 1,..., n . •
We can then re-write the equation (3-18) as n
∑a
i, j
αi = l j ,
j = 1,..., n ,
(3-20)
i =1
with ai , j = a ( wi , w j ) and l j = l ( w j ) .
(3-21)
So the problem is reduced to the solution of the linear system (3-20). We note that the boundary conditions are directly integrated in the formulation. Therefore, an appropriate choice of the base functions allows us to decide in what space the solution will be found. 3.3.2. The mixed FEM
We will here briefly present the Mixed FEM. An extended description of the Mixed FEM can be found in Brezzi and Fortin (1991) or Chavent and Roberts (1991). Like the Finite Element Method, the Mixed FEM solves the weak form of the problem using basis functions. The main difference with the “classical” FEM is that, in MFEM, the set of unknowns is extended by adding the Darcy fluxes at the faces of the gridblocks. To solve for these fluxes, one has to add the divergence-free condition on the velocity field to the weak form of the governing pressure PDE. The natural choice of unknowns for this problem seems to be the pressures at the center of the gridblocks and the fluxes at the faces. Unfortunately, this choice of unknowns yields a system that is not positive definite, hence, hard to invert with “usual” algorithms. This is a major drawback of MFEM: not only the number of unknowns increases significantly but the system is also much harder to invert. 3.3.3. The mixed hybrid FEM
The Mixed Hybrid FEM only differs from the MFEM by the choice of the unknowns of the problem. The equation we will solve will be the same as the one of the MFEM but they will be re-written in terms of the average value of the pressure at the faces of the gridblock. Changing the unknowns will in fact render a system of equations that is positive definite. This will allow for a fast inversion of the matrix of the problem. Hence the use of a Mixed Hybrid Finite Element Method will allow the same accuracy as the MFEM but will provide a system easier to solve for.
24
3.3.4. Advantages of MFEM and MHFEM
The use of MFEM or MHFEM allows the evaluation of the fluxes at the faces of the gridblocks. So both the pressure and the velocity field are simultaneously evaluated. The pressure and velocity fields are different in nature. Indeed, while abrupt changes in pressure can occur in reservoir, the Darcy velocity field will be an order smoother. The idea behind the MFEM approach is then to define appropriate basis functions that will give the correct smoothness properties to the pressure or velocity fields. We will now develop in more details the method we implemented. We will in particular present how the correct smoothness conditions are applied to the pressure and velocity fields. Some details of the following development will be omitted. A full description of the method is available in Chavent and Roberts (1991). The starting point of the method is the set of equations:
∇.u = 0
, (3-22)
u=−
1
µ
k .∇( p)
.
The problem is closed with Neumann boundary conditions. The global velocity field on this domain will have to present a degree of smoothness defined by the space H(div,Ω), given by H (div, Ω) := {u | u ∈ (L2 (Ω)) 2 ; ∇(u ) ∈ L2 (Ω)} .
(3-23)
Where the space L2(Ω) designates the space of square integrable functions on a domain Ω. To be able to obtain this smoothness on the global velocity field, we will force the velocity field within each gridblock to belong to the Raviart-Thomas space on the gridblock. This space will be defined through the set of basis functions {wi | i=1,…,n} defined through:
∫ ∇w
i
=1
,
i = 1,..., n
,
gridblock K
(3-24)
∫ w .v = δ j
ij
.
face Ai
where δij is the Kroenecker delta. So these four basis functions represent a unit flux across a given edge and a zero flux along any other edge. Moreover, these functions will have a unit divergence when integrated over a gridblock.
25
Using these basis functions, the velocity field in a block is expressed as a linear combination of the wi functions. Our choice of basis functions is such that coefficients of this linear combination will be the fluxes across the faces of the gridblock: 4
u K = ∑ Q K , j .w j ,
(3-25)
j =1
where K designates the gridblock and j the face of this gridblock. The term QK,j designates the flux across the edge j of the block K. Three kinds of unknowns can be defined for this problem: PK ∈ R
, pressure at block centers ,
TPK , j ∈ R ,
j = 1,..., n , pressure on faces
,
QK , j ∈ R
j = 1,..., n , flux through faces
.
,
(3-26)
As explained earlier, in a Mixed FEM, the unknowns would be the set {{PK},{QK,j}}, which are the “natural” unknowns of the problem. In a Mixed Hybrid FEM, the set {TPK} will be chosen and this judicious choice of unknowns will yield a system of the form ( M − N ).TP = G ,
(3-27)
where M and N will be symmetric and definite positive matrices, and where G will contain the Neumann boundary conditions of the problem. So even though the number of unknowns compared to a Finite Element or Finite Difference method is much larger, no problem will be encountered to invert the matrices. Once the vector TP of pressures on the faces of the gridblocks is determined, the pressure at the centers as well as the fluxes through the edges can be recovered (Chavent and Roberts (1991)).
3.4. Comparative results
We implemented the Finite Difference method and the Mixed Hybrid Finite Element method to be able to assess the quality of the velocity field in both cases. We will here show comparative results of both methods. It is expected to see identical results as the methods have been proved to reduce to the same order of accuracy in such a simple case as a two dimensional Cartesian grid, the only difference being that the MHFEM is evaluating the fluxes when FD is reconstructing them. To assess the quality of the velocity field obtained from both methods we will be looking at the streamlines traced from both velocity fields. The two methods were implemented in
26
Matlab and are part of the Streamline Simulator “Slsim” developed by Bradley Mallison and Sébastien Matringe with contributions from Mathieu Prévost. The code is using our FD or MHFEM solvers to obtain the pressure and velocity fields. Streamlines are then launched from the left boundary of the reservoir and traced downstream. The quality of the velocity fields obtained was tested on a variety of permeability fields. We have here chosen strongly heterogeneous test cases to be able to see differences between the results of the MHFEM, that is evaluating the fluxes, and the FD, that is reconstructing them. We here tried to present the worst possible cases in terms of flux accuracy, given an orthogonal grid and diagonal permeability tensor. The results presented hereafter will all be in the same format: the first row of plots will be obtained from the FD results, the second line from the MHFEM. Three plots will be provided for each method and each permeability field: the first one on the left will be a contour of the pressure field, the second one represents the streamlines traced and the third one is a saturation map. 3.4.1. Examples of high-permeability streaks
The first kind of permeability field we have tested is composed of a homogeneous domain of low permeability with one (Figure 3-3) or several (Figure 3-4) high permeability streaks.
Figure 3-3: Homogeneous domain of 1md with a permeability streak of 100 mD
27
Figure 3-4: Homogeneous domain of 1md with three permeability streaks of 500 mD
3.4.2. Examples of flow barriers
The second kind of heterogeneity we will present is that of flow barriers. In this case, the domain is of high permeability but flow barriers force the streamlines to bend significantly.
Figure 3-5: Homogeneous domain of 200md with two flow barriers of 1 mD
Figure 3-6: Homogeneous domain of 200md with two flow barriers of 1 mD
28
3.4.3. Examples of real field
Figure 3-7: 40x40 field from slice 37 of SPE10 test case
Figure 3-8: 40x40 field from slice 59 of SPE10 test case
3.4.4. Observations
A first look at the streamlines obtained from both velocity fields is sufficient to realize that, as expected, the difference in accuracy between the two methods is minimal. Most of the streamlines of both methods are almost identical. We do, however, observe a degradation of the streamlines quality of the Finite Difference method for highly heterogeneous cases. In the Mixed Hybrid FEM, the fluxes are computed with the pressure nodes. We are hence guaranteed to obtain the same level of accuracy in the pressure and the velocity field. In FD, the fluxes are reconstructed from the pressure nodes. This reconstruction loses accuracy as the heterogeneity increases.
3.5. When should we use Mixed Hybrid FEM over Finite Difference?
The comparative results showed very little difference between the two methods for reservoir of mild heterogeneity. For more heterogeneous reservoir, like slices of the
29
SPE10 test case, small differences were identified, showing the limitations of the flux reconstruction method used in FD. MHFEM and FD were shown (Durlofski (1994)) to reduce to the same level of accuracy for such a simple grid. Nevertheless, we saw, on strongly heterogeneous domain, the limitations of the flux reconstruction method used in FD. MHFEM results were shown slightly more accurate than FD. Such a method however needs many more unknowns to reach this level of accuracy. Therefore, for a Cartesian grid, with diagonal tensor permeability, the use of mixed finite element methods does not seem necessary. The accuracy of the fluxes may however be improved by using MHFEM when full tensor permeability or unstructured grids are used. Arbogast et al. (1997) and Klausen and Russel (2003) provided an extension of the MFEM for use with full tensor permeability. They showed the improvement in flux accuracy provided by the Expanded Mixed Finite Element Method over a FD 9-point stencil. Durlofsky (1994), Mosé et al. (1994) and Cordes and Kinzelbach (1996) showed that MFEM, FD and CVFEM were providing the same degree of accuracy on triangle-based unstructured grids. The MFEM results however, were shown to maintain their accuracy for strongly heterogeneous reservoir, where the other methods showed their limitations.
30
Chapter 4
4. Improving the interpolation method The previous chapter focused on improving the accuracy of the Darcy fluxes at the faces of the gridblocks. Here, we will propose higher order methods to improve the interpolation made from these fluxes to obtain the velocity field. The tracing method currently used in all commercial simulators was discussed in Section 2.1. It was developed in 1988 by Pollock. Pollock’s algorithm uses a first order interpolation method that assumes a linear velocity field in each direction. This method offers two main advantages: being analytical, it is fast; and constructing a velocity field that satisfies the flux conservation condition. To improve on Pollock’s algorithm, two approaches are proposed. Either method or a combination of the two can be used in a streamline simulator with very little change to the code. The first method is based on an Adaptive Mesh Refinement (AMR) idea. A local grid refinement is used to recover a more accurate streamline. This method will be particularly well suited to problems with upscaled permeabilities. In such a case, the streamlines are traced using the fine scale permeabilities even though the pressure solve is done on the upscaled permeability field. A higher order method will then be proposed. This algorithm uses the theory of Mixed Finite Elements. Here, it will however be presented in a Finite Difference context. This method will relax Pollock’s assumption of a linear velocity field. Morel-Seytoux’s analytical solution of the quarter five spot pattern injection problem will again be used as a reference solution to test the accuracy of the methods proposed.
4.1. Adaptive tracing refinement
This first method is based on an Adaptive Mesh Refinement idea. We refine the gridblocks through which the streamlines are traced. A step-by-step description of the algorithm is given in the next subsection for two dimensions. The extension to the three dimensional case is straight forward. We note that this method fits perfectly in the adaptive mesh refinement framework proposed by Mallison et al. (2003) to solve the transport problem along streamlines. The 31
subgrids allow more accurate calculations of the time of flight and arc length within the gridblocks. Also, the extra points recovered along the streamline by using a subgrid can be used to create a more accurate grid along the streamline for the transport problem. 4.1.1. Basic concept
The algorithm is here described in a few steps: 1. Refine the gridblock for which the tracing needs to be improved. 2. Allocate fluxes at the boundary of the subgrid 3. Solve for the pressure on the subgrid using the chosen boundary conditions 4. Reconstruct the small scale fluxes at the faces of all cells in the subgrid 5. Trace the streamline on the subgrid using Pollock’s method The present algorithm could be implemented in two ways. First, the grid refinement can be done once and for all streamlines at the beginning of the tracing step. If more than one streamline is traced per block, this implementation will be the most effective. We saw in the Chapter 1 that today’s commercial simulators use a mapping algorithm that requires at least one streamline per block. Using the new Kriging mapping proposed by Mallison et al. (2004) decreases the streamline coverage required to ensure a good mapping quality. In this case, a streamline is not necessarily needed in every gridblock. An implementation with refinement on the fly might then be chosen. It will require less memory and be more computationally efficient than the first implementation. The next two subsections will describe the algorithm step by step. We will start by explaining how to obtain the boundary conditions of the fine scale problem. 4.1.2. Defining boundary conditions for the fine scale problem
This subsection describes step 2 of the tracing method proposed in part 4.1.1. At this point in the algorithm, the regular grid and a subgrid are defined. We also know the Darcy fluxes at the faces of the coarse-scale gridblock. These fluxes were defined such that they satisfy the flux conservation condition at the boundary of the gridblock. We will now see how to define the fine scale fluxes from the Darcy fluxes on the coarse gridblocks. A first possibility is simply to use the large scale Darcy fluxes as boundary condition for the subgrid pressure solve. Figure 4-1 is a cartoon representing the large scale Darcy fluxes and the corresponding straight forward boundary conditions for the subgrid problem.
32
Large scale fluxes
Small scale fluxes
Figure 4-1: Boundary conditions for the subgrid are determined from the large scale fluxes.
Allocating the sub-scale fluxes in such a way provides a refined streamline path. However, the exit point of the streamline traced with these boundary conditions will remain the same as the one obtained by Pollock’s algorithm. This is due to the flux conservation constraint on the boundary of the coarse gridblock that defines uniquely an exit point for a given streamline. So, if the flux distribution over the faces of the coarse gridblock is not changed, the Pollock’s exit point will be conserved. Figure 4-2 sketches this behavior.
Figure 4-2: Streamline paths obtained with Pollock’s method (blue) and the proposed method (red) with constant subgrid boundary conditions.
33
But even with such simple boundary conditions, a more accurate streamline path will be recovered, as shown in Figure 4-2. Using a sub grid allows an improvement of the description of the velocity field in the coarse gridblock. The velocity field will then be approximated by a first order method on each sub gridblock thus improving the tracing accuracy. More importantly, the time of flight and arc length will now be known with a higher order of accuracy. The fine scale boundary conditions can however be improved significantly by reconstructing a piecewise constant flux profile using the subgrid. A different fine scale flux can be assigned to each of the subgrid faces. Subsection 4.2.3 will provide a rigorous proof that a piecewise constant flux profile improves the approximation made on the velocity field. However, we can intuitively understand that using such boundary conditions will provide a smoother velocity field within each coarse gridblock, which is a quality criterion. To reconstruct this piecewise constant flux profile, the fluxes of the four neighboring coarse gridblocks will be taken into account. We will hence use an extended local method to recover a higher order description of the velocity profile at the coarse block faces. The figure 4-3 is providing an example of such a reconstruction.
Figure 4-3: Reconstruction of the subgrid fluxes (black arrow) from the neighboring coarse fluxes (blue arrows)
Figure 4-3 depicts how the small scale fluxes are reconstructed for the face (i-½,j). The first step is to compute the slope σ (i−½ , j ) of the flux across the face according to
σ (i −½ , j ) = minmod(V( i , j +1) − V( i , j ) , V(i , j ) − V(i , j −1) ) ,
34
(4-1)
where V(.,.) designates the coarse face velocity and minmod(a,b) is the slope limiter defined by a , if minmod(a, b) = b , if 0 , if
a
and a.b > 0
a>b
and a.b > 0
,
(4-2)
,
(4-3)
or minmod(a, b) =
1 (sgn( a ) + sgn(b)). min( a . b ) 2
where a and b are real constants. A more complete description of this slope limiter is available in LeVeque (1992). This simple slope limiter ensures that the reconstructed subgrid fluxes take values in a range defined by the coarse fluxes.
V(i-1/2,j+1)
V(i-1/2,j+1)-V(i-1/2,j) V(i-1/2,j) V(i-1/2,j)-V(i-1/2,j-1)
V(i-1/2,j-1) Figure 4-4: Defining the slope of the profile from the neighboring fluxes
In Figure 4-4, the coarse velocities used to determine the slope of the small-scale velocity profile are drawn in blue. σ (i−½ , j ) is the minimum of the two differences shown in red, or is set to zero by the slope limiter if both neighboring fluxes are larger or smaller than the center block flux. Figure 4-5 shows the action of the slope limiter in four different cases.
35
Figure 4-5: Ranges (shown in red) of fine scale velocities allowed by the slope limiter for four cases.
The subscale velocities are determined using
v(ki−½ , j ) = V(i −½ , j ) +
(k − n − 1) σ ( i−½ , j ) , k = 1,...,2n + 1 n 2
,
(4-4)
where n is the level of refinement. n=0 defines the coarse scale. For n ≠ 0, the subgrid is of size (2n+1)x(2n+1). Figure 4-6 shows sub-scale velocities as defined by equation 4-4 for one level of refinement.
36
Coarse profile
v(3i −1 / 2, j )
Sub-grid piecewise constant profile
v(2i −1/ 2, j ) v1(i −1 / 2, j )
Figure 4-6: Reconstructed velocity profile at the face (i-½,j) (red profile)
4.1.3. Solving for the fine scale pressure and velocity field
The boundary conditions of the subgrid problem were defined in the previous subsection. The problem is now closed and the pressure on the subgrid can now be determined with a finite difference method. All the boundary conditions being of Neumann type, our problem is singular and the pressure is only determined up to a constant. We note that the velocity, given by the pressure gradient, is well determined. To find the pressure field, we fix the pressure in a point in the domain. A priori this point can be chosen anywhere and the pressure there can take any value, but we decided to set the value of the fine scale pressure at the center of the subgrid to the value of the coarse scale pressure. This choice forces us to define subgrids with odd numbers of gridblocks in each direction. This choice is however arbitrary. The method will work equally well with grids formed from even number of gridblocks in each direction. The pressure equation can now be solved on the subgrid. As the gridblock on which the streamline is traced does not contain a source term, the governing equation of the flow within the gridblock is a Laplace equation of the pressure, given by
∂2 p ∂2 p kx 2 + k y 2 = 0 . ∂y ∂x
(4-5)
A finite difference approach will be used to discretize this equation and solve for the nodal values of the pressure.
37
From the values of the pressure at the nodes of the subgrid, we can now compute the velocity field on the subgrid using Darcy’s law. Once this fine scale velocity field is reconstructed, Pollock’s method can be used on the subgrid to trace streamlines. 4.1.4. Upscaled permeability field
The method we presented here is especially well suited if the permeability field has been upscaled prior to the simulation, in which case, the subgrid permeability is actually known. The proposed tracing method will fit equally well in a multi-scale adaptive Cartesian framework where not only the subgrid permeability but also the subgrid fluxes are known. We will in this subsection describe the application of the proposed method to the case of an upscaled permeability field. The fine scale permeability field being known, it can be used to trace streamlines on the subgrid. The streamlines will therefore be able to follow the fine scale permeability connectivity undetectable at the coarse scale. For the tracing method to work in such a framework, two main modifications will need to be made to the algorithm. First, the fine scale permeability needs to be included in the sub-scale pressure solve. Then, boundary conditions need to be computed using the fine scale permeability information. To determine the fine scale boundary conditions, a downscaling of the velocity field is necessary. We here propose an extended local downscaling method to compute the fine scale velocity profile from the coarse scale fluxes. The method is based on an upscaling approach proposed by Wen et al. (to be published). The following figure sketches the considered situation. In this example, the fine grid will be taken as being three times as fine as the coarse grid in both the x- and y- directions. The face (i+½,j) common to the gridblocks (i,j) and (i+1,j) is represented in the next figure and will here be used as an example to explain how the small scale fluxes can be recovered.
38
Coarse scale pressure
Fine scale fluxes: qmf Fine Coarse grid
∆x c Figure 4-7: Presentation of the downscaling problem
First, M horizontal layers formed by the subgrid (in our example, M=3) are isolated. The layers, represented by the red blocks in the above figure will be taken of length ∆x c , the coarse discretization length in the x- direction. An average permeability per horizontal layer will now be computed. This average permeability will be given by a weighted harmonic average of the fine scale permeability values. In our case, for example k m* =
1 1 2.k(f2,,(mi ,)j )
+
1 k (f3,,m(i ,) j )
+
1 k (1f ,,m(i)+1, j )
+
1
, m = 1,2,3 ,
(4-6)
2.k(f2,,(mi+)1, j )
where m is the index of the layer, k m* is the average permeability of the layer and k (fn,,(mi ,)j ) is the fine scale permeability of the sub gridblock (n,m) in the coarse block (i,j). Using these average permeabilities per layer as weights, the fine scale fluxes can now be computed through f ,( i + 1 , j ) 2
qm
=
k m* M
∑k
* m
.Q(ci+ 1
2
, j)
, m = 1,..., M
,
(4-7)
m =1
39
where Q(ci + 1
2
, j)
is the coarse scale flux through the face (i+½,j) and M is the number of
layers. The boundary conditions defined here meet important quality criteria. First, one can check that the sum of the fine scale fluxes equals the coarse Darcy’s flux. Second, we note that flux continuity is obtained at the edges of the fine gridblocks. If a local downscaling method is chosen, this property is most likely to be violated. This extended local downscaling method to recover the fine scale fluxes is easily extensible to three dimensional problems. The extended local downscaling approach is the simplest one can find that enforces the necessary conditions on the flux continuity. To improve this downscaling problem, one could, for example, use a local-global method inspired from the method proposed by Chen et al. (2003).
4.2. Higher order tracing algorithm
The tracing algorithm proposed in this section is based on the theory of Mixed Finite Element methods. The algorithm was designed to provide a higher order of accuracy than Pollock’s method. It still relies, however, on the same finite difference pressure solve and is applicable to streamline tracing on two-dimensional Cartesian grids. The extension of the method to trace streamlines in curvilinear quadrilateral and unstructured, trianglebased grids is discussed in Subsection 4.2.6. In Pollock’s method, the assumption that most impacts the accuracy of the tracing is the first order variation of the velocity field. The algorithm proposed relaxes this assumption. The first step is to compute the velocity profiles at the faces of the gridblocks. We will here relax Pollock’s assumption that the velocity at a face is constant in space. We will allow for a linear variation of the velocity profile at the edges of the gridblocks, as shown in Figure 4-8.
40
Figure 4-8: Velocity profiles at gridblock faces in Pollock’s method (left) and the higher order method (right)
The method used to recover such a velocity profile from the finite difference pressure solve will be presented in Subsection 4.2.1. From the velocity profiles at the edges, the velocity field within the gridblock will be interpolated. An important design condition is that this velocity field be divergence free, which reflects the flux-conservative nature of the flow. Subsection 4.2.2 provides the theory on which the interpolation method is based and Subsection 4.2.3 describes how to apply this theory to the interpolation of the velocity field in two-dimensional Cartesian grids. Using the expression of the velocity field, the stream function of the flow will be computed. In Subsection 4.2.4 we present how to recover this stream function. The equation of the streamline is then known analytically as we will see in Subsection 4.2.5. From this equation, the exit point of a streamline in a block can be computed analytically. The higher order tracing algorithm is here summarized in the following steps: 1. Compute the velocity profiles at the faces of the gridblock using information from the neighboring blocks. 2. Define a divergence free velocity field within the gridblock that verifies the boundary conditions. 3. Compute the stream function of the flow. 4. Find the equation of the streamline of interest. 5. Find the exit point from the equation of the streamline.
41
6. Numerically integrate the arc length and time of flight.
4.2.1. Defining the velocity profile on a gridblock face
As in Subsection 4.1.2., we will use the face (i − ½, j ) to describe our method. We compute the slope σ (i−½ , j ) of the velocity profile across this face as
σ (i −½, j ) =
2 minmod(V( i , j +1) − V( i , j ) ,V( i , j ) − V( i , j −1) ) , ∆y
(4-8)
where V(.,.) designates the Darcy velocity through a face and minmod(.,.) is the slope limiter defined in Equations 4-2 and 4-3. Using this slope, we can define the velocity profile across this face as v(i−½ , j ) = V(i−½ , j ) + σ (i−½ , j ) .( y − ( yo +
∆y )) , 2
(4-9)
where v(.,.) is the continuous velocity profile, V(.,.) is the discrete Darcy velocity and yo is the vertical origin of the face (i-½,j). In Figure 4-9, we show a schematic of the usual constant velocity profile (blue) and of the recovered linear profile (red). V(i-1/2,j+1)
v(i-1/2,j+1)
V(i-1/2,j)
v(i-1/2,j)
V(i-1/2,j-1)
Kink due to the slope limiter
v(i-1/2,j-1)
Figure 4-9: Linear velocity profile (red) and Pollock’s constant profile (blue)
42
4.2.2. Velocity spaces
Here, we provide the theoretical background necessary to develop our higher order tracing method. A/ Property of the global velocity field
Mixed Finite Element methods were introduced in Chapter 3. There, we saw that the global velocity field should belong to the space H(div,Ω) where Ω is the domain (the reservoir) on which the velocity is defined, with H (div, Ω) := {u | u ∈ ( L2 (Ω)) 2 ; ∇.(u ) ∈ L2 (Ω)} .
(4-10)
This functional space is a modification of a Sobolev space designed so that the normal trace of the velocity field exists on the boundary of the domain. The integral of the normal trace of the velocity field on the boundary of the domain is none other than the volumetric flux through the boundary. Thus, we understand the importance of being able to construct a well-defined normal trace of the velocity. B/ Defining the global velocity field from local velocity fields on finite elements
Finite Element methods solve the global problem on the domain Ω by using a discretization of the domain into non-overlapping elements Ki, so that m
Ω = ∪ Ki
,
(4-11)
i =1
where m is the total number of elements in Ω. The pressure and the velocity fields will now be discretized on the chosen partition of the domain. For the global velocity field to belong to H(div,Ω), we need to enforce certain properties on the local velocity field defined on each element. To find these properties, one can consider the space Y(Ω), defined by Brezzi and Fortin (1991), given by Y (Ω) = {u | u ∈ L2 (Ω), (∀i ∈ [1, m]), u | Ki ∈ H (div, K i )}
.
(4-12)
A function that belongs to this space must satisfy a global smoothness condition on the whole domain ( u ∈ L2 (Ω) ). Also, the restriction of this function to an element Ki should belong to the space H(div,Ki). The trace of a function of Y(Ω) will hence be defined on any element boundary. One can show (see Brezzi and Fortin (1991)) that a sufficient condition for a function of Y(Ω) to be an admissible velocity field is that the normal traces of the function be continuous at the interfaces of the elements. So we have shown that:
43
If, on each element Ki, the local velocity field belongs to H(div,Ki) and if the normal trace of the velocity field is preserved across each element boundary, then, the global velocity field obtained will belong to H(div, Ω).
In the next two subsections, we show how to construct spaces defining local velocity fields that satisfy these two necessary conditions. C/ Conforming method to approximate H(div, Ω)
A conforming method approximates an infinite-dimensional functional space (e.g. H(div,Ω)) by a finite-dimensional subspace. Therefore, the solution will automatically belong to the larger infinite-dimensional space where the continuum problem is posed. The most common way of defining such spaces is by using polynomial spaces. The dimension of this local space will be defined by the number of degrees of freedom per element. In an element, the degrees of freedom can be defined internally (i.e. source term) or at the edges (i.e. flux). The simplest such space one can define on an element K is called the lowest order Raviart-Thomas space RT0(K) (Raviart and Thomas (1977)) which allows one degree of freedom per face. If we restrict RT0(K) to functions having a zero divergence we obtain RT0o ( K ) = {u | u ∈ RT0 (K ), div(u ) = 0}
.
(4-13)
The velocity field used in Pollock’s algorithm belongs to this space. To improve the streamline tracing, we allow the linear velocity profile at the elements edges defined in 4.2.2. To define these velocity profiles, two degrees of freedom per face are needed. The natural choice for the velocity field in this case is the Brezzi-DouglasMarini space of order 1, BDM1(K) (Brezzi et al. (1986)). Since we want the velocity field to be divergence free, we define BDM 1o ( K ) = {u | u ∈ BDM 1 (K ), div(u ) = 0}
.
(4-14)
The space BDM 1o ( K ) is the space of divergence free velocity fields defined by two degrees of freedom per element face. 4.2.3. Constructing the higher order velocity field
In the previous subsection we found the space defining our higher order velocity field. We will now see how to use this information to build the velocity field in the case of a two dimensional Cartesian grid. We impose the velocity to belong to BDM 1o ( K ) . One can show (Brezzi and Fortin (1991)) that, for a two-dimensional Cartesian element, BDM 1 ( K ) can be given by 44
BDM 1 ( K ) = {u | u = p1 ( x, y ) + r curl( x 2 y ) + s curl( y 2 x), p1 ∈ ( P1 ) 2 , r ∈ ℜ, s ∈ ℜ} , (4-15)
where P1 is the space of polynomials of order smaller or equal to 1. So, v x = a1 + b1.x + c1. y − r.x 2 − 2.s.x. y v y = a2 + b2 .x + c2 . y + 2.r.x. y + s. y 2
,
(4-16)
where a1, b1, c1, a2, b2, c2 are real constants. Enforcing the divergence free condition gives b1 + c2 = 0 ,
(4-17)
so that we obtain a velocity field of BDM 1o ( K ) expressed as v x = a1 + b1.x + c1. y − r.x 2 − 2.s.x. y v y = a2 + b2 .x − b1. y + 2.r.x. y + s. y 2
.
(4-18)
We have a total of seven degrees of freedom to determine to obtain the velocity field. Using seven of the eight boundary conditions available, one can determine the parameters needed. We emphasize here that the eighth boundary condition is redundant as it is fixed by seven boundary conditions and the flux conservation equation. 4.2.4. Computation of the stream function
Our method for tracing streamlines differs from Pollock’s in a fundamental way. Instead of integrating directly the velocity field to obtain the path, we first define the stream function of the flow and then use this stream function to define the streamlines. Let us recall that, since the flow is two dimensional and divergence free, a stream function exists. The relationship between the stream function Ψ(x,y) and the components of the velocity field is given by ∂Ψ ∂y ∂Ψ vy = − ∂x vx =
, (4-19) .
Integrating the velocity field defined in (4-18) using (4-19) gives the stream function of a BDM 1o ( K ) velocity field in the form 1 1 Ψ ( x, y ) = a2 .x − a1. y − b1.x. y + b2 .x 2 − c1. y 2 + r.x 2 . y + s.x. y 2 . 2 2
(4-20)
We note that the stream function for Pollock’s method, for which
45
v x = v x ,o + a x .( x − xo )
, ax =
v y = v y ,o + a y .( y − yo ) , a y =
v x ,∆x − v x ,o ∆x v y , ∆y − v y , o ∆y
,
(4-21) ,
is equal to ΨPollock ( x, y ) = v y ,o .x − v x ,o . y − a x .x. y
.
(4-22)
Here, we could have used ay to replace ax because they are related through the flux conservation conditions ax + a y = 0
.
(4-23)
4.2.5. Tracing the streamlines
Once the stream function is defined, the streamline path is easily found. Using the entry point of the streamline in a block (xi,yi), we can write the equation of the path of the streamline as 1 1 a2 .x − a1. y − b1.x. y + b2 .x 2 − c1. y 2 + r.x 2 . y + s.x. y 2 = Ψ ( xi , yi ) . 2 2
(4-24)
To find the exit point of the streamline, we first find the segments along the faces where the velocity vector points outwards. A streamline can only exit the gridblock through one of these boundary segments. We then look for possible solutions of equation (4-24) on these possible exit segments. These solutions can be found analytically as the stream function reduces to a polynomial of order 2 at the faces of the gridblock. Multiple solutions can be found in two cases. First, if the streamline exit point is at the corner of the gridblock, a solution will be found on each of the face sharing the corner. In this case, the two solutions will in fact give the same location, so there is no problem. Another case may occur where the two possible solutions are not equal. An example is illustrated in Figure 4-10. Here the stream function is plotted in red.
46
Figure 4-10: Stream function (variations of red) that yields two possible exit points
A horizontal plane cutting through the entry point (xi,yi,Ψ(xi,yi)) is also plotted in blue. The intersection of the stream surface and the plane defines, in this case, two streamlines. For this case to happen, the stream function needs to present an extremum in the gridblock. A stagnation point is then present in the domain. In such a case, we could either decide to track another streamline (if the coverage is sufficient in the neighborhood) or we could make a choice between the two possible solutions. To do so,
47
one could use Pollock’s method to determine the lower order exit point and then pick the higher order exit point located on the same side of the extremum. 4.2.6. Extension of the method to two dimensional unstructured grids
The higher order method presented here was developed for two-dimensional Cartesian grids. It is, however, based on the theory of Mixed Finite Elements and therefore applicable to more complex grids. The extension of the method to two-dimensional nonorthogonal structured or unstructured grids is now discussed. The problem of tracing streamlines on an unstructured grid is strongly related to the numerical scheme used to solve for the pressure on this unstructured grid. Several numerical schemes can be chosen to find the pressure field on unstructured grids. The three main options are Finite Volume (FV), Control Volume Finite Element (CVFE) and Mixed Finite Element (MFEM) methods. We identified, in Chapter 3, the flux continuity across the element faces as being a key factor in the accuracy of the velocity field obtained from a scheme. We confirmed this in our study of velocity spaces in Subsection 4.2.2, by seeing that this flux continuity was a necessary condition for the velocity field to be in H(div,Ω). CVFE methods do not enforce the flux continuity at element edges (Cordes and Kinzelbach (1992)). To cope with this problem and recover more accurate streamlines, Cordes and Kinzelbach (1992) developed a tracing method using a refinement of each gridblock into a subgrid formed of triangles. The flux continuity is then enforced on the faces of the subgrid. Using a similar concept, Prévost (2000 and 2001) proposed a new refinement technique based on sub-quadrilateral grid. Both methods were shown to improve significantly the velocity field of the CVFEM and to provide accurate streamlines. The FV and MFE methods, on the other hand, provide the flux continuity across the element edges. The higher order tracing method proposed in this section can be used with any of the three numerical scheme presented here. The difficulty is to define the two fluxes per edge necessary to recover a more accurate flux. If a CVFEM is used, the two necessary fluxes per edge can be reconstructed with Cordes and Kinzelbach (1992) or Prévost (2000 and 2001) methods. In a MFEM, they can be solved for. The FV extension is not yet fully developed. If two fluxes per face are defined, our higher order method can be used. The extension of the method for nonorthogonal structured and unstructured grids is provided below. A/ Nonorthogonal quadrilateral grids
The only change required to be able to apply the proposed method to nonorthogonal quadrilateral grids is the application of a Piola transform (Brezzi and Fortin (1991)) to the
48
velocity field and the transformation of the nonorthogonal quadrilateral into a model square gridblock. The streamline is then traced on this reference gridblock. The inverse transform is then applied to recover the streamline on the nonorthogonal quadrilateral gridblock. B/ Triangular grids
On triangular elements, the velocity components of the Brezzi-Douglas-Marini space of order 1 are given by v x = a1 + b1.x + c1. y v y = a2 + b2 .x + c2 . y
,
(4-25)
where a1, b1, c1, a2, b2, c2 are real constants. The divergence free condition imposes b1 + c2 = 0
.
(4-26)
Enforcing the flux conservation condition removes a degree of freedom in the system and provides the sufficient conditions to ensure the existence of a stream function for the flow. The corresponding stream function is Ψ ( x, y ) = a2 .x − a1. y + c2 .x. y −
b2 2 c1 2 .x + . y 2 2
.
(4-27)
4.3. Results
Pollock’s algorithm and our two methods were implemented in Matlab and streamlines were traced against the analytical streamlines obtained from Morel-Seytoux’s solution to the homogeneous quarter five spot. The errors in time of flight and arc length obtained from the three tracing methods were compared. As expected, the improvements provided by the higher order methods are strongly influenced by the slope limiter used in the reconstruction of the velocity profiles at the faces of the gridblocks. Indeed we saw that if the slopes of the fluxes are set to 0, the exit points of both higher order methods will match the exit point of Pollock’s algorithm. Figure 4-11 illustrates this problem. Streamlines are launched from the first gridblock of the domain. At the faces of this gridblock, the slope limiter prevents the higher order methods to provide a more accurate result than Pollock’s method. In the following figure, the black doted lines are the analytical streamlines; the red streamlines are obtained with Pollock’s algorithm, the blue ones with the subgrid algorithm and the green ones with the higher order method.
49
Figure 4-11: The slope limiter reduces the efficiency of the higher order methods on the homogeneous quarter five spot.
In Figure 4-11, the three tracing algorithms give very similar streamlines. In this case, the slope limiter defines a constant profile for 48.75% of the faces crossed by the streamlines. However, we still recover an improvement on the time of flight calculation as shown in the following table. Table 4-1: Time of flight (TOF) error for the three methods and gain in accuracy when using a higher order method over Pollock’s for unfavorable launching locations. TOF Error % % accuracy gained
Pollock AMR HO 7.71 7.15 5.82 11.26 27.57
If the launching location is now changed to reduce the influence of the slope limiter, the results provided by the higher order methods can be improved as shown in figures 4-12 and 4-13.
50
Figure 4-12: Streamlines obtained for launching locations optimizing the influence of the slope limiter.
Figure 4-13: A close-up on analytical streamlines (black doted), Pollock’s (red), the subgrid streamlines (blue) and the higher order streamlines (green)
51
Figures 4-12 and 4-13 show an improvement in the location of the streamlines when higher order methods are used and launching points are chosen to minimize the influence of the slope limiter. In the case shown above, only 21.25% of the velocity slopes were set to zero. The following table shows the corresponding improvement obtained in the time of flight computation. Table 4-2: Time of flight (TOF) error for the three methods and gain in accuracy when using a higher order method over Pollock’s for favorable launching locations. TOF Error % % accuracy gained
Pollock AMR HO 3.34 1.48 1.38 57.06 62.37
4.4. Conclusions
We have in this part presented two methods for tracing streamlines in two-dimensional Cartesian grids. In the first method, the gridblock is refined. A piece wise constant velocity profile is then reconstructed at the faces of the gridblock using the Darcy fluxes of the neighboring coarse blocks. Thus, the smoothness of the velocity field is increased and a more accurate streamline is therefore recovered. If the method is used in an adaptive framework, the sub-scale permeability field is available and can be used to reconstruct a more accurate fine scale velocity field. The second method proposed is based on the theory of Mixed Finite Element methods. A linear velocity profile at the faces of the gridblock is reconstructed from neighboring gridblocks. The velocity field and stream functions are then calculated, yielding a more accurate streamline and time of flight. We saw that this method was extensible to twodimensional nonorthogonal quadrilateral and triangular based unstructured grids. To test the methods and their efficiencies, the analytical solution of a quarter five spot pattern was used. Significant improvements in the time of flight and arc length computation were obtained. The gain of accuracy in the streamline location was found to be dependent on the launching locations because of the influence of the slope limiter.
52
Chapter 5
5. Optimizing the streamline launching location to minimize the tracing errors. In Chapter 2, the location of the launching point of a streamline was found to be a significant source of tracing errors. Here, we will first explain why the location of the launching point of a streamline as such an influence on the tracing errors. Then, methods to improve the location of the launching points will be proposed. We will then use our quarter five spot test case and show that an optimized launching location significantly decreases the errors in time of flight and arc length.
5.1. Launching from the injectors/producers
In streamline simulation, streamlines are usually launched from the injectors or producers and traced downstream or upstream across the field. This offers the advantage of providing a quick way to identify the area of influence of a given well. Such a property is used for example in injection optimization as described in Thiele and Batycky (2003). To our knowledge, today’s commercial simulators all use the approach used by Batycky (1997) of launching streamlines from the faces of the well blocks. But, as discussed in Chapter 2, launching streamlines in the vicinity of the wells as advocated by the current streamline simulators, may lead to large tracing errors. An illustration of the problem was shown on Figure 2-8 and 2-9. In a reservoir, wells are singularities. Around the wells, the streamlines are either converging (producers) or diverging (injectors). So, the spatial density of the streamlines increases as we get closer to these singular points. As streamlines are diverging when traced away from the wells, errors in the location of the streamline near a well are amplified as the streamline is traced. Therefore, launching streamlines from around the wells usually results in large tracing errors. Conversely, if a streamline is launched away from any singularity, this streamline will converge towards wells, high permeability streaks or fractures. Hence, an error introduced in the streamline location around the launching point will not be amplified as the streamline is traced.
53
To limit the tracing errors, streamlines should be launched as far away as possible from the singularities of the reservoir. In the next part, two methods will be proposed to improve the locations of the launching points.
5.2. Optimized launching location
The following study will be treated in two dimensions for simplicity. However, the methods proposed are completely general and the extension to three dimensions is straight forward. 5.2.1. Launch on a circle a distance away from the singularities
The basic idea behind improving the launching location is to start the streamlines away from the singularities. The simplest algorithm one can think of to do so is to draw a circle around the singularities of the domain and launch streamlines from there.
Figure 5-1: The launching points (red dots) are taken on circles around the wells (blue stars)
Figure 5-1 represents three examples of reservoirs with different well configurations. The blue stars represent the wells. In each case, circles around the wells were drawn and examples of possible launching points, are represented by the red dots. The method is simple to implement but a major draw back is that it is very ad-hoc. Indeed, the choice of the radii of the circles for example, is completely up to the user. 5.2.2. Launch on the Voronoi grid built around the singularities of the domain
A more general way of finding the optimal launching locations is to find the locus of the points that are as far away as possible from all the singularities, which is given by the Voronoi grid around the singularities of the reservoir. Figure 5-3 represents Voronoi grids for the same three well patterns as in Figure 5-2.
54
Figure 5-2: Voronoi grid around three different well configurations.
The Voronoi grid around the wells will be the locus of the launching points that will minimize the tracing errors. The choice of the exact number and location of the launching points is left to the user. However, most streamline simulators today are using a traditional mapping algorithm that requires at least one streamline per block. This condition should be met by launching one streamline per block crossing the Voronoi grid. Indeed, the streamlines launched on the Voronoi grid will all converge to the wells and should, if no major flow barrier or fracture is present in the domain, cross all the gridblocks of the domain.
5.3. Comparison on the quarter five-spot
Here, the two launching methods proposed in the previous section were tested on the quarter five spot and compared to the traditional launching strategy. Figure 5-4 shows the streamlines traced using the three methods. The black circles represent the launching points, the blue streamlines are traced with Pollock’s method and the red ones are the analytical streamlines.
Figure 5-3: Pollock’s streamlines (blue) and analytical streamlines (red) launched from the well block (left), from a circle away from the well (middle) and from the Voronoi grid (right).
55
The following tables show the errors in time of flight and arc length obtained when using three launching options shown in Figure 5-3 for two degrees of refinement.
Table 5-1: Average errors in time of flight for the three launching methods shown in figure 5-3 for a 10x10 grid. Launching location 10x10 grid 100x100 grid
Average time of flight error in % Well block Circle Voronoi 7.6 2.6 2.5 5.3 0.02 0.02
Table 5-2: Average errors in arc length for the three launching methods shown in figure 5-3 for a 100x100 grid Launching location 10x10 grid 100x100 grid
Average arc length error in % Well block Circle Voronoi 10.7 8.4 8.1 2.4 0.8 0.8
We see here that, if the streamlines are traced from the well block, the errors in time of flight and arc length are significantly larger than the errors obtained with improved launching locations. Moreover, a grid refinement has only a very small effect on the accuracy of the time of flight and arc length when streamlines are launched from the well block. On the contrary, if streamlines are launched away from the wells, a convergence can be obtained and errors in time of flight and arc length are quickly decreasing. It was found that improving the launching location will have more effect on the tracing accuracy than a grid refinement. In this chapter, we defined two ways to obtain optimal loci of launching points to reduce tracing errors. However, we saw that obtaining the launching points themselves on these loci is still an open question. In the case of a Voronoi grid, we proposed a heuristics based method to provide at least one streamline per gridblock, as needed by the traditional mapping algorithm. However, to our knowledge, no work has been published so far on finding an optimal number of streamlines to trace. We also need to determine where to launch the streamlines from on the launching loci. In Chapter 6, we will look at the repartition of the mapping errors to decide on these issues.
56
Chapter 6
6. Streamline Coverage The first part of this work was focused on improving the accuracy of the tracing of the streamlines. We will now study a problem that affects significantly the efficiency of compositional streamline simulation, which is a current focus in the SUPRI-C group. In a compositional streamline simulator, most of the computation time is spent in the flash calculations along the streamlines. On the other hand, a minimum density of streamlines is needed to provide an accurate enough solution. Hence, an optimal coverage strategy should be designed to find the optimal number and distribution of streamlines that will provide the desired accuracy and avoid unnecessary computations. Batycky (1997) proposed the mapping algorithm currently used in most streamline simulators. The properties of a gridblock are taken as a time of flight weighted average of the properties picked from the streamlines crossing this gridblock. For this method to work, at least one streamline needs to cross each gridblock of the field. Mallison et al. (2004) provided an improved mapping algorithm that uses a Kriging method to map back and forth the saturations from the streamlines to the finite difference grid. In such a mapping, the properties of a gridblock are determined through an interpolation of the nearest streamline data points. This method was shown to significantly reduce the mapping errors over the classical algorithm developed by Batycky. The work of Mallison et al. offers a second advantage: using Kriging, gridblocks properties can be recovered even for blocks not crossed by streamlines. The coverage required to provide an accurate solution is therefore reduced compared to the traditional method. Using Kriging in mappings essentially changes our vision of streamlines. Streamlines can now be seen as an unstructured flow based simulation grid for transport instead of as carriers of fluid. This novel point of view allows for the development of the concept of partial streamline tracing that will be explained in Section 6.3. To efficiently use the new mapping algorithms and design a streamline coverage scheme, an error indicator is needed to determine the optimal number and launching locations of the streamlines. An error indicator is a variable that provides an a priori measure of the error that will be committed during the mappings.
57
6.1. Error indicator 6.1.1. The Kriging variance - Limitations
A complete description of Kriging methods can be found in the books from Deutsch (2002) or Deutsch and Journel (1998). Kriging is a linear interpolation method used to determine the properties at a gridblock from the properties found on the neighboring streamlines. To find the value of a given point, a Kriging algorithm starts by defining a search neighborhood around the point to be interpolated. The streamlines entering this neighborhood will be considered for the interpolation. The algorithm then assumes a model of continuity by defining a matrix of covariance of the data. The interpolation weights will then be computed based on the continuity model and on the distance from the interpolated point to the data points. Naturally, an error is committed when evaluating the properties of a gridblock using Kriging. Commonly, the Kriging variance is used as an error indicator in traditional Kriging applications. It is defined through
σ K2 = E{[ Z ( x) − Z * ( x)]2 } = Var{Z ( x) − Z * ( x)} ,
(6-1)
where Z ( x) is the true value of a property at the location x and Z * ( x) is the Kriging estimate of this property. This property can in our case be, for example, saturation values, or compositions. Z * ( x) is computed through a weighted average of the neighboring data points, from Nα
Z * ( x) = ∑ λα .Z ( xα ) ,
(6-2)
α =1
where λα are the Kriging weights and Z ( xα ) are the known data points located on the streamlines. Nα is the number of data points falling in the search neighborhood. From equations 6-1 and 6-2 we can compute the Kriging variance from: Nα
σ K2 = Var ( Z ( x)) − ∑ λα .C ( x − xα )
,
(6-3)
α =1
where C ( x − xα ) is the covariance corresponding to the distance between the location x and xα . This covariance is given by the chosen variogram model.
58
The error estimator depends on both the distance to the data points and the chosen model of continuity, but does not include information on the data values. Figures 6-1 and 6-2 present an example that explains the limitations of the Kriging variance as an error estimator for our applications.
Interpolated point
Data points
Search neighborhood
Streamline
Figure 6-1: Limitations of the Kriging variance as an error estimator
Here, two saturation fields are presented, four streamlines are traced (in white) and interpolated from four data points on two neighboring streamlines. For both cases shown in figures 6-1 and 6-2, the Kriging variance will have the same value, because the interpolated points are located at the same distances from the data points on the streamlines. However, in Figure 6-1, the saturation front has not yet reached the interpolated point and the uncertainty in the predicted values is low. The streamline coverage is hence sufficient in this case. In Figure 6-2, on the contrary, the saturation 59
front has reached the search neighborhood. In this case, interpolation errors can be high and an extra streamline could be needed to reduce the interpolation errors. This example motivates the need for an error indicator that includes information on the data to be interpolated and not only on location of the interpolation points. 6.1.2. Global variance - Definition
The most commonly used measure of the variation between data values is the data 2 variance σ data . It is defined as the squared deviation of the values to their mean, or
1
Nα
2 .∑ ( Z ( xα ) − Z ) 2 , = σ data N α − 1 α =1
(6-4)
where Z is the mean of the data values of the neighborhood. We can now define our error estimator as the total variance given by a combination of the Kriging variance, and the data variance. So, 2 2 σ total = σ K2 .σ data
.
(6-5)
We will see in the next section how to use this error estimator to decide whether or not to trace additional streamlines.
6.2. Coverage control algorithm using the error indicator
We suggest a coverage control algorithm separated in two main steps. First, a set of streamlines is launched. To reduce tracing errors, this launching can be done from one of the optimized launching loci proposed in Chapter 5. The compositions are then mapped to the streamlines and moved forward in time along these first streamlines. A total variance map can then be computed. We hence have, at this point, an estimation of the error we would commit by using this coarse streamline coverage. The error map will tell us where additional streamlines are needed and where the coarse coverage suffices. In the second step of the algorithm, we will start tracing a second set of streamlines from the points showing the highest total variance. Once this new set of streamlines is traced, the compositions are mapped on to them and moved forward in time and the total variance map is updated by adding the two sets of streamlines. The second step is repeated as long as some gridblocks still show a total variance above a user-defined threshold.
60
Our algorithm is therefore given by the following steps: A. Initial streamline coverage. 1. Trace a first set of streamlines. 2. Map compositions from the background grid to the streamlines. 3. Move compositions forward in time along the streamlines. 4. Compute the total variance map. B. Improved coverage. If some gridblocks present a total variance above a chosen threshold 1. Trace streamlines from locations where the total variance is above the threshold. 2. Map compositions to the streamlines and move them forward in time. 3. Update the total variance map. C. Repeat Step B until the total variance in the whole domain is below the allowed threshold. D. Map compositions from the streamlines to the background grid.
The desired level of accuracy can be set by the user by choosing the total variance threshold. We note that we must be able to compare the total variances from two different search neighborhoods. Therefore a common scaling of the total variance is needed. A first screening of the field can be used to compute the maximum and minimum total variance obtained. These maximum and minimum values can then be used to scale the total variance between 0 and 1. This will allow the user to choose an appropriate level of accuracy.
6.3. Partial streamline tracing
In the previous section, the provided coverage control algorithm computed a total variance map. So every gridblock in the domain was assigned a total variance value. We will here propose a method that limits the coverage in areas such as well regions or high permeability streaks. In these areas, streamlines generally converge and it is not
61
uncommon to find gridblocks with more than ten streamlines. Such a streamline density is obviously not needed to ensure the quality of the solution. The concept of partial streamline tracing was introduced by Mallison and Gerritsen in the 2004 SUPRI-C Affiliates Meeting and has not, as of now, been published. The idea is to trace streamlines only where needed, stop tracing where the coverage is found to be dense enough. There are two difficulties associated with such a technique. First, boundary conditions for the transport problem have to be found at the starting point of the streamline. Indeed, as the streamlines are not traced all the way to the well, the usual method by which boundary conditions are assigned to a streamline cannot be applied here. The second difficulty to address is the decision of when to stop tracing a streamline. A solution to the first problem was found by Mallison and Gerritsen and will be the subject of a future publication. The error indicator proposed in Subsection 6.1.3. provides a solution to the second problem. Indeed, when a streamline is traced through a block, a simple test on the value of the total variance will determine if the streamline needs to be traced further or not. The methods proposed in this chapter have not yet been fully tested on a streamline simulator. Results of the coverage control algorithm and a detailed description of how partial streamlines can be implemented will be the subject of future articles.
62
Chapter 7
7. Conclusions 7.1. Conclusions
SUPRI-C has been, for the past few years, focusing on improving the accuracy and efficiency of streamline simulation methods. Our part in this group effort was focused on two problems: to improve the accuracy of the streamline tracing method and to provide a coverage control algorithm to improve the efficiency of the method. Streamlines traced with Pollock’s method were shown to have a low order of accuracy in terms of location, time of flight and arc length. We identified three possible causes for these tracing errors and proposed methods to improve on each of them. First, the accuracy of the Darcy fluxes at the faces of the gridblock was identified as a source of error. The use of a Mixed Hybrid Finite Element method was studied. The MHFEM was shown to provide the same accuracy than Finite Difference for a twodimensional Cartesian grid with diagonal permeability tensors. MFEM are however expected to yield more accurate results than Control Volume Finite Element methods or Finite Volume methods on unstructured grids or when used with full-tensor permeability. In the next chapter, two methods were proposed to improve on Pollock’s tracing method. First, a method based on an Adaptive Mesh Refinement idea was developed. A subgrid is used to refine the gridblock. The neighboring gridblocks are then used to recover piecewise constant velocity profiles at the faces of the coarse block. These refined boundary conditions were shown to improve the streamline tracing. The method was also shown to be particularly promising in the case of a simulation run on an upscaled permeability field. Using an extended local downscaling of the velocity field, we were able to trace streamlines on the fine scale permeability field even though the pressure was solved on the coarse scale. The second method proposed was an extension of Pollock’s algorithm to a higher order of accuracy. The method is based on the theory of Mixed Finite Element methods and for that reason was promising in its applications. It was also proven to yield a more accurate streamline part, but is for now, available only for two dimensional cases. The last source of tracing errors we were able to isolate was the location of the launching points of the streamlines. It was shown that streamlines had to be launched as far away as possible from the singularities of the domain to limit the tracing errors. Two options were offered to create optimal loci for launching points: trace circles around the wells or create
63
a Voronoi grid around the singularities of the reservoir. Tracing errors were reduced significantly when streamlines were launched from these locations. The problem of streamline coverage was finally studied. The objective of this part was to improve the efficiency of the streamline method by obtaining a measure of the optimal number of streamlines needed to reach a given level of accuracy in the solution. In streamline simulation, most of the computational time is spent in the flash calculations along the streamlines. Therefore, it is crucial to avoid an overcoverage of regions such as the wells neighborhoods or high-permeability streaks. To determine the necessary streamline density in a region, two parameters need to be determined: first the distance between the streamlines, second the variations of the values from which we will interpolate. Defining an error estimator that takes into account both parameters, we were able to design an algorithm that is expected to determine the streamline density needed as a function of the proximity of the composition fronts. We finally saw how promising the concept of partial streamline tracing is to reduce the coverage of some regions.
7.2. Future work
We will pursue our research work in three directions. First, the subgrid tracing method will be tested on upscaled permeability fields. It is expected to provide an increased solution accuracy in very heterogeneous upscaled permeability fields as the streamlines will be able to capture small scale permeability connectivity undetectable at the coarse scale. The higher order tracing method also requires extensions. We will provide the threedimensional version of the algorithm. In three dimensions, the problem will become more difficult as two independent stream functions are needed to define a streamline. Two stream surfaces will be defined by taking each stream function equal to a constant. The three dimensional streamline will then be the intersection of these two stream surfaces. The theory of Mixed Finite Element methods introduced in this work already provides the velocity field within a three dimensional unstructured gridblock. We however still need to define the dual stream function for this velocity field and solve the algebraic problem of finding the intersection of these two stream surfaces to define the streamline. Finally, we will pursue our work on streamline coverage issues. We will implement the algorithm proposed here in a streamline simulator to report on its efficiency.
64
Nomenclature Miscellaneous xi, xe
= Entry and exit abscissa of a streamline in a gridblock
u
= Darcy velocity
vx, vy, vz
= x-, y- and z- components of the velocity field
vx,o, vx,∆x
= Fluxes through the faces of a gridblock in the x- direction
σ (i−½ , j )
= Slope of the reconstructed velocity profile at the face (i-½,j)
k
= Permeability tensor
kx, ky k
* m
= Permeabilities in the x- and y- directions = Upscaled permeability for the layer m of the subgrid
µ
= Dynamic viscosity
φ
= Porosity
q
= Source or sink term (well flow rate)
Τxi−1 / 2, j
= Transmissibility in the x- direction at face (i-½,j)
A
= Transmissibility matrix
b
= Right Hand Side (RHS) vector
p
= Pressure
S
= Volumetric Saturation of the injected fluid
TPK,j
= Pressure at the face j of element K
p
= Vector of unknowns
Functions δij
= Kroenecker delta
a(.,.)
= Bilinear function
l(.)
= Linear function
Φ(x,y)
= Flow potential
Ψ(x,y)
= Stream function
w(z)
= Complex potential
K(.)
= Complete elliptic integral of the first kind
65
F(.|α)
= Elliptic integral of the first kind of angle α
cn(z,ko)
= Elliptic cosine of argument z and modulus ko
cnx, cny
= Elliptic cosine of argument x- and y- and of modulus ko
snx, sny
= Elliptic sine of argument x- and y- and of modulus ko
ko’
= Complementary modulus of ko
minmode(.,.) = Slope limiter Spaces Ω
= Domain of the simulation
∂Ω
= Boundary of the domain Ω
Ωm
= Element of the partitioning of Ω
K
= Element
2
L (Ω)
= Space of square integrable functions on the domain Ω
H(div,Ω)
= Space defining an acceptable global velocity field
RT0(K)
= Lowest order Raviart-Thomas space on the element K
RT0o ( K )
= Restriction of RT0(K) to the divergence free functions
BDM1(K)
= Brezzi-Douglas-Marini space on the element K
o 1
BDM ( K ) = Restriction of BDM1(K) to the divergence free functions
Geostatistical Notations
66
Z ( x)
= true value of a property Z at the location x
Z * ( x)
= Kriging estimate of Z ( x)
Z (x)
= Average value of Z in the neighborhood of x
λα
= Kriging weights
σ K2
= Kriging variance
2 σ data
= Data variance
2 σ total
= Total variance
C ( h)
= Covariance for a distance h
References Arbogast, T., Wheeler, M.F., Yotov, I., 1997, “Mixed Finite Elements for Elliptic Problems with Tensor Coefficients as Cell-Centered Finite Differences”, SIAM J. Numer. Anal., 34(2), pp. 828-852. Batycky, R.P., Blunt, M.J., Thiele, M.R., 1996, “A 3D Mulit-Phase Streamline Simulator with Gravity and Changing Well Conditions”, proceedings of the 17th International Energy Agency Collaborative Project on Enhanced oil Recovery, Sydney, Australia. Batycky, R.P., 1997, “A Three-Dimensional Two-Phase Field Scale Streamline Simulator”, PhD thesis, Stanford University, CA, USA. Batycky, R.P., Blunt, M.J., Thiele, M.R., 1997, “A streamline-based reservoir simulation of the House Mountain Waterflood”, Stanford Center for Reservoir Forecasting annual report, Stanford University, CA, USA. Brenner, S.C., Scott, L.R., 1994, “The Mathematical Theory of Finite Element Methods”, Springer-Verlag, New York. Brezzi, F., Douglas, J., Marini, L.D., 1985, “Two families of mixed finite elements for second order elliptic problems”, Numerische Mathematik, 47, 217-235. Brezzi, F., Douglas, J., Marini, L.D., 1986, “Recent results on mixed finite element methods for second order elliptic problems”, in Vistas in Applied Math., Numerical Analysis, Atmospheric Sciences, Immunology, Optimization Software Publications, New York. Brezzi, F., Douglas, J., Duran, R., Fortin, M., 1987, “Mixed Fininte Elements for second order elliptic problems in three variables”, Numerische Mathematik, 51, 237-250. Brezzi, F., Fortin, M., 1991, “Mixed and Hybrid Finite Element Methods”, SprinerVerlag New York.
67
Byrd, P.F., Friedman, M.D., 1971, “Handbook of Elliptic Integrals for Engineers and Scientists”, Springer-Verlag New York. Chavent, G., Roberts, J.E., 1991, “A unified physical presentation of mixed, mixedhybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems”, Adv. Water Resources, 14(6), pp. 329-348. Chen, Y., Durlofski, L.J., Gerritsen, M.G., Wen, X.H., 2003, “A coupled local-global upscaling approach for simulating flow in highly heterogeneous formations”, Advances in Water Resources, 26, pp. 1041-1060. Cordes, C., Kinzelbach, W., 1992, “Continuous Groundwater Velocity Fields and Path Lines in Linear, Bilinear and Trilinear Finite Elements”, Water Resources Research, 28(11), pp. 2903-2911. Cordes, C., Kinzelbach, W., 1996, “Comment on “Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity?” by R. Mosé, P. Siegel, P. Ackerer and G. Chavent”, Water Resources Research, 32(6), pp. 1905-1909. Deutsch, C.V., Journel, A.G., 1992, “GSLIB: Geostatistical Software Library and User’s Guide”, Oxford university Press. Deutsch, C.V., 2002, “Geostatistical Reservoir Modeling”, Oxford University Press. Durlofsky, L.J., 1994, “Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities”, Water Resources Research, 30(4), pp. 965-973. Durlofsky, L.J., Behrens, R.A., Jones, R.C., Bernath, A., 1996, “Scale up of Heterogeneous Three Dimensional Reservoir Descriptions”, SPE 30709, SPE Journal, Sept. 1996, pp. 313-326. Emanuel, A., Milliken, W., Chakravarty, A., 2000, “Applications of 3D Streamline Simulation to Assist History-Matching”, SPE 63155, proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, TX, USA. Gander, W. and Gautschi, W., 2000, "Adaptive Quadrature - Revisited", BIT, Vol. 40, pp. 84-101. 68
Goude, D.J., 1990, “Particle Velocity Interpolation in Block-Centered Finite Difference Groundwater Flow Models”, Water Resources Research, 26(5), pp. 925-940. Gross, H., 2003, “Optimization Techniques for History-Matching Using StreamlineBased Geostatistical Constraints”, Master’s thesis, Stanford University, CA, USA. Haegland, Hakon, 2003, “Streamline Tracing on Irregular Grids”, Master’s Thesis, Department of Mathematics, University of Bergen. Hughes, T.J.R., 1987, “The Finite Element Method – Linear Static and Dynamic Finite Element Analysis”, Dover Publications, New York, USA. King, M.J., Datta-Gupta, A., 1998, “Streamline simulation: a current perspective”, IN SITU 22(1), pp. 91-140. Klausen, R.A., Russel, T.F., 2003, “Relationships among some locally conservative discretization methods which handle discontinuous coefficients”, Kluwer Academic Publishers. LeVeque, R.J., 1992, “Numerical methods for Conservation laws”, Birkhauser Verlag. Jessen, K., Orr, F.M., 2004, “Gravity Segregation and Compositional Streamline Simulation”, SPE 89448, proceedings of the SPE/DOE 14th Symposium on Improved Oil Recovery, Tulsa, OK, USA. Mallison, B.T., Gerritsen, M.G., Jessen, K., Orr, F.M.Jr., 2003, “High Order Schemes for Two-Phase Multicomponent Flow”, SPE 79691, proceedings of SPE Reservoir Simulation Symposium, Houston, TX, USA. Mallison, B.T., Gerritsen, M.G., Matringe, S.F., 2004, “Improved Mappings for Streamline-Based Simulation”, SPE 89352, proceedings of SPE/DOE 14th Symposium on Improved Oil recovery, Tulsa, OK, USA. Milne-Thomson, L.M., 1965, “Jacobian Elliptic Functions and Theta Functions” in “Handbook of Mathematical Function”, pp.567-588 edited by M. Abramowitz and I. Stegun, Dover Publications, New York, USA.
69
Morel-Seytoux, 1965, “Analytical-Numerical Method in Waterflooding Predicitons”, SPE Journal, Sept 1965, pp. 247-258. Morel-Seytoux, 1966, “Unit Mobility Ratio Displacement Calculations for Pattern Floods in Homogeneous Medium”, SPE Journal, Sept 1966, pp. 217-227. Mosé, R., Siegel, P., Ackerer, P., Chavent, G., 1994, “Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity?”, Water Resources Research, 3(11), pp. 3001-3012. Mosé, R., Siegel, P., Ackerer, P., Chavent, G., 1996, “Reply”, Water Resources Research, 32(6), pp. 1911-1913. Nedelec, J.C., 1980, “Mixed Finite Elements in IR3”, Numerische Mathematik 35, pp. 315-341. Nedelec, J.C., 1986, “A New Family of Mixed Finite Elements in IR3”, Numerische Mathematik 50, pp. 57-81. Pollock, D.W., 1988, “Semianalytical Computation of Path Lines for FiniteDiffference Models”, Ground Water, 26(6), pp. 743-750 Prévost, M., 2000, “The Streamline Method for Unstructured Grids”, Master’s Thesis, Stanford University, CA, USA. Prévost, M., Edwards, M.G., Blunt, M.J., 2001, “Streamline Tracing on Curvilinear Structured and Unstructured Grids”, SPE 66347, proceedings of 2001 SPE Reservoir Simulation Symposium, Houston, TX, USA. Prévost, M., 2003, “Accurate Coarse Reservoir Modeling Using Unstructured Grids, Flow-Based Upscaling and Streamline Simulation”, PhD Thesis, Stanford University, CA, USA. Seto, C.J, Jessen, K., Orr, F.M.Jr., 2003, “Compositional Streamline Simulation of Field Scale Condensate Vaporization by Gas Injection” SPE 79690, in proceedings of SPE Reservoir Simulation Symposium, Houston, TX, USA.
70
Raviart, P.A., Thomas, J.M., 1977, “A mixed finite element method for second order elliptic problems”, In “Mathematical Aspects of the Finite Element Method”, Springer-Verlag, New York. Tee, G.J., 1997, “Continuous Branches of Inverses of the 12 Jacobi Elliptic Functions for Real Argument”, Dept. of Mathematics, University of Auckland, New Zealand. Thiele, M.R., 1994, “Modeling Multiphase Flow in Heterogeneous Media Using Streamtubes”, PhD thesis, Stanford University, CA, USA. Thiele, M.R., Batycky, R.P., Blunt, M.J., 1997, “A Streamline-Based 3D Field-Scale Compositional Reservoir Simulator”, proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA. Thiele, M.R., 2001, “Streamline simulation”, proceedings of the 6th International Forum on Reservoir Simulation, Sept. 3rd-7th, 2001, Schloss Fuschl, Austria. Thiele, M.R., Batycky, R.P., 2003, “Water Injection Optimization Using a Streamline-Based Worlkflow”, SPE 84080, proceedings of Annual Technical Conference and Exhibition, Denver, CO, USA. Thomas, J.M., 1977, “Sur l’analyse numérique des méthodes d’éléments finis hybrides et mixtes”,Thèse d’Etat, Université Pierre et Marie Curie, Paris, France. Wen, X.W., Durlofsky, L.J., Edwards, M.G., Submitted for publication, “Improved Upscaling of Complex Geological Models on Cartesian Grids”.
71
