GTH: A Finite-Difference Simulator for Two-Phase Hydrothermal Fluid Flow in Porous Media

Introduction

Mathematical Model

Fluid Properties

Finite-Difference Model

Numerical Solution

Program Structure

Programming

Example Problems

Appendices

References

GTH is a numerical simulator developed at Georgia Tech to study liquid/vapor two-phase flow and heat transport in porous media, either uniform or fractured. The two-phase hydrothermal systems simulated by GTH may have non-uniform and time-dependent porous medium properties, such as porosity and permeability. Water properties are calculated using Haar et al.'s [1984] equation of state. GTH is capable of handling problems with pressures from 0.05 MPa to 1,000 MPa and temperatures from 0 C to 1200 C.

The mathematical model of GTH is based on Faust and Mercer's [1979] description of two-phase hydrothermal systems. The major assumptions employed by GTH are: (1) potential and kinetic energy are negligible; (2) thermal equilibrium exists among water, steam and porous media; (3) Darcy's law is applicable to individual phases in two-phase hydrothermal systems; and (4) capillary effects other than that on liquid or vapor residual saturations are negligible.

GTH was developed based on a control-volume-based finite difference method [Patankar, 1980]. Because of its conservative nature this method is suitable for numerical modeling of systems described by a differential equation set that involves at least one first-order hyperbolic equation. One such example is the equation set posed for the two-phase hydrothermal systems [e.g., Young, 1993; Xu, 1996].

This document briefly describes the mathematical model and numerical methods applied in developing GTH, and how to install, compile and run the program.

Symbols used in this document are defined in Appendix A.

Two-phase (liquid-vapor) fluid flow and heat transport in a porous medium can be described by a set of conservation equations [Xu and Lowell, 1998]. Conservation of mass and energy are expressed by

,(1)

and

,(2)

respectively. Conservation of momentum for each phase is accounted for by Darcy's Law. Thus:

(3)

and

.(4)

Bulk density r and bulk enthalpy h in equations (1) and (2) are defined as:

,(5)

and

,(6)

respectively. The volume saturation of the individual phases sums to unity. That is:

.(7)

Pressure and enthalpy should be specified for the whole system that is to be simulated. For the pressure equation the boundary condition may be specified as fixed pressure or mass flux or their combination. For the enthalpy equation it may be fixed enthalpy or heat flux or their combination.

In many circumstances a hydrothermal system is subject to certain type of injection or extraction operation. In-situ heat or mass sources may also present. These are represented by source/sink terms in equations (1) and (2). GTH is able to simulate the effect of time-dependent sources and sinks on a two-phase hydrothermal system.

GTH is designed to deal with problems with non-uniform porous medium properties such as porosity, permeability and thermal conductivity. The permeability handled by GTH can also be anisotropic, that is, at a given place permeability has non-identical values at different directions. Permeability may vary with time during the evolution of a hydrothermal system. Certain type of pre-described time-dependent permeability variation can be handled by GTH.

Thermodynamic properties are calculated at given pressure and enthalpy. Haar et al.'s [1984] equation of state, which is valid within a temperature range of 0-1200° C and a pressure range of 0.1-1,000 MPa, is used to calculate pure water properties within sub-cooled liquid, super-heated steam and super-critical regions. The equation of state is given by the Helmholtz free energy as a function of temperature and density

.(8)

Temperature, density, and their derivatives with respect to pressure and enthalpy were calculated as functions of pressure and enthalpy and were then organized to form a data file named TABLES.DATA. During each running of GTH this data file is read first. For a given pair of pressure and enthalpy, the property data at the nearest pressure-enthalpy points are first searched and the properties corresponding to the pressure-enthalpy pair are then calculated through a bilinear interpolation of the data at the nearest points [Press et al., 1986].

The calculation of dynamic viscosity of single-phase pure water is based on Kestin et al.'s [1984] formulation. The viscosity is expressed as a function of temperature and density

.(9)

The formulation is valid for the following range of pressure and temperature:

p £ 500 MPafor150 ° C ³ T ³ 0 ° C,p £ 350 MPafor600 ° C ³ T ³ 150 ° C,(10)p £ 300 MPafor900 ° C ³ T ³ 600 ° C.

