SIMULATION ON THERMAL STRATIFICATION OF THE HUGE-CUBAGE AND DEEP RESERVOIRS*

 

 

Deng Yun1, Zhao Wenqian2, Li Jia3 and Luo Lin4

1Ph.D Candidate, 2,3,4 Professors

Sichuan University, the state key laboratory of high speed flow

Tel: 028-5407780; Fax: 5401558; E-mail: dengyun@fm365.com

 

 

Abstract: In this study the laterally averaged two-dimensional kε turbulent model is firstly tested in Ertan reservoir and the results are compared favorably with the observed measurements in site. Then it is applied to predict the thermal stratification for another reservoir. The forming, growth and erosion of the thermocline are simulated and the development of the thermal structure is analyzed all-around.

 

Keywords: thermal stratification, mathematical simulation, reservoir

1    INTRODUCTION

The thermal stratification is one of important features, which has great impact on the ecological environment in the reservoir area and the downstream. How to simulate and predict the thermal structure for a huge reservoir has become a key problem for the reservoir design and management. The buoyancy due to the temperature difference greatly influences the flow field, which brings great difficulties to the computation. Jiang(2000) calibrated the turbulent viscosity in different direction for different variables with the observed data. Chu(1997) expressed the turbulent viscosity as the function of velocity gradient, thermal stratification parameter and turbulent length scale. The kε turbulent model with buoyancy modification has great advantages in stabilization and few parameters which is adopted by many researchers, such as Chen(1997) and Gerard(1988). For the simulation of a large reservoir in a long period, the cost of three-dimensional computation is too expensive. In the narrow reservoir that the variations of the variables in the transverse direction is not obvious, the laterally averaged two-dimensional model is a good choice which can well simulate the motion of the buoyant flow in the longitudinal and vertical directions and the development of the vertical thermal stratification. In this study, the laterally averaged two-dimensional kε turbulent model is verified with observed field data and applied to predict the thermal stratification in a huge reservoir.

2    MATHEMATICAL MODEL

2.1  The governing equations

2.1.1  The state equation

For the water in nature, the density can be expressed

    (1)

The temperature difference is the buoyancy source. When the temperature difference is small, the buoyancy (gravity term) can be written as the linear function of the temperature difference

,                    (2)

where  and  is the reference temperature and density;  is the coefficient of thermal expansion.

The Boussinesq approximation is used in this study, which means that the buoyant effect is considered only in the gravity term.

2.1.2  The continuity equation

                            (3)

2.1.3  The momentum equations

              (4)

          (5)

2.1.4  The temperature equation

    (6)

In the above equations u and w are the mean longitudinal and vertical velocities respectively; T is the mean temperature; P is the pressure; B is the river width; , where  is the molecular viscosity and  is the kinematic eddy viscosity;  is the turbulent Prandtl number.  is specific heat of water;  is the solar radiation penetrating into the water.

2.2    The turbulent model

The standard kε turbulent model is used with modification by buoyancy forces. The eddy viscosity  is computed from the equation

                               (7)

where k and ε are the turbulent kinetic energy and the kinetic dissipation rate respectively.

The k and ε equations for the buoyant flow are:

  (8)

   (9)

in which  is the production of kinetic energy from the mean flow and is the production or destruction of turbulent kinetic energy by buoyant forces, which are given by the expressions

,          (10)

The buoyancy tends to suppress turbulence and turbulent diffusion in the stably thermal stratification environment and boost up in reverse.

, , the Prandtl number of the kinetic and kinetic dissipation rate, are given for the standard value of 1.0 and 1.3. The other constants are :

2.3    Boundary conditions

The symmetry conditions for velocities are given on the free surface. The reservoir bottom and the dam face are taken as no-slip boundary. At the inflow and outflow boundary u is set as a constant and w as zero.

The bottom of reservoir and the surface of dam are taken as adiabatic. The inflow temperature is given a constant value and a zero gradient condition is set at the outflow boundary. The conditions on the free surface are much more complicated, since the heat exchange is the main heat source of reservoir. The heat fluxes through the water surface mainly include five parts

                     (11)

The net solar radiation                                        (12)

where   is incoming solar radiation; (0.1) is the water reflection coefficient, and (0.65) is the fraction of solar radiation absorbed at surface. The remaining fraction of the solar radiation  is absorbed exponentially with depth

                       (13)

where H is the water depth and (0.5) is extinction coefficient.

The atmospheric longwave radiation                       (14)

