Comparison of DHM Results for One- and Two-Dimensional Flows with Experimental and Numerical Data

Abstract

In this chapter, the performance of DHM for one- and two-dimensional flows is compared with the results of HEC-RAS, HEC-RAS 2D, WSPG, TUFLOW, Mike 21, and OpenFOAM models. The latter four models are currently widely used in industry, and benchmarking their data with DHM can shed more light on the reliability of DHM. As the results indicate, for applications which do not violate the assumptions made in DHM, the results are in agreement.

1. Introduction

Numerical modeling of free surface flow across real-life applications is gaining momentum. These model domains are characterized by thousands of computational cells, and the physical characteristics have varying complexities. Over the last three decades, with the growth of computational and visualization resources, multiple numerical models have been developed for solving the free surface flow equations across one, two, and three dimensions. While some of these models are available for free in the public domain, others are licensed by their respective vendors. Based on the assumptions used in these models, the complexity of flow equations can range from Bernoulli’s energy equation to three-dimensional unsteady Navier-Stokes equations. The models continue to evolve as the physics of flow is better understood, and the need for accurately predicting the flow variables across large spatial and temporal domains as their values is an important factor in the hydraulic design of structures and other related applications. Computational fluid dynamics (CFD) models that focus on solving the complete Navier-Stokes equations, rather than the energy equation or shallow water equations which are used in the hydraulic models, are also gaining popularity among the hydraulic modeling community.

The goal of this chapter is to evaluate the DHM results with a few industry-wide established software and experimental data to underscore the advantages and limitations in the models. To this end, we have chosen one critical application each from one-dimensional and two-dimensional flows.

2. One-dimensional application

2.1 Flows with hydraulic jump

Modeling flows with hydraulic jump where the flow transits from super critical to subcritical has been used by different researchers [1, 2, 3, 4] to test the reliability of their numerical formulations. Hydraulic jump is often created inflows to dissipate the flow energy, which can otherwise among others, erode the channels. They occur in gravity flows and are characterized by a large variation in flow depth and velocity before (Froude number > 1) and after (Froude number < 1) the jump. While capturing the internal flow details like bubble breakup, turbulence characteristics, tracking the water surface, aeration, fluid mixing, and turbulence is not possible using the shallow water equations, these equations can, however, predict the location and the flow depths before and after of the jump at steady state, which are important variables in the design calculations.

2.2 Experimental setup and model variables

A dataset from a series of experiments [5] that were conducted in the hydraulics laboratory at California State University, Fullerton, to simulate the location of steady-state hydraulic jump, was used for validating the models. The rectangular open channel flume was 15.2 m long, 0.46 m wide and 0.6 m in height. The channel sides are of glass, while the bottom interface with water is a metal sheet with a Manning roughness coefficient of 0.01. The bottom slope of the channel can be changed by tilting the flume, and in this investigation, it was set to 0.012. The flow discharge is 0.036 m³/s.

Boundary conditions need to be consistent with the physics of flow and appropriately complement the flow Equations [6]. The number of boundary conditions at the two ends of the flow domain is governed by the local Froude number. From a mathematical perspective, a boundary condition is a constraint imposed at the boundary node to arrive at a unique solution to a well-posed equation set. Specifying more or less than the required number may make the problem “ill-posed” and can lead to incorrect solutions. While one-dimensional supercritical flow requires superimposing two boundary conditions at the upstream end, a subcritical flow requires superimposing one boundary condition at the downstream end. In this simulation, at the upstream end, a flow depth of 0.04 m and flow discharge was specified. At the downstream end, a constant flow depth of 0.24 m was used. For this flow and boundary depth combination, the Froude numbers at the upstream and downstream end of the channel are 3.37 and 0.66, respectively.

2.3 Examined numerical models

The results of the DHM, RAS, WSPG, and TUFLOW models were compared with the experimental data. The other three models are briefly described below.

