In this chapter, the governing flow equations for one- and two-dimensional unsteady flows that are solved in the diffusion hydrodynamic model (DHM) are presented along with the relevant assumptions. A step-by-step derivation of the simplified equations which are based on continuity and momentum principles are detailed. Characteristic features of the explicit DHM numerical algorithm are discussed.
- unsteady flow
- conservation of mass
- finite difference
- explicit scheme
- flow equations
Many flow phenomena of great engineering importance are unsteady in characters and cannot be reduced to a steady flow by changing the viewpoint of the observer. A complete theory of unsteady flow is therefore required and will be reviewed in this section. The equations of motion are not solvable in the most general case, but approximations and numerical methods can be developed which yield solutions of satisfactory accuracy.
2. Review of governing equations
The law of continuity for unsteady flow may be established by considering the conservation of mass in an infinitesimal space between two channel sections (Figure 1). In unsteady flow, the discharge, Q, changes with distance, x, at a rate , and the depth, y, changes with time, t, at a rate . The change in discharge volume through space dx in the time dt is . The corresponding change in channel storage in space is in which . Because water is incompressible, the net change in discharge plus the change in storage should be zero, that is
At a given section, Q = VA; thus Eq. (1) becomes
Because the hydraulic depth D = A/T and , the above equation may be written as
The above equations are all forms of the continuity equation for unsteady flow in open channels. For a rectangular channel or a channel of infinite width, Eq. (1) may be written as
where q is the discharge per unit width.
3. Equation of motion
In a steady, uniform flow, the gradient, , of the total energy line is equal to magnitude of the “friction slope” , where C is the Chezy coefficient and R is the hydraulic radius. Indeed this statement was in a sense taken as the definition of Sf; however, in the present context, we have to consider the more general case in which the flow is nonuniform, and the velocity may be changing in the downstream direction. The net force, shear force and pressure force, is no longer zero since the flow is accelerating. Therefore, the equation of motion becomes
where τ0 is the same shear stress, P is the hydrostatic pressure, h is the depth of water, Δh is the change of depth of water, γ is the specific weight of the fluid, R is the mean hydraulic radius, and ρ is the fluid density. Substituting = into Eq. (7), we obtain
and this equation may be rewritten as
By substituting and the bed slope into Eq. (8), we obtain
Hence Eq. (8) can be rewritten as
This equation may be applicable to various types of flow as indicated. This arrangement shows how the nonuniformity and unsteadiness of flows introduce extra terms into the governing dynamic equation.
4. Diffusion hydrodynamic model
4.1 One-dimensional diffusion hydrodynamic model
where Qx is the flow rate; x,t are spatial and temporal coordinates, Ax is the flow area, g is the gravitational acceleration, H is the water surface elevation, and Sfx is a friction slope. It is assumed that Sfx approximated from Manning’s equation for steady flow by .
where R is the hydraulic radius and n is a flow resistance coefficient which may be increased to account for other energy losses such as expansions and bend losses.
Letting mx be a momentum quantity defined by
then Eq. (13) can be rewritten as
In Eq. (15), the subscript x included in mx indicates the directional term. The expansion of Eq. (13) to two-dimensional case leads directly to the terms (mx, my) except that now a cross-product of flow velocities is included, increasing the computational effort considerably.
where Qx indicates a directional term and Kx is a type of conduction parameter defined by
In Eq. (18), is limited in value by the denominator term being checked for a smallest allowable magnitude (such as > ).
The one-dimensional model of Akan and Yen  assumed = 0 in Eq. (18). The term is assumed to be negligible when combined with the other similar terms—that is, they are considered as a sum rather than as individual directional terms that typically have more significance when examined individually. Additionally, the term “diffusion” routing indicates assuming that several convective and other components have a small contribution to the coupled mass and energy balance equations and therefore are neglected in the computational formulation to simplify the model accordingly. Thus, the one-dimensional DHM equation is given by
where is now simplified as
For a channel of constant width, , Eq. (20) reduces to
Assumptions other than = 0 in Eq. (19) result in a family of models:
4.2 Two-dimensional diffusion hydrodynamic model
The set of (fully dynamic) 2D unsteady flow equations consists of one equation of continuity
and two equations of motion
where and are flow rates per unit width in the x and y directions; and represents friction slopes in x and y directions; H, h, and g stand for water surface elevation, flow depth, and gravitational acceleration, respectively; and x, y, and t are spatial and temporal coordinates.
The above equation set is based on the assumptions of constant fluid density without sources or sinks in the flow field and of hydrostatic pressure distributions.
where represents the sum of the first three terms in Eqs. (25) or (26) divided by gh. Assuming the friction slope to be approximated by the Manning’s formula, one obtains, in the US customary units for flow in the x or y directions,
Eq. (28) can be rewritten in the general case as
The symbol s in Eq. (30) indicates the flow direction which makes an angle of with the positive x direction.
The term is assumed to be negligible [1, 2, 3, 4, 5] when combined with the other similar terms, i.e., they are considered as a sum rather than as individual directional terms that typically have more significance when examined individually. Additionally, the term “diffusion” routing indicates assuming that several convective and other components have a small contribution to the coupled mass and energy balance equations and therefore are neglected in the computational formulation to simplify the model accordingly. Neglecting this term results in the simple diffusion model
If the momentum term groupings were retained, Eq. (32) would be written as
and and are also functions of and , respectively.
5. Numerical approximation
5.1 Numerical solution algorithm
The one-dimensional domain is discretized across uniformly spaced nodal points, and at each of these points, at time (t) = 0, the values of Manning’s n, an elevation, and initial flow depth (usually zero) are assigned. With these initial conditions, the solution is advanced to the next time step (t + ) as detailed below
Between nodal points, compute an average Manning’s n and average geometric factors
Using the flow depths at time t and t + , estimate the mid time step value of selected from Eq. (23)
Recalculate the conductivities using the appropriate values
Determine the new nodal flow depths at the time (t + ) using Eq. (19), and
Return to step (3) until matches mid time step estimates.
The above algorithm steps can be used regardless of the choice of definition for from Eq. (23). Additionally, the above program steps can be directly applied to a two-dimensional diffusion model with the selected (, ) relations incorporated.
5.2 Numerical model formulation (grid element)
For uniform grid elements, the integrated finite difference version of the nodal domain integration (NDI) method  is used. For grid elements, the NDI nodal equation is based on the usual nodal system shown in Figure 3. Flow rates across the boundary Г are estimated by assuming a linear trial function between nodal points.
For a square grid of width δ
In Eq. (35), h (depth of water) and n (the Manning’s coefficient) are both the average of their respective values at C and E, i.e., and . Additionally, the denominator of is checked such that is set to zero if is less than a tolerance such as ft.
The net volume of water in each grid element between time step i and i + 1 is and the change of depth of water is for time step i and i + 1 with interval. Then the model advances in time by an explicit approach
where the assumed input flood flows are added to the specified input nodes at each time step. After each time step, the hydraulic conductivity parameters of Eq. (35) are reevaluated, and the solution of Eq. (36) is reinitiated.
5.3 Model time step selection
The sensitivity of the model to time step selection is dependent upon the slope of the discharge hydrograph () and the grid spacing. Increasing the grid spacing size introduces additional water storage to a corresponding increase in nodal point flood depth values. Similarly, a decrease in time step size allows a refined calculation of inflow and outflow values and a smoother variation in nodal point flood depths with respect to time. The computer algorithm may self-select a time step by increments of halving (or doubling) the initial user-chosen time step size so that a proper balance of inflow-outflow to control volume storage variation is achieved. In order to avoid a matrix solution for flood depths, an explicit time stepping algorithm is used to solve for the time derivative term. For large time steps or a rapid variation in the dam-break hydrograph (such as is large), a large accumulation of flow volume will occur at the most upstream nodal point. That is, at the dam-break reservoir nodal point, the lag in outflow from the control volume can cause an unacceptable error in the computation of the flood depth. One method that offsets this error is the program to self-select the time step until the difference in the rate of volume accumulation is within a specified tolerance.
Due to the form of the DHM in Eq. (22), the model can be extended into an implicit technique. However, this extension would require a matrix solution process which may become unmanageable for two-dimensional models which utilize hundreds of nodal points.
The one- and two-dimensional flow equations used in the diffusion hydrodynamic model are derived, and the relevant assumptions are listed. These equations, which are the basis of the model, are based on the conservation of mass and momentum principles. The explicit numerical algorithm and the discretized equations are also presented. The ability of the model to self-select the optimal time step is discussed.
Akan AO, Yen BC. Diffusion-wave flood routing in channel networks. ASCE Journal of Hydraulic Division. 1981; 107(6):719-732
Hromadka TV II, Berenbrokc CE, Freckleton JR, Guymon GL. A two-dimensional diffusion dam-break model. Advances in Water Resources. 1985; 8:7-14
Xanthopoulos TH, Koutitas CH. Numerical simulation of a two-dimensional flood wave propagation due to dam failure. ASCE Journal of Hydraulic Research. 1976; 14(4):321-331
Hromadka TV II, Lai C. Solving the two-dimensional diffusion flow model. In: Proceedings of ASCE Hydraulics Division Specialty Conference, Orlando, Florida. 1985
Hromadka TV II, Nestlinger AB. Using a two-dimensional diffusional dam-break model in engineering planning. In: Proceedings of ASCE Workshop on Urban Hydrology and Stormwater Management, Los Angeles County Flood Control District Office, Los Angeles, California. 1985
Hromadka TV II, Guymon GL, Pardoen G. Nodal domain integration model of unsaturated two-dimensional soil water flow: Development. Water Resources Research. 1981; 17:1425-1430