The properties of coexisting phases include the saturated density, enthalpy, temperature and viscosity of liquid and vapor. Among them temperature is the same for the two individual phases assuming thermodynamic equilibrium. The temperature data were fitted into a curve represented by a fourth order semi-logarithmic polynomial of pressure. The temperature and its derivative with respect to pressure are calculated along this curve. The other phase boundary data were treated as functions of pressure and made into lookup tables and put into a data file named PHASEBD.DATA. A searching and linear interpolation procedure [Press et al., 1986] is used for calculating them from the lookup tables as functions of pressure.

Pressure and enthalpy are used by GTH as dependent variables since they uniquely define the thermodynamic state of pure water under both single- and two-phase conditions. Within the two-phase region, pressure and temperature are related to each other according to the Clapeyron curve . In the two-phase region, therefore, the independent variables p and T are interchangeable, as are the other set of saturation S, enthalpy h, and density r .

Before the discretization of the governing equations, Equation (1) is re-written, by utilizing equations (3) and (4) and taking derivatives of density with respect to pressure and enthalpy, as

(11)

Equation (2) is similarly re-written as

(12)

so that it becomes explicitly enthalpy-related. The two governing equations (11) and (12) are nonlinearly coupled. They are then discretized to produce two difference equations for pressure and enthalpy respectively using an approach that is similar to the method of Patankar [1980]. The finite difference equations thus obtained have the following forms:

(13)

(14)

In obtaining (13) and (14) from (11) and (12) the following two aspects need to be taken care of.

For the sake of numerical stability, the discretization of the governing equations in time is fully implicit. That is, the discretized equations are solved for pressure and enthalpy iteratively. All coefficients of the discretized pressure and enthalpy equations are calculated according to the newly-obtained pressure and enthalpy after each iteration. This method has been commonly used in the numerical simulation of fluid flow and heat transfer problems [e.g., Patankar, 1980].

The pressure equation is of diffusion type and was discretized by the well-known central-difference method. However, the enthalpy equation is an advection-diffusion equation in the single-phase region and a first-order hyperbolic advection equation in the two-phase region. It needs some special treatment especially for the advection-related term. GTH applies Patankar's [1980] method for all single-phase regions.

Within the two-phase region, the enthalpy equation reduces to a first-order hyperbolic equation because the heat conduction term becomes a function of pressure only and is independent of enthalpy. The calculation of the coefficient of the advection-related term at the interfaces between the adjacent control volumes is crucial in terms of numerical stability. There are many options in calculating the coefficient.

An effective way to approach numerical stability, which is adopted by GTH as one of the options, is to apply the so-called full-upstream weighting to the relative permeabilities [Aziz and Settari, 1979]. That is, the relative permeabilities at an interface between two adjacent control volumes are replaced by those of the control volume that is at the upstream-side of the interface as determined by the flow direction of the individual phases. This method has proved to be quite efficient in many circumstances and in most circumstances it is also sufficiently accurate.

There are several other options offered by GTH. The simplest one is central differencing. One may want to choose this option when the propagation of saturation disturbance is most concerned [Xu and Lowell, 1998]. However it is good only for fine grids and small time steps. A second one was obtained by formally applying Patankar's [1980] method to the two-phase region. While the approach is sometimes much better than the central differencing, it has been observed in some circumstances to suffer from numerical instability. A third way is to calculate the interface properties according to the pressure and enthalpy that are obtained by assuming a local steady-state one-dimensional configuration determined by the pressure and enthalpy of the two adjacent control volumes. This method has proved to be very accurate in obtaining the final steady states for one-dimensional problems. However, it may cause errors near discontinuous interfaces and when applied to strongly transient problems.

In order to save computer space for multi-dimension problems, the two finite difference equations (13) and (14) are solved using a method that is commonly applied to one-dimensional problems [Patankar, 1980]. For a two-dimensional problem, for instance, the finite difference equations are first solved line-by-line in the X-direction using the Thomas algorithm. The same equations are then solved row-by-row in the Y-direction. This procedure is repeated several times or until the solutions are converged. When solving the equations the updated pressure and enthalpy values, whenever and wherever they are available, are utilized. The procedure may be repeated for several times during each iteration between the pressure and enthalpy equations. It is not necessary for the solutions to be convergent because they are the results of an intermediate step during a pressure-enthalpy iteration procedure.

