(1)
P. Brufau, P. Garcfa-Navarro
Fluid Mechanics. CPS. University of Zaragoza, Spain
M.E. Vazquez-Cendon
Applied
Mathematics, University of Santiago de Compostela, Spain A. Mendez Citeec,
University of Coruna, Spain
Abstract: Many efforts have been recently devoted to the development of multidimensional techniques for free surface flows. Among them, those oriented to the resolution of unsteady shallow water flow, have been strongly influenced by the upwind philosophy initially introduced in. These methods are specially applicable for advection dominated problems and their implementation is not straightforward when source terms are relevant. A first order finite volume method and its performance in the simulation of a laboratory experiment are presented.
There has been much research in CFD into the accurate and efficient solution of homogeneous systems of conservation laws. More recently, as numerical models become more complicated and the areas of application of these methods widens, it has become important that other aspects of the discretisation be given due attention. This is certainly true in the field of computational hydraulics where the modelling can be dominated by the effects not only of source terms, but also of wetting/drying fronts.
For many practical applications it is accepted that the unsteady flow of water in a one-dimensional approach is governed by the shallow water or St.Venant equations. These represent the conservation of mass and momentum along the direction of the main flow (Cunge et al., 1980). There exist methods developed to deal with Gas Dynamics problems (Euler equations) able to cope with complex systems of discontinuities and shock waves (Hirsch, 1990, MacCormack, 1971, Roe, 1981). It is well known that the nonlinearity of the hyperbolic shallow water equations may render their solution complicated, in the sense of giving rise to the appearance of discontinuities, which reflect physical phenomena such as hydraulic jumps and surges. The numerical technique adopted determines the quality of the results and second order classical schemes are known to show oscillatory behavior near these discontinuities. New schemes have been reported successful for simple flows but their adaptation to more complex situations is not straightforward. The presence of extreme slopes, high roughness or propagation over dry bed represent a great difficulty that can lead to ^important numerical errors (Bermudez and Vazquez-Cendon, 1994, Vazquez-Cendon, 1999, Brufau 2000)
In this work, an upwind scheme is used and validated with the results from a laboratory model of a 2D dam break. The model is applied to reproduce two field experiments illustrating the numerical method performance.
The equations will be solved by means of a finite volume technique based on an upwind scheme initially designed for systems of conservation laws. Therefore, we are interested in writing the 2D shallow overland flow equations in conservative form,
(1)
where U
represents the vector of conserved variables, h is depth of water, hu
and hv are the unit discharges along
the coordinate directions x, y
respectively with u and v the depth
averaged velocities and g is the
acceleration due to the gravity. F and
G are the fluxes of the conserved
variables across the edges of a control volume. They consist of the convective
fluxes together with the hydrostatic pressure gradients. The source terms
contain the effects of weight of fluid and bed roughness on the mass and
momentum conservation. S0x, S0y
are the bed slopes expressed in terms of the bottom depth and Sfx,
Sfy are the friction terms in the x;
y directions. For the friction term, the Manning equation has been used
expressing the energy grade line in terms of the Manning roughness coefficient n
(2)
A cell centered finite volume method is formulated for the equation over a triangular or quadrilateral control volume where the dependent variables of the system are represented as piecewise constants. The basic method will be first described for the pure conservation law, that is, without source terms. These will be next incorporated. The integral for a fixed cell area S,, applying the divergence theorem to the flux integral, is where C is the boundary of the area .S', and n is the outward unitary normal vector. The contour integral is approached via a mid-point rule and a numerical flux is defined at the midpoint of each edge. The evaluation of the numerical flux is based on the Riemann problem defined by the conditions on the left (L) and right (R) sides of the cell edges. An important feature of the ID Roe's approximate Riemann solver for non-linear systems of equations is exploited here. This is the definition of the approximated flux jacobian, constructed at the edges of the cells. The 2D numerical upwind flux is obtained by applying the expression in a 1 D form to each edge k of the computational cell along the normal- direction to the cell walls, so that represents the approximate Jacobian of the normal flux at that edge with the same shape as A but evaluated at an average state The approximate Jacobian matrix is not directly used in the actual method. Instead, the difference in the vector U across a grid edge is decomposed on the matrix eigenvectors basis as so that the matrix is replaced by its eigenvalues in the product. It has to be stressed at this point that in case of an advancing front over dry bed the average velocities are calculated in a different form (Brufau, 2000) because the velocity values at the right or left cell are zero. However the average value for the celerity is calculated always in the same form, otherwise the balance between the flux and the bed slope is not achieved in steady flow leading to numerical errors.
The source terms that account for the bed slopes are the only ones containing spatial derivatives. For this reason the discretization procedure will follow the flux term discretization as close as possible as suggested in (Bermudez et al, 1998).(refs.). An upwind approach has been adopted to model the bottom variations in order to ensure the best balance with the flux terms at least in steady cases. For the discrete representation, a numerical source H1*has to be defined as an approach of the surface integral over the cell of the source term H1. This numerical source will contain the contribution to cell i from all the edges surrounding the cell. At every cell-edge k the source term is discretized in an upwind manner. It is also decomposed into the eigenvector's basis in each edge so that
(3)
and represents the part of the discrete source term at edge k of cell ; associated to inward normal velocity. Source terms associated with friction were treated in a pointwise semi-implicit form.
Finally, the numerical scheme for every time step is formulated as
(4)
In order to avoid the numerical problems associated to extremely small values of water depth at the front, a procedure has been introduced so that all flow depths smaller than a certain user-defined threshold are ignored. In our case 10-6 m, is the minimum flow depth considered as part of the advancing front. In practical terms, every time the depth of flow is updated for a cell, the program checks whether it is smaller than this threshold. In that case it is reset to zero and this criterion is applied for both the advance and recession fronts.
A non-symmetric dam break laboratory problem is analyzed in a closed square pool (Mendez et al. 2001). The physical model (Fig. 1) was built in the CITEEC Hydraulic Laboratory (Coruna Spain). The pool is divided by a wall containing the dam in a reservoir and the plain to be flooded. The walls are solid and there is not any possibility of outflow from the domain. The bed is flat. Numerical and experimental results are presented with initial conditions:
h=0.5m in the reservoir and h=0.1m in the pool. In Fig. 2 initial conditions are specified in a free surface plot and contour levels at time t=3s are shown. Fig 3 shows a map of velocity vectors at time Is and 6s after the dam breaks for the latter case. In Fig. 4 comparison of the experimental data and numerical results is presented in the gauge points looking at the time evolution of the depth of water during 40s. In a second test case, the same model is used incorporating an obstacle in a side wall (pyramidal shaped). In Fig. 5 the measured time evolution of the depth of water at gauging points: S0, Sl, S2, S3, S4 and P18 and the numerical results are compared.. In all cases it can be noted a good agreement between experiment and numerical results. The gauging points situated inside the reservoir gave a good approximation of the outflow through the dam. The oscillations caused by the closed form of the model are well represented by the numerical model.
A finite volume based upwind scheme has been presented and applied to the numerical simulation of overland flows for irrigation purposes. The procedure is independent of the structured or unstructured character of the mesh. The flow problem is characterized by a thin
|T water layer and the movement is governed by the advection terms as well as by the source terms in the equations. The numerical treatment of both has been emphasized in the work presented. Strong bottom variations and high values of the roughness coefficient made necessary the application of more careful treatment of the corresponding terms in the equations to avoid numerical instabilities in the solution.
Acknowledgements
The authors would like to thank CICYT for providing the funding for this research under research project.
Abbott, M,B., Computational Hydraulics, Ashgate Pub. Comp., 1992. Bermudez A and M.E.Vazquez, Upwind methods for hyperbolic conservation laws with source terms, Computers and Fluids, 23(8)1049-1071, 1994. Bermudez.A, A.Dervieux, J-A.Desideri and M.E.Vazquez, Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes, Comput. Methods Appl.Mech.Engrg., 155 :49-72, 1998.
Brufau, P., Simulacion bidimensional de flujos hidrodinamicos transitorios en geometrias irregulares. PhD thesis, Universidad de Zaragoza, 2000. Chow.V.T, Open channel hydraulics, MacGraw-Hill Book Co. Inc.,1959. jCunge, J.A, F.M.Holly and A.Verwey, Practical aspects of computational river hydraulics, Pitman, 1980.
P.Glaister, Prediction of supercritical flow in open channels, Comput. Math. Applic., 24(7) :69-75, 1992.
Ch.Hirsch, Numerical computation of internal and external flows, John Wiley and Sons, New York, 1990.
RJ.LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm, J. Comput. Phys., 146(1) :346-365, 1998.
A. Mendez, L. Pena, P. Brufau and J. Puertas, An experimental approach to the Dam Break problem, XXIX IAHR Congress Proc., 2001.
P.L.Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys., 43(2) :357-372, 1981.
M.E. Vazquez-Cendon, Improved Treatment of Source Terms in Upwind Schemes for the Shallow Water Equations in Channels with Irregular Geometry, Journal of Computational Physics, 148, 497-526, 1999.
Fig. 1 Plane view of the physical model and gauging points
Fig. 2 Bottom and water depth initial conditions and contour levels at t=3s
Fig. 3 Velocity vectors for the case with obstacle
Fig.
4 Time evolution of the depth at points: P2, P6, P7 and P8. (No obstacle)
Fig. 5 Time evaluation of the dam break with obstacle at points: S0, S1, S2, S3, S4 and P18