- Email: [email protected]
Validation of Finite Element Simulation of the Hydrothermal Behavior of an Artificial Aquifer Against Field Performance
Validation of Finite Element Simulation of the Hydrothermal Behavior of an Artificial Aquifer Against Field Performance H. Daniels Institut fuer Wasserbau und Wasserwirtschaft, Aachen University of Technology, Federal Republic of Germany
ABSTRACT The effect of natural convection and density driven flow is of major importance to the flow field and the temperature field in the artificial aquifer for thermal energy storage that has been built at the Stuttgart University. This paper describes the Finite Element Approximation that has been used to calculate the transient temperature field in the aquifer and shows some validation results against measurement data from field experiments in Stuttgart.
THE TEST SITE At the Institute for Thermodynamics of the Stuttgart University an artificial aquifer has been constructed within the German underground thermal energy storage program. Fig. 1 shows a vertical cross section of the axisymmetric aquifer. Cold water is exchanged through holes near the bottom of the pressure chamber in the central tube while hot water is fed into or released from the aquifer by a drainage ring conduit at a radius of about 9,70 m at the top of the aquifer. A storage layer of gravel 8-16 mm is bounded vertically by two charging layers of gravel with 16-32 mm grain size. This vertical layering as well as the location of the water inlet facilites was emplaced as a result of design studies for different geometries of the aquifer that were done in Aachen. The Stuttgart aquifer contains several hundred temperature measuring probes. The thermal performance of the test aquifer can be observed exactly. Finite Element Analysis of the aquifer's hydrothermal behavior was fairly complicated because of the huge permeability of the gravel filling and temperature differences up to 45OC in the water-saturated aquifer.
319
320
I
Fig. 1: Vertical cross section of the artificial aquifer MATHEMATICS OF ENERGY TRANSPORT IN POROUS MEDIA FLOW The description of transient energy transport in porous media is well understood since more than one decade. Suppose the velocity field vi is given, the convection-diffusionequation can be written:
In eq. (1) local equilibrium between fluid and matrix temperature is assumed. As the velocity field is not known a priori, it has to be calculated before the convective transport can be evaluated. Major difficulties arose and efforts had to be made when calculating the strongly density dependent flow field in the artificial aquifer. Different mathematical formulations were used in order to find an appropriate approximation for the hydraulic behavior of the artificial aquifer. Hvbrid Formulation Assuming DARCYs law to be valid and substituting it into the continuity equation ends up with a hybrid formulation for the pressure distribution p
Multi Arrav Formulation If DARCYs law v
K
i
az
- - A ( s + p f g F ) "f j j
321 is not inserted into the continuity equation
more than one independent variable, i.e. p aqd vi, are to be calculated simultaneously. Ouasi NAVIER-STOKES Formulation Analogous to the laminar viscous flow problem, porous media flow can be written in a NAVIER-STOKES Formulation for isotropic permeability (GARTLING 1 9 8 7 ) . In this case eq. ( 3 ) can be replaced by
Other than eq. (2) and ( 3 ) , eq. (5) takes viscous forces into account and does not necessarily require a zero-order potential field as standard DARCYs law does, Discussion on Buovancv Flow A zero-order potential field F in the way
-
does not exist if the fluid density pf is variable in space, because we obtain rot
=
grad p x grad
1 z 0
pf
(7)
in this case (STOSSINGER 1 9 7 9 ) . Only if the spacial variability of p is small, the flow field can approximately be handled as in eq. ( 2 ) and ( 3 ) . In terms of Finite Element Analysis the element spacing necessary is strongly dependent on grad l/p and at the same time on the value of grad l/p compared to the value of grad p or the matrix permeability K, respectively. NUMERICAL SOLUTION The Finite Element Formulations based on the BUBNOV-GALERKIN method for eqs. (1) to ( 4 ) , are assumed to be familiar. Here the Finite Element form of the quasi NAVIER-STOKES flow eq. ( 5 ) will be considered, Using the GAUSS theorem to reduce the second order diffusion terms in eq. (5) the GALERKIN formulation finally leads to
322
FINITE ELEMENT GRID Because of the domination of buoyancy flow in the artificial aquifer, a Finite Element grid that consists of equilateral triangles proved to be the most appropriate discretization. Fig. 2 shows two ideas for grids that performed about the same accuracy during model validation. Grid a), developed according to KINZELBACH and FRIND (1986) did not only have twice as much nodes than grid b), but also required a much smaller time step , to obtain a numerically stable solution. I
I6
IL
ID
11
7
II
6
10
5
29
'
28
I
27
I5 I&
ldzl
1
I1
2
I2
=
2
26
I
2s
9
17
I
Fig. 2 : Sketch of grids
6
11
Fig. 3 : Centres of mass
Even more obvious is another advantage. Fig. 3 shows the centre of mass for adjacent elements in grid a) and grid b). There are always heavy and light horizontally neighboring elements in grid a), when a vertical layering of water with different temperature is assumed. Once two nodes are warm and one is cold and the element is light, another time one node is warm and two nodes are cold and the element is heavy. A vertical temperature layering and a stable flow field are very hard to obtain with grid a). There are no such problems with grid b). RESULTS OF MODEL VALIDATION Fig. 4 shows the parameters used during validation of the numerical model for the aquifer in Fig. 1. Aquifer permeability is very high. The relation of longitudinal to transversal
323
K [m21 charging layers storage layer
4 days
n
neff
[-I [-I
a~
[ml
I.[ k1 p] [k] aI1
b'
b'
r
m
b'
2.107 0,37 0,34 0,25 0,125
2,2
4.108 0 , 3 7
2,2 1.910 1.994
0,33 0,25 0,125
6 days
Fig. 5a: Two-dimensional temperature field
1.910 1.994
324
Fig, 5b: Temperature over time measured and calculated CONCLUSIONS Permeability and temperature differences are very high in the Stuttgart aquifer test site. The temperature field can be calculated accurately using equilateral triangles according to Fig. 2. A quasi NAVIER-STOKES formulation for the flow field equation may be advantageous compared to standard DARCY as far as natural convection in coarse gravel is concerned. A comparison of accuracy and computational effort between different flow field formulations will be given elsewhere. REFERENCES
1. Fisch N. (1986), Documentation of the Stuttgart University Project, International Agency, Solar Heating and Cooling Program, prepared for Task VII, November 1986 Meeting. 2 . Gartling D.K. (1987), NACHOS I1 - A Finite Element Computer Program for Incompressible Flow Problems PART I-Theoretical Background, SAND 86-1816.UC-32,Sandia Nat. Lab. Albuquerque 3 . Kinzelbach W.K.H., Frind E.O. (1986), Accuracy Criteria for Advection - Dispersion Models, VI International Conference on Finite Elements in Water Resources, Lisboa. 4 . Pinder G.F. and Gray W.G. (1977), Finite Element Simulation
in Surface and Subsurface Hydrology, Acad. Press, New York
5 . St6ssinger W. (1979), Beschreibung der Hydrodynamischen
Dispersion mit der Methode der Finiten Elemente am Beispiel der instationaren Interface zwischen SQB- und Salzwasser in Grundwasserleitern, Mitt. I W , RWTH Aachen, Heft 28.
6. Zienkiewicz D.C. and Morgan K. (1983), Finite Elements and Approximations, A Wiley Interscience Publication, New York.
ABSTRACT The effect of natural convection and density driven flow is of major importance to the flow field and the temperature field in the artificial aquifer for thermal energy storage that has been built at the Stuttgart University. This paper describes the Finite Element Approximation that has been used to calculate the transient temperature field in the aquifer and shows some validation results against measurement data from field experiments in Stuttgart.
THE TEST SITE At the Institute for Thermodynamics of the Stuttgart University an artificial aquifer has been constructed within the German underground thermal energy storage program. Fig. 1 shows a vertical cross section of the axisymmetric aquifer. Cold water is exchanged through holes near the bottom of the pressure chamber in the central tube while hot water is fed into or released from the aquifer by a drainage ring conduit at a radius of about 9,70 m at the top of the aquifer. A storage layer of gravel 8-16 mm is bounded vertically by two charging layers of gravel with 16-32 mm grain size. This vertical layering as well as the location of the water inlet facilites was emplaced as a result of design studies for different geometries of the aquifer that were done in Aachen. The Stuttgart aquifer contains several hundred temperature measuring probes. The thermal performance of the test aquifer can be observed exactly. Finite Element Analysis of the aquifer's hydrothermal behavior was fairly complicated because of the huge permeability of the gravel filling and temperature differences up to 45OC in the water-saturated aquifer.
319
320
I
Fig. 1: Vertical cross section of the artificial aquifer MATHEMATICS OF ENERGY TRANSPORT IN POROUS MEDIA FLOW The description of transient energy transport in porous media is well understood since more than one decade. Suppose the velocity field vi is given, the convection-diffusionequation can be written:
In eq. (1) local equilibrium between fluid and matrix temperature is assumed. As the velocity field is not known a priori, it has to be calculated before the convective transport can be evaluated. Major difficulties arose and efforts had to be made when calculating the strongly density dependent flow field in the artificial aquifer. Different mathematical formulations were used in order to find an appropriate approximation for the hydraulic behavior of the artificial aquifer. Hvbrid Formulation Assuming DARCYs law to be valid and substituting it into the continuity equation ends up with a hybrid formulation for the pressure distribution p
Multi Arrav Formulation If DARCYs law v
K
i
az
- - A ( s + p f g F ) "f j j
321 is not inserted into the continuity equation
more than one independent variable, i.e. p aqd vi, are to be calculated simultaneously. Ouasi NAVIER-STOKES Formulation Analogous to the laminar viscous flow problem, porous media flow can be written in a NAVIER-STOKES Formulation for isotropic permeability (GARTLING 1 9 8 7 ) . In this case eq. ( 3 ) can be replaced by
Other than eq. (2) and ( 3 ) , eq. (5) takes viscous forces into account and does not necessarily require a zero-order potential field as standard DARCYs law does, Discussion on Buovancv Flow A zero-order potential field F in the way
-
does not exist if the fluid density pf is variable in space, because we obtain rot
=
grad p x grad
1 z 0
pf
(7)
in this case (STOSSINGER 1 9 7 9 ) . Only if the spacial variability of p is small, the flow field can approximately be handled as in eq. ( 2 ) and ( 3 ) . In terms of Finite Element Analysis the element spacing necessary is strongly dependent on grad l/p and at the same time on the value of grad l/p compared to the value of grad p or the matrix permeability K, respectively. NUMERICAL SOLUTION The Finite Element Formulations based on the BUBNOV-GALERKIN method for eqs. (1) to ( 4 ) , are assumed to be familiar. Here the Finite Element form of the quasi NAVIER-STOKES flow eq. ( 5 ) will be considered, Using the GAUSS theorem to reduce the second order diffusion terms in eq. (5) the GALERKIN formulation finally leads to
322
FINITE ELEMENT GRID Because of the domination of buoyancy flow in the artificial aquifer, a Finite Element grid that consists of equilateral triangles proved to be the most appropriate discretization. Fig. 2 shows two ideas for grids that performed about the same accuracy during model validation. Grid a), developed according to KINZELBACH and FRIND (1986) did not only have twice as much nodes than grid b), but also required a much smaller time step , to obtain a numerically stable solution. I
I6
IL
ID
11
7
II
6
10
5
29
'
28
I
27
I5 I&
ldzl
1
I1
2
I2
=
2
26
I
2s
9
17
I
Fig. 2 : Sketch of grids
6
11
Fig. 3 : Centres of mass
Even more obvious is another advantage. Fig. 3 shows the centre of mass for adjacent elements in grid a) and grid b). There are always heavy and light horizontally neighboring elements in grid a), when a vertical layering of water with different temperature is assumed. Once two nodes are warm and one is cold and the element is light, another time one node is warm and two nodes are cold and the element is heavy. A vertical temperature layering and a stable flow field are very hard to obtain with grid a). There are no such problems with grid b). RESULTS OF MODEL VALIDATION Fig. 4 shows the parameters used during validation of the numerical model for the aquifer in Fig. 1. Aquifer permeability is very high. The relation of longitudinal to transversal
323
K [m21 charging layers storage layer
4 days
n
neff
[-I [-I
a~
[ml
I.[ k1 p] [k] aI1
b'
b'
r
m
b'
2.107 0,37 0,34 0,25 0,125
2,2
4.108 0 , 3 7
2,2 1.910 1.994
0,33 0,25 0,125
6 days
Fig. 5a: Two-dimensional temperature field
1.910 1.994
324
Fig, 5b: Temperature over time measured and calculated CONCLUSIONS Permeability and temperature differences are very high in the Stuttgart aquifer test site. The temperature field can be calculated accurately using equilateral triangles according to Fig. 2. A quasi NAVIER-STOKES formulation for the flow field equation may be advantageous compared to standard DARCY as far as natural convection in coarse gravel is concerned. A comparison of accuracy and computational effort between different flow field formulations will be given elsewhere. REFERENCES
1. Fisch N. (1986), Documentation of the Stuttgart University Project, International Agency, Solar Heating and Cooling Program, prepared for Task VII, November 1986 Meeting. 2 . Gartling D.K. (1987), NACHOS I1 - A Finite Element Computer Program for Incompressible Flow Problems PART I-Theoretical Background, SAND 86-1816.UC-32,Sandia Nat. Lab. Albuquerque 3 . Kinzelbach W.K.H., Frind E.O. (1986), Accuracy Criteria for Advection - Dispersion Models, VI International Conference on Finite Elements in Water Resources, Lisboa. 4 . Pinder G.F. and Gray W.G. (1977), Finite Element Simulation
in Surface and Subsurface Hydrology, Acad. Press, New York
5 . St6ssinger W. (1979), Beschreibung der Hydrodynamischen
Dispersion mit der Methode der Finiten Elemente am Beispiel der instationaren Interface zwischen SQB- und Salzwasser in Grundwasserleitern, Mitt. I W , RWTH Aachen, Heft 28.
6. Zienkiewicz D.C. and Morgan K. (1983), Finite Elements and Approximations, A Wiley Interscience Publication, New York.