The pressure and enthalpy equations are nonlinearly coupled. An iteration procedure was designed to obtain pressure and enthalpy solutions for each time step by solving (13) and (14) for pressure and enthalpy respectively. The pressure equation (13) is solved first for pressure and then the enthalpy equation (14) for enthalpy, and thus complete one solution of the iteration procedure. The coefficients of (13) and (14) are then calculated using the new pressure and enthalpy values and the equations are then solved for pressure and enthalpy again. Before solving the enthalpy equation (14) one may choose to or not to calculate coefficients b using new pressure values. The iteration procedure is repeated till it converges. A relaxation procedure is designed to help the convergence of the above iteration. The relaxation procedure is used by GTH to reduce the rate of pressure and enthalpy changes when strong variations occur and to accelerate the rate when it is slow during the iteration. The change rate is also slowed down by the relaxation procedure when phase change is taking place within some control volumes.

The convergence criteria of the iteration procedure for the coupled pressure and enthalpy equations are given by

(15)

and

,(16)

where p and h are the solutions obtained during the current pressure-enthalpy iteration, p₀ and h₀ are those obtained during the last iteration, p⁰ and h⁰ are the values obtained at the last time step, and E_p and E_h are sufficiently small numbers. Usually E_p and E_h are set to be equal or less than 10^-2. The iteration procedure is said to be convergent when both (15) and (16) are satisfied. Sometimes when phase change occurs within one or a few control volumes, the status of the control volume(s) tends to switch forth and back between two different phases during an iteration and the convergence is difficult to reach. An averaging algorithm is designed to help the convergence of the iteration procedure. Basically, when pressure and enthalpy variations during the switching are sufficiently small, averaged values of the pressure and enthalpy over the last several iterations, instead of the numbers obtained by the last iteration, are chosen as the values of the current iteration. This averaging approach has been proved to work well in improving the iteration convergence in most circumstances.

A maximum time-step interval is pre-set and may be changed during a GTH run. When the pressure-enthalpy iteration converged at the last time step and the time-step interval is smaller then the maximum interval the time-step interval for the next time step is automatically increased. If the pressure-enthalpy iteration does not converge then the time-step interval will be cut to a smaller value and redo the iteration.

Program Structure

The main program is named GTH. A typical GTH run starts with reading tables of pre-calculated thermodynamic properties (tables.data) and phase boundary data (phasebd.data). GTH then read parameters of the physical problem, numerical model, initial conditions and program control parameters from gxin. Boundary conditions and source/sink information are read in separately from bc.data and bgsosi.data. When data input work are finished, GTH proceeds to solve the problem at each time step using a DO-loop.

For each time step, an input file named answers.data is read first. The data file contains parameters that control the program process and can be changed by the user during the programming to improve the program process. Finite difference pressure and enthalpy equations are solved using an iteration procedure that is represented by the pressure-enthalpy iteration DO-loop.

For each pressure-enthalpy solution a subroutine named PROPTY is called first to calculate thermodynamic and thermophysical properties of water at given pressure and enthalpy values and to insure phase equilibrium if liquid and vapor coexist in the problem domain. The subroutine INTERFACE is then called to calculate properties for each interface between adjacent control volumes. Flow rates are calculated applying Darcy’s law and various types of upstream weighting, if they are chosen, may be conducted for advection-related terms such as relative permeabilities. Results of above calculations, combined with pre-introduced model parameters, are used to form finite difference equations for pressure and enthalpy by calling EQTNP and EQTNH, respectively. Finally SOLVEQ is called to solve the equations for pressure and enthalpy, respectively. A relaxation procedure is employed to improve the convergence of the pressure-enthalpy iteration.

At the end of a successful pressure-enthalpy iteration the results are stored temporarily in G3pickup in high precision formats for restarting and continuing the program in cases of unexpected program termination as the result of a power failure, a computer down, or the difficulty in program continuation due to convergence problems. At times preset by the user results are also written into output files G3out1, G3out2, G3out3 and sometimes also G3out4.

GTH stops after it successfully completed the works for all time steps.

• BC – read boundary conditions and boundary indices that describe the location and type of boundary conditions from bc.data. It is called from DATA.

• BGSOSI – read from bgsosi.data the source/sink information for mass and energy within the simulation domain and parameters that control the permeability variation with time if this is the case. It is called from DATA.

• COMMON.GEOTH – store arrays and parameters that used by GTH, BC, BGSOSI, DATA, DENDRIV, DIF, EQTNH, EQTNP, GFILES, INITIAL, INTERFACE, INTERPHASE, INTFPROP, LBYL, MYANSWERS, PDATA, PHASEINDEX, PICKUP, PRESGRAD, PRINTOPT, PROPTY, RD, READ, RELPERM, RESET, SOLVEQ, SOLVLINE, TIMEPARA, TTOH, UPSTREAM, VARIATION, VELOCITY, and WR.

