|
|
MODELLING OF FLOW IN FRACTURED AQUIFERS USING AUTOMATIC GRID ADAPTATION
Institut
für Strömungsmechanik und Elektronisches Rechnen im
Bauwesen
Universität Hannover
Appelstrasse 9A, 30167 Hannover, Germany
email: kaiser@hydromech.uni-hannover.de
ABSTRACT
This paper deals with the utilization of automatic grid adaptation
(h-adaptive method) on fluid flow processes in fracture networks. The
h-adaptive method is used to improve the numerical solutionfor fractured
aquifer simulations. As an application, a fault system near
Soultz-sous-For\^ets in France is investigated, where a geothermal research
site (Hot Dry Rock System) is under development.
The finite element program Rockflow-3 has been developed to simulate
complex flow and transport processes of one and more fluid phases in subsurface
hydrosystems, combining coupled modules with the method of automatic grid
adaptation (Kolditz et al. 1998, Kaiser et al. 1998, Habbar et al. 1998,Thorenz
1998).
The different flow zones in the fractured aquifer, such as surroundings
of sinks and sources, the interaction area between fracture and matrix, and the
area between zones of different material parameters are discretized for
3d-problems by an automatically adapted grid. The velocity field can now be
modelled with high accuracy for each point of time, and the new grid can be
used later for an adaptive tracer and/or heat transport simulation.
The numerical solution of the flow equation is based on a semi-discrete
finite element method for unstructured grids where the porous rock matrix is
discretized by 3d-volumes, fractures by 2d-planes and main pathways within the
fractures by 1d-line-elements.
Keywords: fracture network, network meshing, automatic grid adaptation
1 INDROTUCTION
In recent years, fractured subsurface systems have received increasing
attention in several geoscientific and geotechnical disciplines dealing with
water supply from karst aquifers, waste deposition in the geologic subsurface,
as well as in the context of utilization of geothermal energy (hydrothermal or
hot dry rock reservoirs). Near Soultz-sous-For\^ets (France) in Northern Alsace
a European HDR (Hot Dry Rock) project has been under progress since 1987. An
artificial underground heat exchanger has been created in the granitic basement
of the Rhine Valley (Pribnow et al. 1998).
In a first step, fluid flow processes in natural fault zones were
computed using automatic grid adaptation. The fault planes are represented by a
2.5d-fracture network, i.e. intersected planes in a 3d-space (Figure 1,
right). A numerical model has been
developed to compute coupled fluid flow and heat transport processes. In this
case the 3d-start meshes have to be constructed by conformly coupled
1d,2d,3d-elements, where the 3d-elements represent the rock matrix (Figure 1,
left). In order to mesh the fracture network, unstructured grids are required.
Automatic grid adaptation controlled either by heuristical (differences,
gradients, curvatures) or by analytical (Babu\v{s}ka) indicators is suitable
for problem-related discretization. Refinement is limited to zones where the
finite element solution, based on a coarse grid, cannot approximate the exact
solution satisfactorily. As a result, coarse meshes are allowed in regions with
small changes of the computed field variable, and meshes are refined in areas
containing large errors.
For the calculation of the velocity field an incompressible flow is
assumed. A stationary velocity field will usually be required for long-time
heat transport simulations.
Figure 1: Conform coupled elements and a small system of intersecting
fractures (2d-elements)
The unstructured start mesh needs contain only the basic structural
information on the reservoir geology. Nevertheless, the generation of this mesh
is non-trivial, especially for 3d-problems. During the simulation process, the
mesh is automatically adjusted where necessary. The aim of the proposed
automatic grid adaptation method is to achieve a more effective simulation of
strongly heterogeneous reservoirs.
2 GEOMETRICAL MODEL
Geometrical
description of fracture systems:
Fracture networks are modelled by using planes in the 3d-space. Each
plane can be described by its center and two angles (strike and dip) to define
the location and the orientation of the fracture. In addition, a finite
diameter of the planes must be defined. The planes can be modelled by using a
CAD system as well. In this case each fracture is defined by a polygon
describing the plane boundary. The planes can be oriented arbitrarily in the
3d-space.
Generating single fractures (e.g. AutoCAD):
The CAD system AutoCAD creates a DXF file which has to be converted to
the internal format of the geometrical mesh generator RF-HGM/Meshview (Kasper
1998). For the geometrical description of the planes, the data of the drawing
elements LINE, POLYLINE and 3DPOLYLINE are converted. In addition, topological
information can be given by using different LAYERs (AutoCAD).
Meshing of fracture systems:
After creating and converting the boundaries of the fracture planes,
each plane has to be triangulated (Kasper et al. 1998). Some of the resulting
Triangulated Irregular Networks (TIN) penetrate each other. Each fracture is
treated as a single TIN. In order to create a fracture network the lines of
intersection between the different TINs have to be calculated. At the beginning
two TINs are coupled by appending the triangles of the second TIN to the first
one. Two different results can be expected. In the first case the triangles of
the TINs do not penetrate each other. The result is a first network (also a
TIN) with two independent fractures. For penetrating triangles, the lines of
intersection are calculated. Each triangle involved in the intersection has to
be modified. The resulting network (TIN) can contain two up to four fractures
(e.g. three fractures, if one TIN is divided completely by the other) separated
by the intersection lines. In a following step, the other TINs are coupled with
the existing network as described above.
The mesh
quality has to be improved in order to get a grid for numerical analysis. Several
functions for inserting or removing nodes on the intersection lines or inside
the fractures are available. The resulting mesh can be smoothed using the
Laplacian method as well. Finally a conversion from triangles into
quadrilaterals can be realized.
3 MATHEMATICAL MODEL
With the
assumption of an imcompressible flow the flow equation becomes linear. The
partial differential equation of groundwater flow is then given by
(1)
including extented Darcy's law (non-linear flow behaviour)
(2)
with S: storage coefficient, t: time; 
:
fluid source term, r: density, p: pressure, q: seepage velocity tensor, m: viscosity and k: permeability tensor.
The
parabolic flow equation is treated with the Bubnov-Galerkin method in space and
the implicit Euler method in time.
The flow
equation and the heat transport equation wich reads as
(3)
(with T: temperature, cr: heat capacity of the porous
medium, 
:
specific heat capacity of water, 
:
density of water, D:
dispersion-diffusion tensor of the porous medium and 
:
heat source term of the porous medium)
are connected by the velocity field. A
stationary velocity field will normally be required for the computation of
long-time transport processes , so that a storage coefficient equal to zero can
be selected. The adaptive method equally works for unsteady fluid flow
problems. In this case the Neumann criterion must be considered.
The
implemented h-adaptive algorithm (Schulze-Ruhfus 1996) for mixed element types
(Figure 2) requires irregular nodes in order to establish new elements with a
shape identical to that of the original elements. These irregular nodes are
eliminated on the element level.
The advantage of this procedure lies in its strongly hierarchical
character. In this way it is possible to reverse the refinement of any element
if necessary, e.g. in the case of tracer or heat transport (Barlag 1997).
Critical areas can be found by heuristical refinement indicators, e.g.
differences, gradients or curvatures. For each element type, different
indicators can be applied. In the following application the curvature (jump)
indicator
(4)
is used for the 2d-elements (fractures), where U denotes the numerical solution and a the area of the element E.
The same indicator can also be used for 3d-elements. This indicator calculates
(because of the linear shape functions) the maximum of the jumps across the
element sides t in the normal component of the gradient of U.
The geometry of the reservoir model is based on observations in the
Soultz-sous-For\^ets HDR site, France. After successive intersection of single
fractures, the resulting triangular mesh (Figure 3, left) has been refined,
smoothed and converted into a quadrilateral mesh (Figure 3, right). This
2.5d-fracture network represents the fault planes observed at the Western flank
of the Rhine Valley (Pribnow et al. 1998) and has been used to compute fluid
flow processes in these natural fault zones. The fractures are assumed to be of
constant thickness, and they are penetrated by two boreholes (sink and source)
with constant discharge and recharge rate.These boreholes can be located in
Figure 4 where the automatic grid adaptation had occurred.
Figure 3: Triangular mesh and quadrilateral initial grid
The adaptive grid controlled by the curvature indicator is refined in
areas of large changes of pressure gradients, in particular near the sinks and
sources. Figure 4 shows an automatically refined fracture network (left) and
the velocity
field in one of the fractures (right).
Figure 4: Automatically refined grid
This work was partly supported by Geowissenschaftliche
Gemeinschaftsaufgaben (GGA). The authors appreciate the support of D. Pribnow
and C. Clauser, in particular their help in data acquisition.
REFERENCES
Barlag,
C., Zielke, W.,
Dynamic Adaptive Transport Simulation in Fractured
Subsurface,
In: Proc. XI Int. Conf. on Computational Methods in
Water Resources, Cancun, Mexico, July 22-26, 1996.
Barlag,
C.,
Adaptive
Methoden zur Modellierung von Stofftransport im Kluftgestein,
Dissertation,
Institut für Strömungsmechanik und Elektronisches
Rechnen
im Bauwesen, Universität Hannover, 1997.
Habbar,
A., Kolditz, O.,
Numerical modelling of reactive transport processes in
fractured porous media,
In: Proc. Int. Conf. on Groundwater Quality GQ'98, pp
313-319, Tübingen, Germany,
September
21-25, 1998.
Kaiser,
R., Kolditz, O., Zielke, W.,
Automatic grid adaptation for subsurface fluid flow
problems - Application
to fractured-porous reservoirs,
In: Proc. XII Int. Conf. on Computational Methods in
Water Resources, pp 125-132, Crete, Greece, June 15-19,
1998.
Kasper,
H.,
Ein 3-D
CAD-System zur numerischen Grundwassermodellierung,
Dissertation,
Institut für Strömungsmechanik und Elektronisches
Rechnen
im Bauwesen, Universität Hannover, in progress.
Kasper, H., Kosakowski, G., Rother, T., Thorenz, C.,
Kolditz, O. & Taniguchi, T.,
Development of a 3-D CAD System for Numerical Analysis
of Subsurface Flow \& Transport,
In: Proc. 6th Int. Conf. on Numerical Grid Generation in
Computational Field
Simulation,
pp 683-694, Greenwich, UK, July 06-09, 1998.
Kolditz,
O., Habbar, A., Kaiser, R., Kasper, H., Schulze-Ruhfus, M.,
Rother, T., Thorenz, C., Zielke, W.,
Software Concept of Simulating Coupled Processes in
Subsurface Hydrosystems,
In: Proc. Hydroinformatics 98, pp 613-618, Copenhagen,
August 24-26, 1998.
Pribnow, D., Clauser, C.,
Heat- and Fluid-Flow in the Rhine Graben: Regional and
Local Models for
a Hot-Dry-Rock System, 4th HDR Forum, Strasbourg,
France, September 28, 1998.
Schulze--Ruhfus,
M.,
Adaptive
Verfeinerung und Vergröberung gekoppelter 1D/2D/3D--Elemente,
Diplomarbeit,
Institut für Strömungsmechanik und Elektronisches Rechnen im Bauwesen,
Universität Hannover, 1996.
Thorenz, C.,
Numerical modelling of multiphase flow and tracer
transport,
In: Rockflow-Manual version 3.3, 1998. \\
http://www.hydromech.uni-hannover.de/w3-pages-rockflow/Rockflow.html