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

## Abstract

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.

### Keywords

• 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

Qxdxdt+Tdxytdt=Qxdxdt+dxAtdt=0

Simplifying

Qx+Tyt=0E1

or

Qx+At=0E2

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

VAx+Tyt=0E3

or

AVx+VAx+Tyt=0E4

Because the hydraulic depth D = A/T and A=Ty, the above equation may be written as

DVx+Vyx+yt=0E5

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

qx+yt=0E6

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

γAΔhτ0PΔx=ρAΔxVVx+Vt

that is

τ0=γRhx+VgVx+1gVt
γRHx+1gVtE7

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

Hx+1gVt+V2C2R=0E8

and this equation may be rewritten as

Se+Sa+Sf=0E9

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.

By substituting H=V22g+y+zand the bed slope So=zxinto Eq. (8), we obtain

Hx=zx+yx+VgVx=So+yx+VgVx=1gVtSfE10

Hence Eq. (8) can be rewritten as

E11

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

Qxx+Axt=0E12
Qxt+Qx2/Axx+gAxHx+Sfx=0E13

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

Qx=1.486nAxR2/3Sfx1/2E14

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

mx=Qxt+Qx2/Axx/gAxE15

then Eq. (13) can be rewritten as

Sfx=Hx+mxE16

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

Qx=KxHx+mxE17

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

Kx=1.486nAxR2/3Hx+mx1/2E18

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

XKXHX+mX=AXtE19

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

XKXHX=AXtE20

where KXis now simplified as

Kx=1.486nAxR2/3HX12E21

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

XKXHX=WXHtE22

Assumptions other than mX= 0 in Eq. (19) result in a family of models:

mx={(Qx2/AX)X/gAX(convectiveaccelerationmodel)QXt/gAX(localaccelerationmodel)QXt+(Qx2/AX)X/gAX(fullydynamicmodel)0(DHM)E23

### 4.2 Two-dimensional diffusion hydrodynamic model

The set of (fully dynamic) 2D unsteady flow equations consists of one equation of continuity

qxx+qyy+Ht=0E24

and two equations of motion

qxt+xqx2h+yqxqyh+ghSfx+HX=0E25
qyt+yqy2h+yqxqyh+ghSfy+Hy=0E26

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

mZ+Sfz+HZ=0,z=x,yE27

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,

qZ=1.486nh5/3Sfz1/2,z=x,yE28

Eq. (28) can be rewritten in the general case as

qZ=KZHZKZmZ,z=x,yE29

where

KZ=1.486nh5/3HS+mS1/2,z=x,yE30

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

qZ=KZHZ,z=x,yE31

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

XKxHX+yKyHy=HtE32

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

xKxHx+yKyHy+S=HtE33

where

S=xKxmx+yKxmy

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

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.

For a square grid of width δ

qГE=KXГEHEHCδE34

where

Kx|ГE1.486nh5/3ГE/HEHCδcosθ1/2;HEHc0;HEHC<εE35

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

HCi+1=HCi+HCiE36

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.

compact PDF
Citations in RIS format
Citations in bibtex format

## More

© 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

### 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:

### Related Content

Next compact

#### Verification of Diffusion Hydrodynamic Model

By Theodore V. Hromadka II and Chung-Cheng Yen