• COMMON.PHASEBD – store arrays for phase boundary data. It is included in GTH, D2PHASE, DENDRIV, H2PHASE, INTERFACE, INTERPHASE, INTFPARA, PHASEINDEX, PROPTY, READTABLE, and VISC2PHASE.

• COMMON.TABLET – store arrays for pre-calculated thermodynamic data. It is included in GTH, INTERFACE, INTERPHASE, INTFPARA, PROPTY, and READTABLE.

• D2PHASE – calculate liquid and vapor density values at a given pressure. It is called from PROPTY, DENDRIV, INTFPARA, and UPSTREAM.

• DATA -- read parameters of the physical problem, numerical model, initial/boundary conditions and program control parameters from gxin, bc.data and bgsosi.data. It is called from GTH.

• DENDRIV – calculate density derivatives with respect to pressure and enthalpy for two-phase mixture. It is called from PROPTY.

• DENPT – calculate enthalpy and density at given pressure and temperature for single-phase pure water. It is called from TTOH.

• DIF – calculate the difference of pressure and enthalpy values at adjacent iteration levels for relaxation of the pressure-enthalpy iteration. It is called from GTH.

• EQTNH – form finite difference enthalpy equations based on the energy conservation equation. It is called from GTH.

• EQTNP – form finite difference pressure equations based on the mass conservation equation. It is called from GTH.

• GEOTH.FTN – convert time unit either from second to year or from year to second. It is included in GTH, BGSOSI, DATA, PDATA, PICKUP, PRINTOPT, and TIMEPARA.

• GEOTH.PAR – specify the maximum numbers of control volumes in X-, Y- and Z-dimensions. It is included in COMMON.GEOTH.

• GFILES – open input and output files for reading and writing. It is called from DATA and PICKUP.

• GTH – is the main body of the program.

• H2PHASE – calculate liquid and vapor enthalpy values at a given pressure. It is called from PHASEINDEX, INTERPHASE, INTFPARA, and UPSTREAM.

• HELMHOLTZ – calculate Helmholtz free energy and its derivatives as functions of temperature and density using the equation of state for pure water [Haar et al., 1984]. It is called from DENPT, PROPTY and INTFPARA.

• INITIAL – refresh variables and properties using the values obtained at the last time step. They are used as the initial values for the current time step. It is called from GTH.

• INTERFACE – calculate variables and properties at each interface between adjacent control volumes. It is called from GTH.

• INTERPHASE – determine calculate variables and properties at each interface between adjacent control volumes using analytic solutions for 1-D steady states. It is called from INTERFACE.

• INTFPARA – calculate properties and parameters at each interface between adjacent control volumes as a function of given pressure and enthalpy. It is called from INTERFACE and INTERPHASE.

• INTFPROP – determine porous medium properties for at each interface between adjacent control volumes. It is called from GTH.

• LBYL – solve finite different equations line by line, row by row, and column by column. It is called from SOLVEQ.

• LOCATE – search property data tables. It is called from PROPTY, D2PHASE, H2PHASE, and VISC2PHASE.

• MYANSWERS – read data contained in answers.data. It is called from GTH and RESET.

• PDATA – write simulation results into output files. It is called from GTH.

• PHASEINDEX – determine phase status for each control volumes. It is called from GTH and PROPTY.

• PHDERIV – calculate pressure, temperature and their derivatives following Haar et al. [1984]. It is called from DENPT, PROPTY and INTFPARA.

• PHREG – assign a phase index the specific control volume: -1 for liquid, 1 for vapor, and 0 for two-phase mixture. It is called from PHASEINDEX and INTFPARA.

• PICKUP – write simulation results for each time step into G3pickup. It is called from GTH.

• PRESGRAD – calculate pressure gradient for each interface between adjacent control volumes. It is called from VELOCITY.

• PRINTOPT – read parameters that control what data are to be written into output files. It is called from DATA.

• PROPTY – calculate thermodynamic and thermophysical properties, determine phase equilibrium and relative permeability, as a function of pressure and temperature for each control volume. It is called from GTH.

• RD – read arrays of input data. It is called from READ.

• READ – read parameters for the geometry model and arrays of input data. It is called from DATA.

• READTABLE -- read pre-calculated thermodynamic properties from tables.data and phase boundary data from phasebd.data. It is called from GTH.

