WaterSubject

Many codes are available for the modeling of fluid flow and heat transport. A recent overview is given by Anderson. We used the program SHEMAT for the modeling of deep groundwater and heat transport in a sedimentary basin. This program uses a finite volume method (FVM) on an orthogonal grid. Geometrically complex three-dimensional flow domains often cannot be represented well by a rectangular spatial discretization, because the geometry is approximated by steps. Therefore, a more flexible, non-orthogonal discretization was implemented into this well-tested program, preserving as much of the original code’s structure as possible. Only the code which deals directly with the set-up of the system of equations had to be modified.

Typically arbitrary grids are used in codes based on the finite element method (FEM). Here, we wanted to improve an existing FVM code. The code we developed (called sno) is programmed in Fortran, is designed in a modular way, optimized for parallel computing (OpenMP), and uses structured, hexahedral grids. A graphical user interface is also available to assist in setting up input files. The computed numerical results can be stored in different formats for visualization (among others for Tecplot and VTK).

1.2. Computational methods for non-orthogonal grids

Different approaches are available for the computation of fluid flow with the FVM or finite differences method (FDM) on non-orthogonal grids. The integrated finite differences method (IFDM) approach by Narasimhan and Witherspoon may possibly be the earliest approach to use a FVM on non-orthogonal grids in groundwater modeling. It is easy to apply but has the disadvantage that it is based on the criteria that the cell-faces are perpendicular to the connections between the cell-midpoints. Further, the cell-faces should be located approximately in the center. These are significant restrictions for the structure of the grid, in particular with respect to complex geologic structures.

In general, two types of curvilinear grids are distinguished in CFD: grids which are based on curvilinear orthogonal coordinates and grids which are based on general curvilinear coordinates. The first type assumes that the grid lines are perpendicularly arranged in the points of intersection. One of the latter methods has been developed by Demirdži? and Peri?. It is based on a tensorial coordinate transformation. Until now this method has not yet been used for the calculation of flow and heat transport in porous media. It is most appropriate for this purpose – which we will demonstrate in the following.

Lately, Loudyi et al. have presented a code which uses also a FVM on non-orthogonal grids. In contrast to our approach they use an improved least squares gradient technique to take the non-orthogonality into account. Their approach is implemented for groundwater flow in porous media on structured 2D-grids only.

A coordinate transformation method permitting the use of 3D-non-orthogonal grids, has been implemented into an existing code using the finite volume method. The code is applicable to a variety of single phase flow and heat transport problems. As non-orthogonal structured grids allow a better spatial representation of the modeling domain, eventually with a smaller number of nodes, the new code may be helpful and can increase the correctness of results and of predictions deduced from them. Additionally it is possible – often without increasing the number of nodes – to refine areas with higher groundwater velocities for instance along fault zones. However, as the code is restricted to structured hexahedral grids, it is not possible to apply adaptive mesh refinements as in some finite element codes. The coordinate transformation method applied here is straightforward, increasing the flexibility of finite difference and finite volume modeling for complex geometries considerably. For geothermal basin analysis it is advantageous to model the basin geometry as accurately as possible. The applicability of 3D body-fitting grids is favorable for this purpose, of course.

Sufficiently accurate results have been obtained from analytical and benchmark problems. In contrast to the computational error which is discussed in this paper, Fletcher gives some information with regard to the mathematical errors associated with the use of generalized coordinates.

In general, appropriate mesh quality criteria are important when non-orthogonal grids are used. Numerical errors can occur in case of highly skew grids. This is due mainly to the dominance of the explicitly solved cross-derivatives.

A restriction for the presented approach is that buoyancy dependent models can be only computed correctly, if the vertical grid lines are parallel to the gravity vector. Although most hydrogeologic subsurface models fulfil this requirement, in some rare cases more flexible grids may be necessary. Circumventing this shortcoming will be work for the future. As discussed earlier two solutions could be possible: a different formulation of the basic equation or a least-squares method.

A simple graphical user interface is available which uses GOCAD and Tecplot/VTK for pre- and post-processing, respectively. It is coded in Java and therefore applicable at most computer systems.

In the near future, the code related to the coordinate transformation will be released under the LGPL open-source licence. The main-numerical routines are closed source and therefore unavailable. However, executables of the code can be obtained upon request from the first author.