Open access peer-reviewed compact

Diffusion Hydrodynamic Model Theoretical Development

By Theodore V. Hromadka II and Chung-Cheng Yen

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

DOI: 10.5772/intechopen.93207

Downloaded: 398


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=dxAtdtin 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= V2C2Rinto 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+zand the bed slope So=zxinto 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), Kxis 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 mXterm 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 KXis 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 qXand qyare flow rates per unit width in the x and y directions; Sfxand Sfyrepresents 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 mZrepresents 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/qxwith the positive x direction.

The mzterm 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 Kxand Kyare also functions of mxand 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 mXselected from Eq. (23)

  4. Recalculate the conductivities KXusing the appropriate mXvalues

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

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

The above algorithm steps can be used regardless of the choice of definition for mXfrom 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/2and n=nC+nE/2. Additionally, the denominator of KXis checked such that KXis set to zero if HEHCis less than a tolerance εsuch as 103ft.

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

© 2020 The Author(s). Licensee IntechOpen. This compact is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited.

How to cite and reference

Link to this compact Copy to clipboard

Cite this compact Copy to clipboard

Theodore V. Hromadka II and Chung-Cheng Yen (September 9th 2020). Diffusion Hydrodynamic Model Theoretical Development, A Diffusion Hydrodynamic Model, Theodore V. Hromadka II, Chung-Cheng Yen and Prasada Rao, IntechOpen, DOI: 10.5772/intechopen.93207. Available from:

compact statistics

398total compact downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next compact

Verification of Diffusion Hydrodynamic Model

By Theodore V. Hromadka II and Chung-Cheng Yen

Related Book

First chapter

Strategies for Testing the Impact of Natural Flood Risk Management Measures

By Barry Hankin, Peter Metcalfe, David Johnson, Nick A. Chappell, Trevor Page, Iain Craigen, Rob Lamb and Keith Beven

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us