• RELPERM – calculate relative permeability for each control volume. It is called from PROPTY.

• RESET – re-initialize variables and properties if a pressure-enthalpy iteration does not converge. It is called from GTH.

• TIMEPARA – read parameters for different time intervals. It is called from GTH.

• SOLVEQ – solve finite different equations for pressure or enthalpy. It is called from GTH.

• SOLVLINE – solve finite different equations for a specific line, row or column. It is called from LBYL.

• T2PHASE – calculate two-phase equilibrium temperature and its derivatives at a given pressure. It is called from DENPT, PROPTY and INTFPARA.

• TDMA – apply the Thomas algorithm to solve finite different equations. It is called from SOLVLINE.

• TSTEPSET – write information on time step setting. It is called from GTH.

• TTOH – convert temperature to enthalpy when an initial state does not involve any two-phase region and temperature rather than enthalpy values are read. It is called from DATA.

• UPSTREAM – apply upstream weighting on advection-related properties, especially relative permeability, at each interface. It is called from INTERFACE.

• VARIATION – track the variation of mass and energy sources/sinks and permeability. It is called from GTH.

• VELOCITY – calculate the velocity and mass flow rate of fluid flow according to Darcy’s lay, also calculate saturation wave speed if the corresponding interface is located in a two-phase region. It is called from INTERFACE.

• VISC2PHASE – calculate liquid and vapor viscosity values at a given pressure. It is called from PROPTY, INTFPARA and UPSTREAM.

• VISCOSITY – calculate fluid viscosity as a function of temperature and density. It is called from PROPTY and INTFPARA.

• WR – print what read by READ and RD. It is called from PDATA.

Programming

GTH was written in FORTRAN-77. It has been run on PC, UNIX workstations and the SGI Origin 2000 supercomputer at Georgia Tech.

Users should first download all files listed in the section of Brief Description of Subroutines to a subdirectory of their computer. You may want to make changes to file GEOTH.PAR according to the size of the numerical model before compilation to save computer space. These FORTRAN files are then compiled using a FORTRAN compiler to create an executable file named GTH or whatever name you like. If you have a UNIX workstation you may want to use the provided file MAKEFILE to do compilation work. Before the compilation a subdirectory "run" should be created for running the program. The executable file created by compilation should be put in the "run" subdirectory as well all the data files needed for a program run.

Input files are prepared for establishing look-up tables of thermodynamic properties, reading problem specifications, and adjusting program-control parameters.

• tables.data – It is read by READTABLE. These are two-dimensional data of temperature, density, viscosity and density derivatives with respect to pressure and enthalpy calculated as functions of given pressure and enthalpy values for single-phase, either liquid or vapor, pure water using Haar et al.’s [1984] equation of state. The whole single-phase region is divided into five sub-regions. For each of the sub-region, pressure and enthalpy entries are read first and then the properties.

• phasebd.data – It is read by READTABLE. These are one-dimensional density, enthalpy and viscosity data for saturated liquid and vapor calculated as functions of given pressure values. Pressure entries are read first and then the properties.

• gxin – It is read by DATA and TIMEPARA. These are the data specified for the problem to be solved and the parameters controlling the program run. They are discussed in greater detail in Appendix B.

• bc.data – It is read by BC. These are boundary conditions set for the problem to be solved. They are discussed in greater detail in Appendix C.

• bgsosi.data – It is read by BGSOSI. These are the mass and energy sources/sinks of the problem to be solved. It also includes the information on time-dependent variation of permeability and those sources/sinks. They are discussed in greater detail in Appendix D.

There are two output files, G3out0 and G3pickup, that are automatically produced by each GTH run. Users may also choose to have up to four other output files, namely G3out1, G3out2, G3out3, and G3out4.

• G3out0 – saves the information of a GTH run process. This file is produced for each GTH run.

• G3pickup – saves simulation results after each time step. It may be used to check the status of a GTH run or continue a run that was terminated for some reason. This file is produced for each GTH run.

• G3out1 – saves simulation results of thermodynamic properties and phase equilibrium calculations at times preset by parameters read from gxin.

• G3out2 – saves liquid flow rate results at times preset by parameters read from gxin.

• G3out3 – saves vapor flow rate results at times preset by parameters read from gxin.

• G3out4 – saves information for specified control volumes at times preset by parameters read from gxin.

Example Problems

Appendices

Appendix A. List of Symbols

Symbols:

Other Subscript:

Superscript:

References