HEC-River Analysis System (RAS): HEC-RAS (steady state) model is based on the solution of the one-dimensional energy equation between two sections with energy losses given by Manning’s equation. The momentum equation can be used in situations where the water surface profile is rapidly varied as in hydraulic jump, hydraulics of bridges, and evaluating profiles at river confluences (stream junctions). RAS model also has modules to solve unsteady flows, sediment transport, and water quality analysis. In this work, the steady-state model was used. The model was developed by the US Army Corps of Engineers and can be downloaded for free [7].

Water Surface Pressure Gradient (WSPG): WSPG is one of the first models in computational hydraulics that was developed by the Los Angeles County Department of Public Works. The model solves the Bernoulli energy equation between any two cross sections, using the standard step method. The program computes uniform and nonuniform steady flow water surface profiles. As part of the solution, it can automatically identify a hydraulic jump in the channel reach. The model is currently distributed for a fee by civil design [8].

Two-dimensional unsteady flow (TUFLOW): The two-dimensional depth-averaged shallow water equations are solved in TUFLOW using a structured grid system with an alternating direction implicit scheme. The algorithm can capture flow transitions from supercritical to subcritical. TUFLOW incorporates the 1-D component (ESTRY software) or quasi-2D modeling system based on the full one-dimensional free surface flow Equations [9]. The model was developed by BMT WBM and can be downloaded for a fee.

2.4 Results

Figure 1 is a plot of the steady-state depth profile from the four models together with the experimental data. But for DHM, all other models satisfactorily predict the location and the flow depth before and after the jump. The reason as to why DHM could not capture the jump is because of the number of boundary conditions that the DHM permits from the end user. At the upstream end, DHM allows for only the flow discharge to be specified (and not two boundary conditions). The model had 18 grid elements in the computational domain. The upstream element is #1, and the end downstream element is #18. Element #14 corresponds to the end of the channel (length = 15.2 m). At element #1, the input discharge is 0.036 m³/s. At element #18, critical depth condition is specified, and the grid elevation was progressively raised such the depth at element #14 equals 0.24 m.

Figure 1.

Steady-state results for one-dimensional rapidly varying flow.

Because of this boundary condition limitation in the DHM, the jump is smoothened out in the solution. Although the downstream depth is consistent with other models, at the upstream end, the DHM predicted depth is higher than actual depth. DHM computed the flow transitioning from supercritical to subcritical without going through a hydraulic jump as required by theory and observed in the flume model. It can be concluded that DHM cannot be used in applications which require the prediction of the location of hydraulic jump.

3. Two-dimensional applications

3.1 Overland flow on a sloping domain

Overland flow is a dynamic response of the watershed to excess rainfall. Overland flow typically occurs as sheet flow on the land surface, and when the flow joins a channel, it is known as streamflow. The spatial and temporal distribution of two-dimensional overland flow variables is driven by the topographical characteristics of the domain and the boundary conditions. The nonhomogeneous surface characteristics and varying width of the natural watercourse path in the direction of flow make it an ideal case for comparing the results of different models.

Modeling overland flow has drawn the attention of many researchers. Since excess rainfall can cause flooding and mudslides which has the potential to cause loss of human life and disrupt the local economy, reliably predicting the flow variables for different precipitation scenarios can assist decision-makers and emergency personnel. Researchers have modeled these flows by solving the two-dimensional fully dynamic shall water Equations [10, 11] or their diffusion [12, 13] or kinematic approximations [14, 15].

3.1.1 Examined numerical models

The results of the DHM, extended DHM (EDHM), RAS, WSPG, and the CFD (OpenFOAM) models were compared. The back engine for EDHM is identical to DHM, with the primary difference between the two models being the larger array size of the variables in EDHM. When DHM was originally developed, the maximum array dimension was limited to 250, largely because of the available computational resources in 1980s. In EDHM, the dimension of all the arrays was increased to 9999. Since background information relating to RAS and WSPG has been given earlier, characteristic features of the CFD model, OpenFOAM, are summarized below.