where Ta is air temperature (Celsius degrees), (=5.67×10-8W/m2·K4) is Stefan-Boltzmann constant, and  is atmospheric emissivity 

K (0.17) is the cloud height coefficient and C is the cloud cover ratio.

The back radiation from water to the air follows the same formulation with air temperature replaced by the water surface temperature Ts

                       (15)

and the emissivity  is fixed at 0.975

The evaporative loss (latent heat)                              (16)

where es and ea are the saturation vapor pressure above the water surface and the vapor pressure; W is the air velocity 10 meters above the water surface. The function of W is expressed

The conductive loss (sensible heat)                       (17)

Zero gradient conditions for k andε are given on the free surface, the reservoir bottom, the outflow boundary and the dam surface. At the inflow boundary, k andε are computed by

,                       (18)

where u0 and H0 are the water velocity and the water depth at the inflow boundary.

2.4    Numerical computation

SIMPLE algorithm with staggered grid is used in this study (Patankar, 1980). Integration of governing equations over the control volume gives the discretized equations that are solved by TDMA (Tri-Diagonal Matrix Algorithm).

3    APPLICATION AND RESULT DISCUSSION

In order to verify the kε turbulent model for prediction of thermal stratification in the huge and deep reservoir, Ertan Reservoir is chosen as a test example. Then the model is applied to another huge reservoir, which is ongoing, for prediction its thermal structure.

3.1   Model verification

Ertan reservoir locates in Southwest of China. The length is about 136 km, width about 400 m and maximum water depth about 190m with a huge storage capacity of 4.478×109 m3. It is a long, narrow, deep and huge reservoir, and is suitable to apply laterally averaged two-dimensional stratified model. The bottom of dam is at 1010m above the sea level. The maximum water level in the reservoir is 1200m in dry season. The intake of power station is located at 1128m with 10m high. In flood season, the superfluous water is discarded from the spillways and discharge orifices.

The temperature distribution in the reservoir is measured in January 14th of 2000, which is taken as initial conditions of temperature. Since the simultaneous velocity is not known, the velocity of non-buoyant flow is computed as the initial conditions. For lack of the simultaneous meteorological data, the monthly averaged meteorological data in 1999 at Shengli station is used and the long-time average water temperature at the end of reservoir upstream is taken as inflow boundary condition. The flow rate and water level are all monthly average shown in Fig. 1.

The reservoir area is discretized into 368×67 rectangle elements, the size varying from 10 to 400 meters in longitudinal direction and 2 to 3 meters in vertical direction. The fixed mesh is produced for each month with different water level.

Fig.1    Variation of discharge and water level with time

Numerical simulations from January to July of 2000 are performed and compared with the field observed data. Figure 2 shows the initial temperature distribution in January 14th. Figure 3(a)-(f) compares the computed and measured temperature distribution in the vertical cross section 700 meters away from the dam site upstream. It is surprised that they coincide fairly well especially in March and April under the conditions that the meteorological data and boundary conditions are not so accuracy. The results from May to July are not so good, due to the frequent reservoir operation (as shown in Fig. 1).

Fig.2    Temperature Distribution in January 14th

Figure 4 compares the computed outflow temperature and the temperature measured at the downstream hydrometric station. The result is similar with Fig. 3 that the lines of the computed and measured coincide very well from January to April and the computed temperature curve is much more smoothly than the measured. Before May, though the water level decreases fast, the flow rate varies smoothly and the water is all discharged from the intake of power station, so the outflow temperature rises steadily and is simulated well. The temperature downstream rises abruptly in late of May which conforms with the flood discharging from the tunnel spillway and the sequent frequent reservoir operation leads to fluctuating of the outflow temperature. The reservoir operation has great impact on outflow temperature, which can be used to control the outflow temperature.

Fig.3    Comparison of temperature vertical profile

¿ measured      ——computed

Fig.4    Comparison of outflow temperature

In this model, only the solar radiation fraction and the extinction coefficient are needed to be determined. Through several trails, 0.65 and 0.5 are chosen for them. In fact they just appreciably influences the temperature of the water in 10 meters below the surface. Few parameters make the model very suitable for prediction.

3.2  Application

The temperature model is applied to a huge reservoir, which is ongoing to construct. Because there are important fisheries downstream, the temperature of discharged water from reservoir is greatly concerned.

The reservoir is 200 km long and about 600 m wide. It is a typical river-like reservoir. The maximum water depth is 220 m with a total storage capacity of 1.16×1010 m3. The intake of its power station is located at 67 meters bellow the normal pool level with 10 m high, from which water are almost discharged to the downstream in normal year.

