Open access peer-reviewed chapter

Diffusion Hydrodynamic Model Theoretical Development

Written By

Theodore V. Hromadka II and Chung-Cheng Yen

Reviewed: June 17th, 2020 Published: September 9th, 2020

DOI: 10.5772/intechopen.93207

From the Compact

A Diffusion Hydrodynamic Model

Authored by Theodore V. Hromadka II, Chung-Cheng Yen and Prasada Rao

Chapter metrics overview

780 Chapter Downloads

View Full Metrics


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

1. Introduction

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 Qx, and the depth, y, changes with time, t, at a rate yt. The change in discharge volume through space dx in the time dt is Qxdxdt. The corresponding change in channel storage in space is Tdxytdt=dxAtdt in which A=Ty. Because water is incompressible, the net change in discharge plus the change in storage should be zero, that is

Figure 1.

Continuity of unsteady flow.






At a given section, Q = VA; thus Eq. (1) becomes




Because the hydraulic depth D = A/T and A=Ty, 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, dHdx, of the total energy line is equal to magnitude of the “friction slope” Sf=V2/C2R, 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


that is


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 τ0γR = V2C2R into Eq. (7), we obtain


and this equation may be rewritten as


where the three terms of Eq. (9) are called the energy slope, the acceleration slope, and the friction slope, respectively. Figure 2 depicts the simplified representation of energy in unsteady flow.

Figure 2.

Simplified representation of energy in unsteady flow.

By substituting H=V22g+y+z and the bed slope So=zx 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

The mathematical relationships in a one-dimensional diffusion hydrodynamic model (DHM) are based upon the flow equations of continuity (2) and momentum (11) which can be rewritten [1] as


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 [1].


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.

Rewriting Eq. (14) and including Eqs. (15) and (16), the directional flow rate is computed by


where Qx indicates a directional term and Kx is a type of conduction parameter defined by


In Eq. (18), Kx is limited in value by the denominator term being checked for a smallest allowable magnitude (such as HX+mX1/2 > 103).

Substituting the flow rate formulation of Eq. (17) into Eq. (12) gives a diffusion type of relationship


The one-dimensional model of Akan and Yen [1] assumed mX = 0 in Eq. (18). The mX 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 KX is now simplified as


For a channel of constant width, WX, Eq. (20) reduces to


Assumptions other than mX = 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 qX and qy are flow rates per unit width in the x and y directions; Sfx and Sfy 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.

The local and convective acceleration terms can be grouped, and Eqs. (25) and (26) are rewritten as


where mZ 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 θ=tan1qy/qx with the positive x direction.

The mz 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


The proposed 2D DHM is formulated by substituting Eq. (31) into Eq. (24)


If the momentum term groupings were retained, Eq. (32) would be written as




and Kx and Ky are also functions of mx and my, 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 + t) as detailed below

  1. Between nodal points, compute an average Manning’s n and average geometric factors

  2. Assuming mX = 0, estimate the nodal flow depths for the next time step (t + ∆t) by using Eqs. (20) and (21) explicitly

  3. Using the flow depths at time t and t + t, estimate the mid time step value of mX selected from Eq. (23)

  4. Recalculate the conductivities KX using the appropriate mX values

  5. Determine the new nodal flow depths at the time (t + t) using Eq. (19), and

  6. Return to step (3) until KX matches mid time step estimates.

The above algorithm steps can be used regardless of the choice of definition for mX from Eq. (23). Additionally, the above program steps can be directly applied to a two-dimensional diffusion model with the selected (mX, my) 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 [6] 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.

Figure 3.

Two-dimensional finite difference analog.

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., h=hC+hE/2 and n=nC+nE/2. Additionally, the denominator of KX is checked such that KX is set to zero if HEHC is less than a tolerance ε such as 103 ft.

The net volume of water in each grid element between time step i and i + 1 is qCi=qrE+qrw+qrN+qrS and the change of depth of water is HCi=qCit/δ2 for time step i and i + 1 with t 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 (Qt) 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 Qt 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.


6. Conclusions

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.


  1. 1. Akan AO, Yen BC. Diffusion-wave flood routing in channel networks. ASCE Journal of Hydraulic Division. 1981;107(6):719-732
  2. 2. Hromadka TV II, Berenbrokc CE, Freckleton JR, Guymon GL. A two-dimensional diffusion dam-break model. Advances in Water Resources. 1985;8:7-14
  3. 3. 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
  4. 4. Hromadka TV II, Lai C. Solving the two-dimensional diffusion flow model. In: Proceedings of ASCE Hydraulics Division Specialty Conference, Orlando, Florida. 1985
  5. 5. 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
  6. 6. 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

Written By

Theodore V. Hromadka II and Chung-Cheng Yen

Reviewed: June 17th, 2020 Published: September 9th, 2020