Abstract
This work brings as its main contribution the study of the phenomenon of density currents in nonstratified reservoirs, with flows in steady regimes. Flows are analyzed for a wide range of Reynolds and Richardson numbers in the entrance of the reservoir. Based on this hypothesis, a hybrid numeric model is presented considering the Reynolds Transport's Equation  focusing on the conservation for volume, mass and momentum  with the intention of obtaining threedimensional components of velocities, reduced acceleration of gravity and geometric characteristics of currents along the reservoir. It is possible to notice in the numeric simulations, mainly, the need of complementation of the model that refers to the inclusion of the entrainment coefficient and the analysis in the unsteady regime.
density current; laser velocimetry; stratified flowsyword
TECHNICAL PAPERS
Density currents at steady regime
Alexandre Augusto Barbosa^{I}; Frederico Fábio Mauad^{II}; Carlos Eduardo Sato^{III}; Ana Paula Moni Silva^{IV}
^{I}Universidade Federal de Itajubá  UNIFEI Instituto de Recursos Naturais 37500903 Itajubá, MG, Brazil. barbosa@unifei.edu.br
^{II}Universidade de São Paulo  USP Escola de Engenharia de São Carlos  EESC Centro de Recursos Hídricos e Ecologia Aplicada 13560970 São Carlos, SP, Brazil. mauadffm@sc.usp.br
^{III} carlos@setorverde.org.br
^{IV}Universidade Federal de Itajubá  UNIFEI Instituto de Recursos Naturais 37500903 Itajubá, MG, Brazil. anamoni@unifei.edu.br
ABSTRACT
This work brings as its main contribution the study of the phenomenon of density currents in nonstratified reservoirs, with flows in steady regimes. Flows are analyzed for a wide range of Reynolds and Richardson numbers in the entrance of the reservoir. Based on this hypothesis, a hybrid numeric model is presented considering the Reynolds Transport's Equation  focusing on the conservation for volume, mass and momentum  with the intention of obtaining threedimensional components of velocities, reduced acceleration of gravity and geometric characteristics of currents along the reservoir. It is possible to notice in the numeric simulations, mainly, the need of complementation of the model that refers to the inclusion of the entrainment coefficient and the analysis in the unsteady regime.
Keywords: density current, laser velocimetry, stratified flowsyword
Introduction
According to Alavian et al. (1992) within the area of Environmental Fluid Mechanics density currents represent an important subarea of the stratified turbulent flows. Those are nothing else than the insertion of the flow of a fluid into a reservoir that contains another fluid of lower or higher density than the first one, i.e. density currents or gravity currents are the result of the interaction of two or more fluids of different density  a widely common natural phenomenon.
These density differences may be the result of fluids having different temperatures, concentrations of sediments in suspension or concentration of dissolved solids, of different salinity or distinct fluids (Birman and Meiburg, 2006).
The subject can be studied in several areas of science and engineering, cited by Tsihrintzis and Alavian (1997): geophysics, hydraulics, limnology, heat and mass transfer and weather forecast. Its application within the environmental area has been drawing significant attention in recent times.
We may see density currents resulting from dense and not very dense inflows, whether they occur in stratified lakes or reservoirs or not.
Density currents are governed by the following equations: momentum, continuity and volume conservations, mentioned by Hauenstein and Dracos (1984). They also include forcing functions, boundary conditions and mixing locating regions. A combination of laboratory and field experiences and analytical and/or numerical approximations are necessary for a better understanding of these phenomena.
Field experiments can be carried out for situations when the geometric features of the reservoirs and inflows have little complexity. Such experiments must attain measures of the general condition of the reservoir and of the behavior of the density currents and, simultaneously, detailed measures for the definition of mixture turbulent processes (according to Barbosa, 1999; Alavian et al., 1992). It is important to evolve to a greater understanding of how the phenomenon can influence the transport and the mixture of materials, on the heat transfer and on the dissolved and suspended substances.
It is extremely urgent to carry out a reasonable number of laboratory experiments. This offers an important advantage regarding controlled visual inspections, programmable conditions and rapid flow variations for the simulation of density currents. The input data and those attained by these experiments are essential for the modeling that may be used for a better arrangement of the reservoirs and environmental forecasts that are, at least, minimally reliable.
The development of currents in channels with variable inclination and in stratified water environments (with different densities) demands more studies, given that they are localized events. In these cases, the modern optical devices for twodimension instantaneous observations offer new opportunities for data collections, studied by Barbosa (1999). For example, it is possible to mention the LDV (Laser Doppler Velocimetry), the PIV (Particle Image Velocimetry) and the LIF (Laser Induced Fluorescence), among others.
In addition, new advances in the development of analytical and/or numerical modeling of the currents inside the reservoirs, conciliating the large scale behavior and the appearance of localized turbulent events, are still necessary (Akiyama et al., 1994).
This study may serve as a basis to connect density current models and global dynamic behavior models of lakes and reservoirs.
Nomenclature
A = crosssection area of density current, cm^{2} b = width of density current, cm Cƒ = drag coefficient, dimensionless d = infinitesimal quantities E = entrainment coefficient, dimensionless F = force, N Fr = Froude number, dimensionless ƒsi = proportionality factor for the implicit simulation g' = reduced gravity acceleration, m/s ^{2} h = thickness (depth) of density current, cm L = length of reservoir, cm Q = volumetric flux, cm^{3}/s Re = Reynolds number, dimensionless Rh = hydraulic radius, cm Ri = Richardson number, dimensionless S = declivity, dimensionless SC = closed surface t = time, s u = xcomponent velocity, cm/s v = ycomponent velocity, cm/s V = Volume, cm ^{3} Ve = Vendernikov number, dimensionless V = velocity field, cm/s VC = control volume w = zcomponent velocity, cm/s m = mass, g s = direction, surfaceGreek Symbols
ρ = density, kg/m ^{3} θ = slope angle, degreesSubscripts
ms relative to mass ss relative to surface r reservoir sub submersion point o entrance of reservoir 1 relative to the internal diameter of the internal tube 2 relative to the external diameter of the internal tube 3 relative to the internal diameter of the external tubeExplaining the Phenomenon
The main studies on the insertion of affluents in reservoirs that have been developed so far deal with density currents that spread in the inclined bottom of the reservoir, as shown in Fig. 1. It is possible to observe that for a determined depth the flow of the density current gets separated in the bottom. This takes place due to the fact that the density of the fluid of the reservoir is the same or higher than the density of the current. It is important to highlight that many times not all of the five flow zones manifest themselves in all of the cases. It is possible, for example, the occurrence of only two of them (zones 2 and 3).
The approximation region (zone 1) consists of a density flow ρ_{o} and can be analyzed by using the hydraulic approach of channels with free surface, using experimental data and empirical relations, according to Chow (1959). Zones 4 and 5 will not be portrayed by this study.
After establishing the density current, there is a flow steady regime (expected when there is amplification and overlapping of vortical waves).
The most complex cases regarding density currents can be separated into two domains, which are based on the Richardson number at the entrance. They are:
momentum domain;
density difference domain.
For extremely critical regimes at the entrance, these two domains will exist and for subcritical regimes, just the latter will be present.
Figure 2 shows both domains for a three dimensional density current.
The boundaries of the submersion point, h_{sub}, are given by Eq. (1), which refers to data of laboratory experiments:
This equation is valid to h_{sub} > h_{o}.
Both domains will be treated by onedimensional integral forms of the equations of volume, continuity and momentum applied to a prismatic control volume of infinitesimal length and finite width and thickness.
Phenomenon Governing Equations
Based on Reynolds Transport's Equation, the referred equations in their general forms are:
Volume conservation:
Mass conservation:
Momentum equation:
where:
VC is the control volume;
SC is closed surface A, which recovers the control volume;
is the field of flow velocity;
us is the velocity component of the referred direction;
Δρ is the density difference between the density current and the receptor reservoirs;
Fms and Fss are the mass and surface forces acting on the control volume towards the referred direction
Simplifying Hypotheses
Permanent and onedimensional flow;
Hydrostatic distribution of the pressures;
Isotropic and homogeneous turbulence;
Similar velocity and density profiles;
Incompressible flow;
Velocity and density uniform profiles at the entrance and exit of the control volume;
Drag at the liquidliquid interface is considerably small;
There is only moment flow through the perpendicular sections towards x.
Control Volume
The configuration of the control volume, shown in Fig. 3, is valid for both domains.
The additional variables in Fig. 3 are:
x, y and z refers to coordinate system;
u, v and w are velocities components of density current in directions x, y and z, respectively;
the subscript "r" refers to reservoir conditions;
infinitesimal quantities are denoted for differential "d".
In the first domain, the upper face of the control volume is under the action of the atmospheric pressure, and in the second, of the water column of the reservoir (the flow is submerged).
According to Barbosa (1999), solving the Eqs. (2), (3) and (4) for both domains, the differentials equations are:
Equations for the Control Volume  First Domain
Equations for the Control Volume  Second Domain
In order to simplify the numerical solutions of differentials equations, the explicit derivatives are established.
Explicit Differential Equations  First Domain
Explicit Differential Equations  Second Domain
In the equations above:
h is the depth;
b is the width;
g' is the reduced gravity;
E is the entrainment coefficient;
θ is the reservoir slope angle;
C_{} is the drag coefficient of current density.
Numerical Solution of the Explicit Differential Equations
According to Akyiama et al. (1994) and Barbosa (1999), the equations of the topic above are numerically solved through the regressive and central finite differences of first and second orders with implicit approximation.
The program used for the calculation of the previous equations was developed using Visual Basic WindowsExcel^{®}.
The program is selfexplanatory regarding the entrance variables:
fluid viscosity (assumed as constant and invariable);
geometric features of the entrance section, ho and bo;
velocity and the reduced acceleration of the inflow at the entrance, uo and g'o;
slope angle, θ;
length of the reservoir, L;
calculation grid (towards x);
length between the calculations points.
The calibration parameters are:
entrainment coefficient, E;
the calculation grid;
and the proportionality factor for the implicit simulation of the differential equations, si.
The output variables are presented in a table with values of the desired calculation distances for each calculation section:
longitudinal reduced acceleration g';
three components of the velocity (u, v and w) along x;
thickness h;
width b;
hydraulic radius, Rh;
and the Richardson (Ri), Reynolds (Re) and Vendernikov (Ve) numbers of the flow.
In addition, there is also the coefficient and the forces that can be calculated, but they were not presented.
Results
The simulations that were carried out were compared to experimental data developed by the authors using an experimental apparatus of the Engineering School of São Carlos  University of São Paulo.
All the simulations were carried out with longitudinal spacing (x direction) of the control volume of 0.01 cm. Some simulated quantities such as velocity components v and w, the hydraulic radius of the current and the Richardson, Reynolds and Vendernikov numbers, as shown in charts (c) and (d) of Figs. 4, 5, 6 and 7, do not present experimental proof, either because they are indirectly attained from basic variables or because this work does not deal with their experimental measures (in this case they fit as velocity components v and w);
Other variables could have also been plotted such as forces and drag coefficient. They were not, but the program allows their calculation.
Each simulation in the laboratory is identified by a number that has five digits representing:
the 1st digit shows the type of declivity: 1 for angles of 12.5 degrees and 2 for those ranging from 4 to 12.5 degrees;
the 2nd digit shows the type of roughness: 1 for smooth bottoms and 2 for rough ones;
the 3rd digit represents the density of the inflow mixture: 1 for d = 1.005, 2 for 1.015, 3 for 1.025 and 4 for 1.035;
the 4th digit represents the width of the channel or the hydraulic depth at the entrance: 1 for b = 4.9 cm, 2 for b = 2 cm and 3 for ho ≠ 0.7 cm;
the 5th digit shows the amount of inflow: from the smallest ones (1) to the quantity defined for each series (the maximum value of the last digit corresponds to the highest attained flow). In order to vary the flow of the inflow, an option was made towards the valve control. The flow can be varied from 4 to 160 cm^{3}/s.
Figures 4 to 7 display the simulations carried out and their comparisons with some experimental data.
Comments and Conclusions
Besides being second order for finite differences (first order for the first calculation station), the numerical treatment method still brings a better accuracy to most of the calculations by means of implicit approximations of the second member of the ordinary differential equations;
iterative calculations were carried out for all of the equations, aiming at refining the value of the referred variable for each calculation station;
the chosen basevariable to attain the refinement of the convergence is the longitudinal velocity u, for it is calculated through the most complex differential equation. Generally, the number of necessary iterations for such convergence is no higher than four;
the RungeKutta method was tested for some cases, but in relation to the experimental data it did not present good results, neither regarding convergence nor accuracy;
depending on the values adopted for the entrainment coefficient (E) and on the characteristics of the implicit calculation (si) of the forcing function, there are instabilities and not convergences for certain simulations. However, it is possible to work within a range of values for si and E where the convergence is attained without problems, and whose accuracy is in charge of the refinement of such parameters;
for the cases that are supercritical at the entrance, it is possible to work with si values between 0.45 and 0.50 and E values between 0.10 and 0.16. In such cases, the accuracy of the result does not vary significantly. It is also possible to work with si values between 0 and 0.05 and E between 0.065 and 0.08 in subcritical cases.
for cases of approximated forecast, it is possible to work with si = 0 and E = 0.075 for the subcritical cases and with si = 0.47 and E = 0.13 for the supercritical ones;
experiments 21433 and 21434 do not present submerging points. On the other hand, experiments 21435 and 21436 present them, and the program identifies them with satisfactory precision. The validity of the program does not lie on the fact regarding the identification of the submerging point, given that it is carried out at the beginning of the program through Eq. (1), but by the fact that the behavior of the features of the density currents follows what was experimentally observed, that is:
1. a growth until the submerging point and then a drop in the thickness of the current;
2. a constant drop and a reversion of the direction of component w of the velocity so that, right afterwards, it will come close to zero at large distances from the entrance of the reservoir;
for the cases where there is a hydraulic jump (the submersion itself), it is possible to notice from Figs. 6 and 7 that the tendency of velocity components w is to reduce at the area presenting a large amount of movement, become negative at the submerging area and afterwards reduce intensity, tending towards a value close to zero;
the simulation also confirms the rise in component v at the same time there is a drop in the value of component u. Also, the transversal velocity tends toward zero values when the current moves away from the reservoir entrance;
all of the simulated cases portray situations when the Reynolds and Richardson numbers tend to constant values when the longitudinal distance is increased, confirming the suppositions considered by other authors: Ellison and Turner (1959), Hauenstein and Dracos (1984), Alavian (1986), Akiyama et al. (1994), Tsihrintzis and Alavian (1996), Barbosa (1999).
the cases in which the entrances are extremely supercritical  experiments 21435 and 21436  show a better agreement (and also convergence) between the experiment and the simulation; the subcritical cases, experiments 21433 and 21434, are considerably more unstable and, sometimes, present simulated results that show fluctuations that do not exist in a real context. Such behavior is explained by the slight variation of the variables along each calculation station, i.e., the rises in x places the calculation points very close to each other, causing the truncation errors to be the same as the differences itself;
for the same reasons explained above, it is not advisable to reduce the increase in the calculation significantly for the subcritical cases of the previous item, because this may cause simulation instability and the impossibility of attaining results;
for the regions where the quantity of movement present considerable values, mainly close to the entrance, the simulated behavior of the density current width lies on a level that is far lower than what was experimentally observed. This happens because the equations used for the calculation do not take the evolution of the current spreading caused by the diffusive effects into account; however, what is observed in terms of width  the convective and diffusive effects  is not portrayed by the modeling;
based on what was explained in the previous item, the inclusion of such terms in the modeling for both calculation domains in future studies would be interesting;
in order to keep the same density current global features, there must be an interdependence of the variables. In the subcritical cases, it is possible to notice that the longitudinal velocities are generally underdimensioned. On the other hand, the thickness of the current usually becomes higher than the real value, so that certain compensation can be carried out.
As the contributions of this work it can be considered:
there are few studies published about evaluations of crosssectional component of velocity (ycomponent);
for the evaluation of zcomponent, this is the first study;
there is reasonable agreement between model and experimental results;
the possibility of modeling with changes in slope bottom of the reservoir (e.g. 4.5 to 12.5 degrees at position x = 150 cm).
Paper accepted February, 2010.
Technical Editor: Francisco Ricardo Cunha
 Akiyama, J., Ura, M., Wang, W., 1994, "Physicalbased numerical model of inclined starting plumes", Journal of Hydraulic Engineering, Vol. 120, No. 10, pp. 11391158.
 Alavian, V., 1986, "Behavior of density currents on an incline", Journal of Hydraulic Engineering, ASCE, Vol. 112, No. 1, pp. 2742.
 Alavian, V., Jirka, G.H., Denton, R.A., Johnson, M.C., Stefan, H.G, 1992, "Density currents entering lakes and reservoirs", Journal of Hydraulic Engineering, Vol. 118, No. 11, pp. 14641489.
 Barbosa, A.A., 1999, "Density currents in reservoirs (In Portuguese)", Ph.D. Thesis, University of São Paulo, São Carlos, 278 p.
 Birman, V.K., Meiburg, E., 2006, "Highresolution simulations of gravity currents", J. Braz. Soc. Mech. Sci. & Eng, Vol. 28, No. 2.
 Chow, V.T., 1959, "Open Channel Hydraulics", 1st Ed., McGrawHill Book Co., New York, N.Y., 680 p.
 Ellison, T.H., Turner, J.S., 1959, "Turbulent entrainment in stratified flows". Journal of Fluid Mechanics, Vol. 6, pp. 423448.
 Hauenstein, W., Dracos, T.H., 1984, "Investigation of plunging density currents generated by inflows in lakes". Journal of Hydraulic Research, IAHR, Vol. 22, No. 3, pp. 157179.
 Tsihrintzis, V.A., Alavian, V., 1996, "Spreading of threedimensional inclined gravity plumes", Journal of Hydraulic Research, IAHR, Vol. 34, No. 5, pp. 695711.
Publication Dates

Publication in this collection
01 Dec 2010 
Date of issue
Sept 2010