The reservoir area is divided into 527×71 rectangle elements, the size of which is the same as Ertan. It is observed from Ertan reservoir that the temperature tends to be uniform in February. So the inflow temperature in February is taken as the initial temperature in the whole reservoir area, and the velocity of non-buoyant flow as initial velocity. The corresponding meteorological data, inflow rate and water level are adopted. In order to eliminate the error of initial conditions, after the one year computation, the computation are repeated from the newly computed temperature and velocity in January until the temperature difference between the adjacent iterations is small enough, and then the temperature computed in the newest iteration is taken as the final result.

Figure 5 shows the development of the temperature structure in four months, representing four seasons, in the normal year. The temperature tends to be uniform in February with the largest difference about 2 in vertical profile. In May with the increasing of inflow temperature and heat flux from surface, the upper water got warmer fast with the great temperature gradient forming at upper. But the water temperature still keeps low near the bottom. The apparent thermal stratification forms. In August the inflow and air temperature reach ultimate; there exists double temperature gradient in the vertical profile. Due to the heat exchange between water and air, a relatively small temperature gradient forms in a thin layer at surface. Beneath it, the increasing flow rate enforces the water convection on vertical section, so a warm and thermally uniform layer forms near the outlet. The cold water in deep bottom is slowly warmed up and a thermocline layer exists between the warm layer and the cold layer. In November, the air and the inflow temperature decreases soon, the water loses heat to air and gets cold soon. With cold-water gets heavier and goes down, the thermocline near surface disappears. At the same time, the inflow gets cold and goes down along the bottom, which makes the water near the bottom warmed up. The temperature difference declines soon. With the weather gets colder, the temperature declines in whole vertical profile until it tends to be uniform in the second February.

Figure 6 shows the courses of air temperature, inflow temperature and the predicted water temperature at surface and bottom at the dam site, and compared the outflow temperature before and after reservoir built. The temperature at surface varies soon, due to the direct heat exchange with the air, and the temperature variation is small at bottom less than 4. Thermal stratification is apparent in the huge reservoir. On the other hand, the temperature deferment of the huge reservoir is obvious. February and October are the cross junctions. From March to September after reservoir built, the outflow temperature declines 1.5 averagely, with the maximum of 3.5 in April. From November to the next January, the temperature rises 2.8 than before with the maximum of 3.5 in December. The thermal storage capacity also makes the maximum temperature difference in the whole year decline from 12.4 to 9.9. Though the mean temperature in one year is almost not changed, the temperature course changes much, which possibly influences the ecological environment at downstream.

Fig.5    Temperature structure in different time

xDistance from upstream(m)          ywater level(m)

Fig.6     Temperature variations with time()

4    CONCLUSION

The mathematic model described in this paper could give realistic simulations of development, growth and erosion of the thermal stratification in the reservoir. Few parameters need to be calibrated. They are insensitive to the thermal structure for the deep reservoirs, which makes the model suitable for prediction.

The computed results in the two huge and deep reservoirs show that the thermal stratification exists almost in the whole year and displays the different temperature profiles in different seasons. The surface temperature is influenced directly by the climate conditions and varies fast. Reversely, the bottom temperature keeps cold in whole year and reaches the warmest in the autumn. The temperature deferment is obvious for the huge reservoir especially in the April and December. As to the outflow temperature, the temperature difference in one-year declines and its variation changes greatly, which may influence the ecological environment at downstream.

References

[1]  Jiang Chun-bo, Zhang Qing-hai, Gao Zhong-xin. A 2-D unsteady flow model for predicting temperature and pollutant distribution in vertical cross section of a river. Journal of Hydraulic Engineering, Chinese. 2000, (9): 20-24.

[2]  C.R. CHU and C.K. SOONG. Numerical simulation of wind-induced entrainment in a stably stratified water basin. J Hydraulic Research, Vol. 35, 1997, No.1.

[3]  Gerard J. Farrell, Heinz G. Stefan. Mathematical modeling of plunging reservoir, J. Hydraulic Research. Vol. 26, No. 5, 1988.

[4]  Chen Xiaohong, Liu Meinan, Lin Yanshan. Study on two-dimensional water quality distribution in reservoirs. Journal of Hydraulic Engineering, Chinese. 1997, (4): 9-16.

[5]  PATANKAR, S.V.(1980). Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.



* Supported by the Key Program of the National Natural Science Foundation of China (Granted No. 59639240)