OpenFOAM (Open-source Field Operation and Manipulation) is a freely available open-source software [16] that is gaining popularity across CFD applications. Its versatile C++ toolbox for the Linux operating system enables developing customized, efficient numerical solvers, and pre-/post-processing utilities for all kinds of CFD flows [17]. This code uses a tensorial approach following the widely known finite volume method (FVM), first used by McDonald [18]. Both structured and unstructured meshes can be used in the computational domain. The time integration can be done through backward Euler, steady-state solver, and Crank-Nicholson. The available gradient, divergence, Laplacian, and interpolation schemes are second-order central difference, fourth-order central difference, first-order upwind, and first-/second-order upwind. The turbulence models that can be used in OpenFOAM are LES, k-ɛ, and k-ɯ. The available solvers, options in specifying the boundary conditions, mesh generation tools, flow visualization software, and extensive documentation are making OpenFOAM popular among the CFD modeling community [2, 13, 19].

3.1.2 Study area and model variables

The study area is shown in Figure 2. Overland flow generated by a storm down a steep slope hits the flat main street after which a significant portion of the flow continues flowing North. The flow is supercritical from the inlet boundary location to the downstream of the main street. At the street downstream, there is a wall which reduces the flow velocity, thus forcing the flow to be critical. In this analysis, lateral flow on the main street was neglected.

Figure 2.

Map of the study area. The inflow and outflow boundary locations are identified by and .

In the models, the Manning roughness coefficient was set to 0.015 for the street portion and 0.03 for the earth. A uniform grid size of 15 ft was used, which resulted in a total of grids in the domain. The grids were oriented with the natural watercourse (NWC) path. The NWC ranged between 27 and 35 ft in the vicinity of the upstream end (southwest corner), and its width ranged between 45 and 60 ft in the vicinity where the water hits the main street. Having a 15 ft grid enabled us to cover the entire NWC path (Figure 3). The elevation at the center of the grid in the DHM was obtained from the topography map of the area. For the RAS and WSPG models, the required cross-sectional data was obtained from the top map. The rest of the input variables were consistent with the DHM data.

Figure 3.

DHM computational domain. The domain is aligned with the natural watercourse path and had 248 cells, which are 15 ft squares.

The intensity of rainfall and the bottom slope enhances the power of gravity-driven overland flow to make it supercritical (Froude number > 1). The available power in the water near the street can potentially push any moving or stationary automobiles. At the upstream end center cell, the flow hydrograph (Figure 4) was used as the boundary condition. The peak discharge at t = 0.5 hours is 755 cfs. At the downstream end, a critical depth boundary was specified. To keep the effect of downstream boundary minimal on the solution, the domain was extended by about 250 ft north of the main street.

Figure 4.

Input hydrograph at the upstream end.

Our focus was on predicting the flow depth at multiple probe locations on the main street. To conserve space, these results are plotted at two of the thirteen probe locations. These probe locations are 9 and 13 (Figure 5).

Figure 5.

Location of the probes along the center of the main street. The flow depth results were compared at probe locations 9 and 13.

3.1.3 Results

The DHM computational results include a variety of hydraulics relationships that are useful in further detailed analysis, such as flow velocity, flow depth, Froude number, and so forth. Of course, the code can be readily included in the DHM or as a post-processor routine, which enhances the DHM outcome. Of particular interest are the computational results from the DHM in comparison with the computational results produced by the CFD application. To display these computational results, hypothetical “probes” are inserted into the computational mesh where computational results are assembled and collated into a form suitable for visualization. The results from the visualization assessment are depicted below.

Figures 6 and 7 are plots comparing the flow depth, at probe location 9 and 13 on the main street. The data from DHM, EDHM, WSPG, RAS, and CFD (OpenFOAM) are plotted, for a specific CFD simulation time period. It is noted that the computational results are of high similarity, yet the computational effort ranges vastly. The DHM takes approximately 1 hour of CPU for the indicated simulation. In comparison, because of the varying small grid sizes in the domain, the CFD model required 2 weeks of CPU time using a parallel processor.

Figure 6.

Comparison of flow depth at probe 9 location.

Figure 7.

Comparison of flow depth at probe 13 location.

3.2 Open channel flow with a constriction

Numerically predicting the characteristics of flow through a channel with symmetric abrupt constriction (Figure 8) in the form of reduced channel width has drawn the attention of many researchers and has been part of any standard text book in hydraulics. The idea of developing DHM for this application was inspired after reading a recent paper [20] who tested a series of 2D models for multiple applications, one of which is flow through a constriction. Our focus was to estimate the DHM head loss at steady state and compare it with the published data.

Figure 8.

Definition sketch of the test problem along with the location of the two points (P1 and P2).

3.2.1 Examined numerical models

The results of the extended DHM (EDHM), Mike 21, TUFLOW, and HEC-RAS 2D models were compared, with the bench mark data from the equations provided by the Federal Highways Administration (FHA). Mike 21 and HEC-RAS 2D are briefly outlined here. Mike 21 solves the two-dimensional free surface flows where stratification can be neglected. It was originally developed for flow simulation in coastal areas, estuaries, and seas. The various modules of the system simulate hydrodynamics, advection-dispersion, short waves, sediment transport, water quality, eutrophication, and heavy metals [21].

HEC-RAS 2D (5.0.1) solves the two-dimensional Saint-Venant Equations [7] for shallow water flows using the full momentum computational method. The equations can model turbulence and Coriolis effects. For flow in sudden contraction, which is accompanied by high velocity, using the full momentum method in RAS 2D is recommended. The model uses an implicit finite volume solver.

3.2.2 Model variables

Figure 8 is the definition sketch of the test problem. The rectangular channel is 3100 ft long and 320 ft wide. The constriction is 60 × 60 ft. The channel length before constriction is 310 ft, and its length after constriction is 2730 ft. The computational domain in EDHM had 10 ft square cells, and the total number of cells was 9920. The longitudinal slope is 1%, the transverse slope is zero, and the model was run for a total of 1 hour. The upstream inflow is 1000 cfs. Since there are 30 cells at the upstream end, a uniform steady inflow of 33.3 cfs was specified at each of the cells. At the downstream end, a free overall boundary was specified. Constricting the flow area results in loss of energy. This loss of energy is reflected in a rise in energy gradient line and energy upstream of the constriction. Of interest is to estimate the head loss that occurs between points 1 and 2 (shown in Figure 8). The head loss (H_L) equals WSE₂ − WSE₁, where WSE is the water surface elevation [20].

3.2.3 Results

Table 1 shows the comparison of the head loss value obtained from the models along with the published data of other models. It is noted that in this effort, the computational models are compared with respect to head loss (as given in [20]) through the constriction, and this is the primary form of assessment. The DHM WSE change value is within the range of other model predictions, although all the model predictions are above the FHA value.

	WSE change (ft)
HEC-RAS 2D	1.27
TUFLOW	0.99
Mike 21	1.28
EDHM	1.16
FHA equation	0.8

Table 1.

Comparison of change in water surface elevation at constriction between EDHM and published data*.

^*

Except for EDHM all other data were obtained from the literature [20].

4. Conclusions

Results from multiple computer models are compared with those of DHM for one- and two-dimensional flows. The considered one-dimensional flow was a mixed flow with a hydraulic jump. All the model results (DHM, WSPG, RAS, TUFLOW) were compared with the benchmark experimental data. Because of the way the boundary conditions are specified in the DHM, the model cannot simulate the hydraulic jump. For the two-dimensional overland flow, the model results (DHM, EDHM, TUFLOW, MIKE 21, WSPG, RAS, RAS2D, and the CFD model, OpenFOAM) were compared between themselves. The agreement of the predicted flow variables reinforces the reliability of the current model.