A Short Review on Computational Hydraulics in the Context of Water Resources Engineering

Shiblu Sarker

Department of Civil, Environmental, & Construction Engineering, University of Central Florida, Florida, USA.

DOI: 10.4236/ojmsi.2022.101001                

Abstract

The term “hydraulics” is concerned with the conveyance of water that can consist of very simple processes to complex physical processes, such as flow in open rivers, flow in pipes, the flow of nutrients/sediments, the flow of groundwater to sea waves. The study of hydraulics is primarily a mixture of theory and experiments. Computational hydraulics is very helpful to quantify and predict flow nature and behavior. The mathematical model is the backbone of the computational hydraulics that consists of simple to complex mathematical equations with linear and/or non-linear terms and ordinary or partial differential equations. Analytical solution to these mathematical equations is not feasible in the majority of cases. In these consequences, mathematical models are solved using different numerical techniques and associated schemes. In this manuscript, we aim to review hydraulic principles along with their mathematical equations. Then we aim to learn some commonly used numerical techniques to solve different types of differential equations related to hydraulics. Among them, the Finite Difference Method (FDM), Finite Element Method (FEM) and Finite Volume Method (FVM) have been discussed along with their use in real-life applications in the context of water resources engineering.

Keywords

Computational Hydraulics, Finite Difference Method, Finite Element Method, Finite Volume Method

1. Introduction

Hydraulics is a subfield of fluid mechanics that studies the behavior of water at rest and in motion. This term is frequently used in connection with hydrodynamics to refer to the forces acting between water and its boundaries and their effects on subsequent flow patterns. To measure and predict flow behavior, it is critical to apply physical conservation rules, which will be discussed in detail in the subsequent section. Although they attempted to solve numerous fluid motion problems strictly theoretically, assuming an ideal fluid (frictionless or non-viscous), their results were still extremely effective in specific scenarios and conditions. On the other hand, hydraulics developed mathematical formulas to handle real-world fluid flow problems using experimental data from numerous laboratory experiments and field observations. Ludwig Prandtl’s essential idea of boundary layer combined hydrodynamics and hydraulics. Additionally, the illustrious research of Reynolds, Froude, Prandtl, and von Karman persuades us that studying fluids requires a combination of theory and experiment. Here, we can trace the evolution of hydraulics all the way back to 200 BCE [1] (Figure 1).

Hydraulics originated around the time of Archimedes, who established the now-famous rules of floating bodies. Following Galileo (1564-1642), Gastelli (1577-1644), Torricelli (1600-1647), and Guglielmini (1655-1710) developed mathematical models for barometers, flow from containers, steady flow, and resistance in rivers. When Newton (1642-1727) published his famed laws of motion, they expanded our quantitative understanding of fluid resistance in terms of velocity gradient and drag on spheres. Following Newton, Daniel Bernoulli, Leonhard Euler, Clairaut, d’Alembert, Lagrange (1736-1813), Laplace (1749-1827), and Gerstner (1756-1832) contributed significant concepts on fluid flow and surface waves to hydrodynamics. Additionally, De Pitot invented a tube for measuring speed; Chezy devised a formula for open channel resistance; Borda conducted orifice-related experiments; Bossut installed a towing tank; and Venturi

Figure 1. Contributors in hydrodynamics and hydraulics (sources: Google image).

experimented with the flow in varying cross-sections, all of which contributed to the field of hydraulics’ enrichment [1].

Coulomb (1736-1806) conducted flow resistance studies concurrently with Ernst Brothers (1795-1878) and Weber (1804-1891), all of whom made significant contributions to the field of hydraulics. As a result, the Buron (1790-1873), Fourneyman (1802-1867), Coriolis (1792-1843), Francis (1815-1892), and Russel (1808-1882) all made significant contributions to engineering. The Hagen (1797-1889), Poiseuille (1799-1869), and Sazon Weisbach (1806-1871) all contribute significantly to pipe flow, whereas Saint-Venant (1797-1886) contributes to open channel flow. On the other hand, Dupuit (1804-1866), Bresse (1822-1883), Basin (1829-1917), Darcy (1803-1858), and Manning (1816-1897) conducted research on both open channel hydraulics and pipe flow. William Froude (1810-1879) and Robert Froude (1846-1924) in particular provide a valuable criterion for classifying the flow based on their ship model tests. Additionally, Osborne Reynolds did a successful experiment based on streamline analysis to discern between the laminar and turbulent flow. Navier (1785-1876), Cauchy (1789-1857), Poisson (1781-1840), Saint-Venant (1797-1886), Boussinesq (1842-1929), Stokes (1819-1903), Lord Rayleigh (1842-1919), Lamb (1849-1934), Helmholtz (1821-1894), and Kirchoff (1824-1887) all made significant contributions to the development of theoretical and applied hydrodynamics during the nineteenth century. Additionally, Euler’s equation of motion for an ideal (non-viscous) fluid had advanced significantly. However, it does not account for certain significant observed results, such as the decrease in pipe pressure. As a result, engineers in practice created their own empirical hydraulics. Recent computational hydraulics is one of several domains of science in which computer technologies enable an intermediate mode of operation between pure theory and experiment. This discipline, which is frequently referred to as computational hydraulics, is more of a synthesis of hydraulics and hydrodynamics than an independent development. The objective of computational hydraulics is to use computers to mimic various physical processes involving seas, estuaries, rivers, canals, and reservoirs. We will explore fundamental hydraulic processes and their quantification in this article, focusing on water resource applications such as open channels, pipes, groundwater, and coastal waves [1].

2. Hydraulics of Open Channels

Open channel flows deal with the flow of water in open channels, where the pressure at the surface of the water is ambient and the pressure at each portion is proportional to the depth of the water. Pressure head is the ratio of pressure and the specific weight of water ( $P/\gamma$ ). Elevation head (Z) is the height of the section under consideration above a datum and Velocity head ( $V^2/2g$ ) is attributed to the average velocity of flow in that vertical section. Hence, the total head can be expressed by the Equation (1).

$\mathcal{H}=Z+\frac{P}{\gamma }+\frac{V^2}{2g}$ (1)

In an open channel, water flow is primarily due to the gradient of the head and gravity. The Open Channels are primarily used for irrigation supply, industrial and domestic water supply [2].

2.1. Conservation Laws

Mass Conservation, Momentum Conservation and Energy Conservation are the main conservation laws used in open Channel. A control volume commonly considered under the Conservation of Mass, the overall mass change in the control volume due to inflow and outflow, is equal to the net rate of weight change in the control volume. This results to the classical continuity equation balancing the inflow, out flow and the storage change in the control volume. Since we only find water that is regarded as incompressible, that is, it is possible to disregard the impact on density [2].

$AV=Q$ (2)

The rate of change of momentum in the control volume under the Conservation of Momentum is equal to the net forces acting on the control volume. Since the water under consideration is flowing, external forces are acting upon it. This leads ultimately to the Newton’s second law of motion [2].

$\frac{Q^2}{gA}+\bar{y}A=Constant$ (3)

Energy conservation states that neither energy can be generated nor destroyed, it only changes its form. The energy would be in the form of potential energy and kinetic energy, primarily in open channels. The elevation of the water parcel is due to potential energy, while the kinetic energy is due to the movement of the parcel. In the context of open channel flow the total energy due these factors between any two sections is conserved, that leads to the classical Bernoulli’s equation.

$Z+y+\frac{V^2}{2g}=Constant$ (4)

This equation must account for the energy loss between the two sections when used between two sections, which is due to the resistance of bed shear etc. to the flow [2].

2.2. Types of Open Channel Flows

The flow in an open channel is classified as Sub critical flow, Super Critical flow, and Critical flow depending on the Froude number ( $F_r$ ), where the Froude number can be described as 5:

$F_r=\frac{V}{\sqrt{gy}}$ (5)

Types of open channel flow can be further categorized on the basis of time and space parameters (see Figure 2(a)). Flow is said to be steady when discharge

Figure 2. (a) Types of open channel flows; (b) Specific energy curve for a given discharge.

does not change along the course of the channel flow. On the other hand, as the discharge varies over time, flow is said to be unsteady. When both the depth and discharge are the same at any two sections of the channel, flow is said to be uniform. Furthermore, flow is said to be gradually varied whenever the depth changes gradually along the channel, whenever the flow depth changes rapidly along the channel the flow is termed rapidly varied flow. Whenever the flow depth gradually varies due to the change in discharge, the flow is called a spatially varied flow [2]. Steady uniform flow (kinematic wave), steady non-uniform flow (diffusion wave) and unsteady non-uniform flow (dynamic wave) can be seen as types of possible flow based on the above classification.

2.3. Specific Energy and Specific Force

Specific energy is defined as the energy acquired by a section of water due to its depth and the velocity with which it is flowing, specific Energy E is given by,

$E=y+\frac{V^2}{2g}$ (6)

where y is the depth of flow at that section and v is the average velocity of flow. Specific energy is minimum at critical condition (see point C in Figure 2), and the corresponding depth at point C is termed as critical depth y_c. This y_c can be measured for a specific section under steady flow [2].

Specific force is defined as the sum of the momentum of the flow passing through the channel section per unit time per unit weight of water and the force per unit weight of water [2].

$F=\frac{Q^2}{gA}+yA$ (7)

The specific forces of two sections are equal provided that the external forces and the weight effect of water in the reach between the two sections can be ignored. The term specific force is very useful to quantify the property of rapidly varied flow such as hydraulic jump [2]. At the critical state of flow the specific force is a minimum for the given discharge. When the specific energy is minimal, flow is critical. Also, the flow has to go through critical conditions if the flow changes from sub critical to super critical or vice versa. Sub-critical flow-the depth of the flow is higher, although the velocity is lower and the super-critical flow-the depth of the flow is lower, but the velocity is higher. Critical flow is the

flow over a free over-fall. For Specific energy to be a minimum $\frac{\text{d}E}{\text{d}y}=0$ :

$\frac{\text{d}E}{\text{d}y}=1-\frac{Q^2}{g{A}^3}\frac{\text{d}A}{\text{d}y}$ (8)

However, $\text{d}A=T\ast \text{d}y$, where T is the width of the channel at the water surface, then applying $\frac{\text{d}E}{\text{d}y}=0$, will result in following:

$\frac{Q^2{T}_c}{g{A}_c^3}=1;\frac{A_c}{T_c}=\frac{Q^2}{g{A}_c^2};\frac{A_c}{T_c}=\frac{V_c^2}{g}$ (9)

For a rectangular channel $\frac{A_c}{T_c}=y_c$. Following the derivation for a rectangular channel,

$F_r=\frac{V_c}{\sqrt{g{y}_c}}=1$ (10)

The same notion applies to trapezoidal and other cross-sections [2]. The critical flow state describes an unique relationship between depth and discharge that is very useful in the design of flow measurement structures.

2.4. Uniform Flows

This is one of the most important principles of open channel flow. The most important uniform flow equation is Manning’s equation stated as:

$V=\frac{1}{n}\ast R^{2/3}\ast S^{1/2}$ (11)

where R = the hydraulic radius = A/p and p = wetted perimeter = $f\left(y\mathrm{,}{S}_0\right)$, y = depth of the channel bed, S₀ = bed slope (same as the energy slope, S_f) and n = the Manning’s dimensional empirical constant. Uniform Flow concept is used in most of the open channel flow design. The uniform flow implies that there is no acceleration of the flow due to the weight portion of the flow being balanced by the resistance of the shear of the bed. In terms of discharge the Manning’s Equation (11) is given by:

$Q=\frac{1}{n}\ast A\ast R^{2/3}\ast S^{1/2}$ (12)

This is a non linear equation in y the depth of flow for which most of the computations will be made. Derivation of uniform flow equation is given below, where $W\sin \theta$ = weight component of the fluid mass in the direction of flow, ${\tau }_0$ = bed shear stress and $P\Delta x$ = surface area of the channel. The force balance equation can be written as:

$W\sin \theta -{\tau }_{0}p\Delta x=0\to \gamma A\Delta x\sin \theta -{\tau }_{0}p\Delta x=0\to {\tau }_0=\gamma \frac{A}{p}\sin \theta$ (13)

NowA/p is the hydraulic radius,R, and $\sin \theta$ is the slope of the channel S_o. The shear stress can be expressed as:

${\tau }_0=\rho C_f\frac{V^2}{2}$ (14)

where $c_f$ is resistance coefficient, V is the mean velocity $\rho$ is the mass density. Therefore the Equation (14) can be written as:

$\rho C_f\frac{V^2}{2}=\gamma R{S}_o\to V=\sqrt{\frac{2g}{C_f}}\sqrt{R{S}_o}\to V=\mathcal{C}\sqrt{R{S}_o}$ (15)

where C is Chezy’s constant. For Manning’s equation,

$\mathcal{C}=\frac{1.49}{n}{R}^{\frac{1}{6}}$ (16)

2.5. Gradually Varied Flow

Flow is said to be gradually varied whenever the depth of flow changed gradually. The governing equation for gradually varied flow is given by:

$\frac{\text{d}y}{\text{d}x}=\frac{S_o-S_f}{1-F_r^2}$ (17)

where the variation of depth y with the channel distance x is shown to be a function of bed slope $S_o$, Friction Slope $S_f$ and the flow Froude number $F_r$. This is a non linear equation with the depth varying as a non linear function (Figure 3).

Figure 3. (a) Steady uniform flow in an open channel; (b) Total head at a channel section.

Gradually varied flow can be derived from the conservation of energy at two sections of a reach of length $\Delta x$, can be written as:

$y_1+\frac{V_1^2}{2g}+S_o\Delta x=y_2+\frac{V_2^2}{2g}+S_f\Delta x$ (18)

Now, let $\Delta y=y_2-y_1$ and $\frac{V_2^2}{2g}-\frac{V_1^2}{2g}=\frac{\text{d}}{\text{d}x}\left(\frac{v^2}{2g}\right)\Delta x$,

Then the above equation becomes:

$\Delta y=S_o\Delta x-S_f\Delta x-\frac{\text{d}}{\text{d}x}\left(\frac{v^2}{2g}\right)\Delta x$ (19)

Dividing through $\Delta x$ and taking the limit as $\Delta x$, approaches zero gives us:

$\frac{\text{d}y}{\text{d}x}+\frac{\text{d}}{\text{d}x}\left(\frac{v^2}{2g}\right)=S_o-S_f$ (20)

After simplification, can be done in terms of Froude number and the general differential equation can be written as Equations (21) to (23):

$\frac{\text{d}y}{\text{d}x}=\frac{S_o-S_f}{1+\frac{\text{d}}{\text{d}y}\left(\frac{V^2}{2g}\right)}$ (21)

$\frac{\text{d}}{\text{d}y}\left(\frac{V^2}{2g}\right)=\frac{\text{d}}{\text{d}y}\left(\frac{Q^2}{2g{A}^2}\right)=\frac{-2Q^2}{2g{A}^3}\left(\frac{\text{d}A}{\text{d}y}\right)=\frac{-Q^2}{g{A}^2}\left(\frac{1}{D}\right)=-F_r^2$ (22)

$\frac{\text{d}y}{\text{d}x}=\frac{S_o-S_f}{1-F_r^2}$ (23)

Numerical integration of the gradually varied flow equation will give the water surface profile along the channel. Depending on the depth of flow where it lies when compared with the normal depth and the critical depth along with the bed slope compared with the friction slope different types of profiles are formed such as M (mild), C (critical), S (steep) profiles. M (mild)-If the slope is so small that the normal depth (Uniform flow depth) is greater than critical depth for the given discharge, then the slope of the channel is mild. C (critical)-if the slope if the slope’s normal depth equals its critical depth, then we call it a Critical slope, denoted by C. S (steep)-if the channel slope is so steep that a normal depth less than critical is produced, then the channel is Steep, and water surface profile designated as S (see [1] [3] ).

2.6. Rapidly Varied Flow

This flow has very pronounced curvature of the streamlines. It has pressure distribution that cannot be assumed to be hydrostatic. The rapid variation in flow regime often takes place in short span. When rapidly varied flow occurs in a sudden-transition structure, the physical characteristics of the flow are basically fixed by the boundary geometry of the structure as well as by the state of the flow. Channel expansion and channel contraction, Sharp crested weirs and Broad crested weirs can be seen as examples. Specific force before and after the flow regime can be considered same i.e., $F_1=F_2$ [2].

2.7. Unsteady Flows

When the flow conditions vary with respect to time, we call it unsteady flows. Some terminologies such as wave is defined as a temporal or spatial variation of flow depth and rate of discharge. Wave length is the distance between two adjacent wave crests or trough. Amplitude is the height between the maximum water level and the still water level [4]. Wave celerity ( $\mathcal{C}$ ) is the relative velocity of a wave with respect to fluid in which it is flowing with V. Absolute wave velocity ( $V_w$ ) is the velocity with respect to fixed reference as below:

$V_w=V\pm \mathcal{C}$ (24)

Plus sign indicates the wave is traveling in the flow direction and vice versa. For shallow water waves $\mathcal{C}=\sqrt{g{y}_0}$, where y₀ = undisturbed flow depth. Unsteady flows occur due to following reasons: surges in power canals or tunnels, surges in upstream or downstream channels produced by starting or stopping of pumps and opening and closing of control gates. Furthermore, navigation looks can also generate waves in the navigation channel. Flood waves in streams, rivers, and drainage channels due to rainstorms and snow-melt, tides in estuaries, bays and inlets. Unsteady flow commonly encountered in an open channels and deals with translatory waves. A translatory wave is a gravity wave that propagates in an open channel and results in appreciable displacement of the water particles in a direction parallel to the flow. For purpose of analytical discussion, unsteady flow is classified into two types, namely, gradually varied and rapidly varied unsteady flow. In gradually varied flow the curvature of the wave profile is mild, and the change in depth is gradual. In the rapidly varied flow the curvature of the wave profile is very large and so the surface of the profile may become virtually discontinuous. Continuity equation for unsteady flow [5] in an open channel becomes:

$D_f\frac{\partial V}{\partial x}+V\frac{\partial y}{\partial x}+\frac{\partial y}{\partial t}=0$ (25)

For a rectangular channel of infinite width, the Equation (25) may be written as:

$\frac{\partial q}{\partial x}+\frac{\partial y}{\partial t}=0$ (26)

When the channel is to feed laterally with a supplementary discharge of $q^{\prime }$ per unit length, for instance, into an area that is being flooded over a dike, then the equation becomes:

$\frac{\partial q}{\partial x}+\frac{\partial y}{\partial t}+q^{\prime }=0$ (27)

The general dynamic equation for gradually varied unsteady flow is given by (see [5] [6] ):

$\frac{\partial y}{\partial x}+\frac{\alpha V}{g}\frac{\partial V}{\partial x}+\frac{1}{g}\frac{\partial V^{\prime }}{\partial t}=0$ (28)

3. Hydraulics of Pipe Flows

Pipe flows are mainly due to pressure difference between two sections. The total head here also consists of the pressure head, elevation head and velocity head. The principle of continuity, energy, momentum is also used in this type of flow. For example, when designing a pipe, we use continuity and energy equations to obtain the appropriate pipe diameter. Then, applying the momentum equation [5] [6], for a given discharge, we get the forces that act on bends. The key factors in the design and operation of a pipeline are head losses, pressures and stresses acting on the material of the pipe, and discharge. Head loss for a given discharge relates to flow efficiency; i.e., an optimum size of pipe will yield the least overall cost of installation and operation for the desired discharge. Choosing a small pipe results in low initial costs, but due to high energy costs from significant head losses, subsequent costs may be excessively high. The design of conduit should be such that it needs least cost for a given discharge. The hydraulic aspect of the problem requires applying the one dimensional steady flow form of the energy equation. Energy equation can be written as:

$\frac{P_1}{\gamma }+\frac{{\alpha }_1{V}_1^2}{2g}+Z_1+h_p=\frac{P_2}{\gamma }+\frac{{\alpha }_2{V}_2^2}{2g}+Z_2+h_t+h_L$ (29)

where, $\frac{P}{\gamma }$ = pressure head, $\frac{\alpha V^2}{2g}$ = velocity head, Z = elevation head, h_p =

head supplied by a pump, h_t = head supplied by a turbine and h_L = head loss between 1 and 2. The kinetic energy correction factor is given by $\alpha$, and it is defines as, where u = velocity at any point in the section.

$\alpha =\frac{{\displaystyle \int }\left(u^3\text{d}A\right)}{V^{3}A}$ (30)

$\alpha$ has minimum value of unity when the velocity is uniform across the section. $\alpha$ has values greater than unity depending on the degree of velocity variation across a section. For laminar flow in a pipe, velocity distribution is parabolic across the section of the pipe, and $\alpha$ has value of 2.0. However, if the flow is turbulent, as is the usual case for water flow through the large conduits, the velocity is fairly uniform over most of the conduit section, and $\alpha$ has value near unity (typically: $1.04<\alpha <1.06$ ). Therefore, in hydraulic engineering for ease of application in pipe flow, the value of $\alpha$ is usually assumed to be unity, and the

velocity head is then simply $\frac{V^2}{2g}$.

The head supplied by a pump is directly related to the power supplied to the flow as $P_p=Q\gamma h_p$, similarly, if head is supplied to turbine, the power supplied to the turbine will be $P_t=Q\gamma h_t$. The head loss term $h_L$ accounts for the conversion of mechanical energy to internal energy (heat), when this conversion occurs, the internal energy is not readily converted back to useful mechanical energy, therefore it is called head loss. Head loss results from viscous resistance to flow (friction) at the conduit wall or from the viscous dissipation of turbulence usually occurring with separated flow, such as in bends, fittings or outlet works. Head loss is due to friction between the fluid and the pipe wall and turbulence within the fluid. The rate of head loss depends on roughness element size apart from velocity and pipe diameter. Further the head loss also depends on whether the pipe is hydraulically smooth, rough or somewhere in between. In water distribution system, head loss is also due to bends, valves and changes in pipe diameter. Head loss for steady flow through a straight pipe can be written as below:

${\tau }_0{A}_w=\Delta P{A}_r$ (31)

$\Delta P=4\mathcal{L}{\tau }_0/D$ (32)

${\tau }_0=f\rho V^2/8$ (33)

$h=\frac{\Delta P}{\gamma }=f\frac{\mathcal{L}}{D}\frac{V^2}{2g}$ (34)

This is known as Darcy-Weisbach equation. $h/L=S$, is slope of the hydraulic and energy grade lines for a pipe of constant diameter. Head loss in laminar flow can be calculated using Hagen-Poiseuille equation, which gives:

$S=\frac{32V\mu }{D^2\rho g}$ (35)

Combining above equations with Darcy-Weisbach equation provides:

$f=\frac{64\mu }{\rho VD}$ (36)

That can be further written in terms of Reynolds number:

$f=\frac{64}{R_e}$ (37)

This relation is valid for $R_e<1000$. In turbulent flow, the friction factor is a function of both Reynolds number and pipe roughness. As the roughness size or the velocity increases, flow is wholly rough and f depends on the relative roughness. Where graphical determination of the friction factor is acceptable, it is possible to use a Moody diagram. This diagram gives the friction factor over a wide range of Reynolds numbers for laminar flow and smooth, transition, and rough turbulent flow. The quantities shown in Moody Diagram (see details in [2] ) are dimensionless so they can be used with any system of units. Minor losses caused by valves, bends and changes in pipe diameter. This is smaller than friction losses in straight sections of pipe and for all practical purposes ignored. Minor losses are significant in valves and fittings, which creates turbulence in excess of that produced in a straight pipe. A minor loss coefficientK may be used to give head loss as a function of velocity head.

$h_m=K\frac{V^2}{2g}$ (38)

A flow coefficient C_v which gives a flow that will pass through the valve at a pressure drop of 1psi may be specified. Given the flow coefficient the head loss can be calculated as:

$h=18.5\times {10}^6\frac{D^4{V}^2}{2g{C}_v^2}$ (39)

The flow coefficient can be related to the minor loss coefficient by:

$K=18.5\times {10}^6\frac{D^2}{C_v^2}$ (40)

Major losses are due to friction between the moving fluid and the inside walls of the duct. The Darcy-Weisbach formula is generally considered to calculate major losses in pipes. This method is generally considered more accurate than the Hazen-Williams method. Additionally, the Darcy-Weisbach method is valid for any liquid or gas. Moody Friction Factor f can be calculated as below:

$f=\{\begin{array}{l}\frac{64}{R_e}{R}_e=\frac{VD}{\nu }\text{for}{R}_e\le 2100\to \left(\text{i}\text{.e}\text{.}\text{laminar}\text{flow}\right)\\ \frac{1.325}{\ln \left(\frac{e}{3.7D}+\frac{5.74}{R_e^{0.9}}\right)}\text{for}5000\le R_e\le {10}^8\to \left(\text{i}\text{.e}\text{.}\text{turbulent}\text{flow}\right)\to \text{i}\text{.e}\text{.}\to {10}^{-6}\le \frac{e}{D}\le {10}^{-2}\end{array}$ (41)

Major loss in pipes can also be calculated Using Hazen-Williams friction loss equation:

$V=k{C}_{HW}{R}^{0.63}{S}^{0.54}\text{where}S=\frac{h_f}{\mathcal{L}}\text{and}R=\frac{D}{4}f\text{or}\text{circular}\text{pipe}$ (42)

Hazen-Williams is only valid for water at ordinary temperatures (40˚ - 75˚F). The Hazen-Williams method is very popular, especially among civil engineers, since its friction coefficient (C_HW) is not a function of velocity or duct (pipe) diameter. Hazen-Williams is simpler than Darcy-Weisbach for calculations where one can solve for flow-rate, velocity, or diameter. When the flow conditions are changed from one steady state to another, the intermediate stage flow is referred to as transient flow. This occurs due to design or operating errors or equipment malfunction. This transient state pressure causes lots of damage to the network system [4] [7]. Pressure rise in a close conduit caused by an instantaneous change in flow velocity. If the flow velocity at a point does vary with time, the flow is unsteady. The terms fluid transients and hydraulic transients are used in practice [4]. Consider a pipe length of length L Water is flowing from a constant level upstream reservoir to a valve at downstream. Assume valve is instantaneously closed at time $t=t_0$ from the full open position to half open position. This reduces the flow velocity through the valve, thereby increasing the pressure at the valve. The increased pressure will produce a pressure wave that will travel back and forth in the pipeline until it is dissipated because of friction and flow conditions have become steady again. This time when the flow conditions have become steady again, let us call it $t_1$. So the flow regimes can be categorized into: 1) Steady flow for $t<t_0$ ; 2) Transient flow for $t_0<t<t_1$ ; and 3) Steady flow for $t>t_1$. Transient Transient-state pressures are sometimes reduced to the vapor pressure of a liquid that results in separating the liquid column at that section; this is referred to as liquid-column separation. If the flow conditions are repeated after a fixed time interval, the flow is called periodic flow, and the time interval at which the conditions are repeated is called period. The analysis of transient state conditions in closed conduits may be classified into two categories: lumped may be classified into two categories: lumped-system approach and distributed system approach. In the In the lumped system lumped system approach the conduit walls are assumed rigid and the liquid in the conduit is assumed incompressible, so that it behaves like a rigid mass, other way flow variables are functions of time only. In the distributed system approach the liquid is assumed slightly compressible. Therefore flow velocity varies along the length of the conduit in addition to the variation in time. The $1-D$ form of momentum equation for a control volume that is fixed in space and does not change shape may be written as:

${\displaystyle \sum \mathcal{F}}=\frac{\text{d}}{\text{d}t}{\displaystyle \int \rho AV\text{d}x}+{\left(\rho A{V}^2\right)}_{out}-{\left(\rho A{V}^2\right)}_{in}$ (43)

If the liquid is considered to be incompressible and the pipe is rigid, then at any instant the velocity along the pipe will be same, i.e., ${\left(\rho A{V}^2\right)}_{out}={\left(\rho A{V}^2\right)}_{in}$. Substituting for all the forces acting on the control volume we get:

$\mathcal{P}A+\gamma A\mathcal{L}\sin \alpha -{\tau }_0\pi D\mathcal{L}=\frac{\text{d}}{\text{d}t}\left(\rho AV\mathcal{L}\right)$ (44)

where, $\mathcal{P}=\gamma \left(h-\frac{V^2}{2g}\right)$, $\alpha$ = pipe slope, D = pipe diameter, $\mathcal{L}$ = pipe length,

$\gamma$ = specific weight of fluid and ${\tau }_0$ = shear stress at the pipe wall. Frictional force is replaced by $\gamma h_{f}A$, and $H_0=h+\mathcal{L}\sin \alpha$ and $h_f$ from Darcy-Weisbach friction equation. The resulting equation yields:

$H_0-\frac{f\mathcal{L}}{D}\frac{V^2}{2g}-\frac{V^2}{2g}=\frac{\mathcal{L}}{g}\frac{\text{d}V}{\text{d}t}$ (45)

When the flow is fully established $\frac{\text{d}V}{\text{d}t}=0$. The final velocity $V_0$ will be:

$H_0=\left(1+\frac{f\mathcal{L}}{D}\right)\frac{V_0^2}{2g}$ (46)

We use the above relationship to get the time for flow to establish:

$\text{d}t=\frac{2\mathcal{L}D}{D+f\mathcal{L}}\frac{\text{d}V}{V_0^2-V^2}$ (47)

If the flow changes are rapid, it is important to take into account fluid compressibility [4]. Changes are not instantaneous in the system, because in the piping system, pressure waves travel back and forth. The walls of the pipe must be rigid and the liquid slightly compressible [2] [4]. Assume that the flow velocity at the downstream end is changed from V to $V+\Delta V$, thereby changing the pressure from P to $P+\Delta P$. The change in pressure will produce a pressure wave that will propagate in the upstream direction. By assuming that the velocity reference system travels with the pressure wave, the unsteady flow condition can be translated into steady flow. Using momentum equation with control volume approach to solve for $\Delta p$ (see details in [2] [4] ). The system is now steady, the momentum equation now yield:

$PA-\left(P+\Delta P\right)A=\left(V+a+\Delta V\right)\left(\rho +\Delta \rho \right)\left(V+a+\Delta V\right)A-\left(V+a\right)\rho \left(V+a\right)A$ (48)

4. Hydraulics of Ground Water Flow

Velocity, discharge and head are the essential quantities used to characterize the flow of groundwater. They are followed by the Darcy’s law and some other metrics related to soil physical characteristics. The specific discharge is defined as the volume of water flowing through a unit area of soil per unit time, which has the unit of velocity. Specific discharge usually denoted by q and is defined as below:

$q=\frac{Q}{A}$ (49)

where Q is the flow rate of water through a cross sectional area A. The average velocity of fluid at a certain point of the porous medium is called Seepage velocity. The seepage velocity is represented by v, and is defined as:

$v=\frac{q}{n_e}=\frac{Q}{A{n}_e}$ (50)

where n_e is the effective porosity of the soil. The effective porosity, n_e is defined as ratio of the volume of continuous pore spaces (open to groundwater flow) of a soil sample to the total volume of the sample. For ground water flow, the total head denoted by h is defined as the sum of the pressure head, and the elevation head, z as follows:

$h=\frac{P}{\rho g}+z=h_p+z$ (51)

where P is the pressure at any point, $\rho$ is the density of water at the prevailing temperature. If streamlines are not perpendicular to the vertical piezometer, then h_p may include the contribution of the vertical component of the velocity head for most groundwater flows, the velocity head constitutes a very small part of the total head and is usually neglected. The important relationship is that the specific discharge at a certain point is proportional hydraulic gradient at that point and the discharge occurs in the direction of decreasing gradient, which is known as Darcy’s law:

$q_x=K\frac{\text{d}h}{\text{d}x}$ (52)

where q_x is the specific discharge along the x direction, $\frac{\text{d}h}{\text{d}x}$ is the gradient of

head causing the flow along the x direction, and K is known as hydraulic conductivity. The hydraulic conductivity can be further expressed as:

$K=\frac{\rho g}{\mu }$ (53)

If the hydraulic conductivity K is independent of position within a geologic formation, the formation is homogeneous. In this case, K is a constant. On the other hand, if K is dependent of position then it is called heterogeneous formation. If the hydraulic conductivity K is independent of the direction of measurement at a point in the geologic formation, the formation is isotropic at that point. On the other hand, if K varies with the direction of measurement at a point, the formation is anisotropic at that point. In an aquifer with layered heterogeneity, flow may occur parallel to the layers or across the layers or both. To simplify the flow calculations, it is possible to use estimates of equivalent K. If the layers are horizontal then the equivalent K for a flow parallel to the layers, K_x, and the same for a flow across the layers, K_z, are given by the following:

$K_x={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{K_i{D}_i}{\mathcal{D}}$ (54)

$K_z=\frac{\mathcal{D}}{{\displaystyle {\sum }_{i=1}^n\frac{D_i}{K_i}}}$ (55)

where $\mathcal{D}={\displaystyle {\sum }_{i=1}^n{D}_i}$ ; $D_i$ and $K_i$ are thickness and hydraulic conductivity of the i^th homogeneous layer. The Darcy law is for a homogeneous and isotropic medium. Assuming that the principal axes of K are aligned with the coordinate axes then the governing equation for a three dimensional anisotropic confined aquifer can be written as:

$\frac{\partial }{\partial x}\left(K_x\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_y\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial h}{\partial z}\right)=S_s\frac{\partial h}{\partial t}+\mathcal{R}$ (56)

Here, $\mathcal{R}$ is source or sink term. If the medium is homogeneous and isotropic, $K_x=K_y=K_z=K$. In addition, if the aquifer is horizontal and has a constant thickness b, then $\mathcal{T}=Kb$, and $S=b{S}_S$, in two dimensions the equation becomes,

$\frac{{\partial }^{2}h}{\partial x^2}+\frac{{\partial }^{2}h}{\partial y^2}=\frac{S}{\mathcal{T}}\frac{\partial h}{\partial t}$ (57)

where $\mathcal{T}$ is called the transmissivity of the aquifer, and S is called the storage coefficient. The above equation called the diffusion equation. For the steady state condition this equation becomes the well known Laplace equation:

$\frac{{\partial }^{2}h}{\partial x^2}+\frac{{\partial }^{2}h}{\partial y^2}=0$ (58)

The two dimensional flow equation for an unconfined aquifer is given by:

$\frac{\partial }{\partial x}\left(K_{x}h\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{y}h\frac{\partial h}{\partial y}\right)=S_y\frac{\partial h}{\partial t}+\mathcal{R}$ (59)

Here, $S_y$ is called the specific yield, which is the quantity of water that can be drained out of a saturated volume of porous medium due to unit lowering of the water table. The above equation is also known as the Boussinesq equation. It is a nonlinear, second order partial differential equation but linear in $h^2$.

5. Water Wave Dynamics

Small-amplitude or linear wave theory is the most elementary wave theory developed by Airy (1845). With a large range of wave parameters, this theory is easier to apply, offering a reasonable approximation of wave characteristics. It is necessary to assert the key assumptions of this theory as below. Homogeneous, ideal and incompressible fluid, irrotational flow, surface tension and Coriolis effect is neglected, pressure at the free surface is uniform and constant, small amplitude and bed is a horizontal, fixed, impermeable boundary. Two dimensional irrotationality can be stated by Equations (60) and (61):

$u=\frac{\partial \varphi }{\partial x}$ (60)

and:

$w=\frac{\partial \varphi }{\partial z}$ (61)

where, $\varphi$ is known as velocity potential is a scaler function [8]. In addition, incompressible flow implies that there is another function termed as the stream function $\psi$ (62-63).

$\frac{\partial \varphi }{\partial x}=\frac{\partial \psi }{\partial z}$ (62)

and:

$\frac{\partial \varphi }{\partial z}=-\frac{\partial \psi }{\partial x}$ (63)

$\psi$ is orthogonal to the potential function $\varphi$. The lines of constant values of the potential function (i.e., equipotential lines) are mutually perpendicular to the lines of constant values of the stream function. Both $\varphi$ and $\psi$ satisfy the Laplace equation which governs the flow of an ideal water (see Equations (64) and (65) waves).

$\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial z^2}=0$ (64)

and:

$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial z^2}=0$ (65)

The speed at which a wave form propagates is termed the phase velocity or wave celerity $C=L/T$. Where, the length L is the horizontal distance between corresponding points on two successive waves and the period T is the time for two successive crests to pass a given point. The most commonly used terms to describe water wave dynamics are depicted in Figure 4, where, a progressive

Figure 4. Definition of elementary terms used to define a progressive sinusoidal wave [8].

wave represented by the spatial variable x, temporal variable t and $\eta$ denotes the displacement of the water surface relative to the SWL. As a result their phase can be defined as $\theta =kx-\omega t$, where k and $\omega$ are described as wave number

( $k=\frac{2\pi }{L}$ ) and angular frequency ( $\omega =\frac{2\pi }{T}$ ) respectively. In this consequence, the wave profile can be written as following Equation (55):

$\eta =a\cos \left(\frac{2\pi x}{L}-\frac{2\pi t}{T}\right)=\frac{H}{2}\cos \left(kx-\omega t\right)$ (66)

Therefore, the amplitude of wave can be depicted as $a=\frac{H}{\mathrm{<mo>\; @\; </mo>}}$, where, H is the

wave height. Wave motion can be often defined in terms of some dimensionless

parameters such as wave steepness ( $\frac{H}{L}$ ), relative depth ( $\frac{d}{L}$ ) and relative height ( $\frac{H}{d}$ ). It is important to point out that sometimes relative depth and relative

height quantified as kd and ka since they differ only by a constant factor of 2π. Water waves usually classified into three following categories based on relative

depth: shallow water wave ( $0<\frac{d}{L}<\frac{1}{20}$ i.e., $0<k\ast d<\frac{\pi }{10}$ ), transitional water wave ( $\frac{1}{20}<\frac{d}{L}<\frac{1}{2}$ i.e., $\frac{\pi }{10}<kd<\pi$ ) and deep water wave ( $\frac{1}{2}<\frac{d}{L}<\infty$ i.e., $\pi <kd<\infty$ ). In general, wave phase speed can be expressed as:

$\mathcal{C}=\frac{gT}{2\pi }\tanh \left(kd\right)$ (67)

Therefore, the measurement of wave length can be:

$L=\mathcal{C}T=\frac{g{T}^2}{2\pi }\tanh \left(kd\right)$ (68)

In addition to different wave parameter, the local flow velocities and accelerations during the passage of a water wave must often be calculated using Equations (69)-(72).

$u=\frac{H}{2}\frac{gT}{L}\frac{\cosh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}\cos \left(kx-\omega t\right)$ (69)

$w=\frac{H}{2}\frac{gT}{L}\frac{\sinh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}\sin \left(kx-\omega t\right)$ (70)

${\alpha }_x=\frac{\partial u}{\partial t}=\frac{g\pi H}{L}\frac{\cosh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}\sin \left(kx-\omega t\right)$ (71)

${\alpha }_z=\frac{\partial w}{\partial t}=-\frac{g\pi H}{L}\frac{\sinh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}\cos \left(kx-\omega t\right)$ (72)

Water particle displacements from mean position for shallow/transitional-water (i.e., $\frac{d}{L}<\frac{1}{2}$ ) and deep-water waves (i.e., $\frac{d}{L}>\frac{1}{2}$ ) can be tracked by following equations:

$\xi =-\frac{H}{2}\frac{\cosh \left[k\left(z+d\right)\right]}{\sinh \left(kd\right)}\sin \left(kx-\omega t\right)$ (73)

$\zeta =\frac{H}{2}\frac{\sinh \left[k\left(z+d\right)\right]}{\sinh \left(kd\right)}\cos \left(kx-\omega t\right)$ (74)

Assuming, $\mathcal{A}=\frac{H}{2}\frac{\cosh \left[k\left(z+d\right)\right]}{\sinh \left(kd\right)}$ and $\mathcal{B}=\frac{H}{2}\frac{\sinh \left[k\left(z+d\right)\right]}{\sinh \left(kd\right)}$, we can rewrite Equation (73) and 74 in the form of $\frac{{\xi }^2}{{\mathcal{A}}^2}+\frac{{\zeta }^2}{{\mathcal{B}}^2}=1$. Therefore, water particle

displacements follows elliptical orbit for shallow/transitional-water while circular orbit for deep-water waves. In other words, for shallow/transitional-water

$\mathcal{A}=\frac{H}{2}\frac{L}{2\pi d}$ and $\mathcal{B}=\frac{H}{2}\frac{z+d}{d}$ ; for deep-water waves $\mathcal{A}=\mathcal{B}=\frac{H}{2}\exp \frac{2\pi z}{L}$.

Water pressure waves are classified into three main components: dynamic component due to acceleration, static component of pressure and atmospheric pressure. Mathematically expressed using Equation (75).

$P^{\prime }=\rho g\frac{H}{2}\frac{\cosh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}\cos \left(kx-\omega t\right)-\rho gz+P_a$ (75)

Therefore, relative pressure will be $P=P^{\prime }-P_a$. If pressure response factor is denoted by $K_z=\frac{\cosh \left[k\left(z+d\right)\right]}{\cosh \left(kd\right)}$ then the Equation (75) becomes

$P^{\prime }=\rho g\left(\eta K_z-z\right)$. Besides, pressure force, it is possible to divide water wave energy into traditional potential ( $E_p=\frac{1}{16}\rho g{H}^{2}L$ ) and the kinetic

( $E_k=\frac{1}{16}\rho g{H}^{2}L$ ) form of energy. Therefore, the total energy per wave length and unit width can be written as $E=\frac{1}{8}\rho g{H}^{2}L$. Furthermore, the rate of energy

transfer by water wave is known as energy flux, that can be quantified using the Equation (76).

$\bar{P}=\frac{1}{8}\rho g{H}^2\frac{1}{2}\left[1+\frac{2kd}{\sinh \left(2kd\right)}\right]\frac{L}{T}$ (76)

The Equation (76) can simply written as the form of $\bar{P}=En\mathcal{C}=E{C}_g$, where $C_g$ is commonly known as group velocity. Therefore, group velocity can be quantified using Equation (77).

$C_g=n\mathcal{C}=\frac{1}{2}\left[1+\frac{2kd}{\sinh \left(2kd\right)}\right]\frac{L}{T}$ (77)

Coastal hydraulics usually considers problems near the shoreline and the key processes that can affect a wave as it propagates from deep into shallow water includes: shoaling, refraction, reflection, diffraction, breaking and damping. The transformation of the wave form due to interaction with bathymetry is known as wave shoaling. Wave refraction is the changes in the direction of wave propagation due to differences in wave velocity along the crest, which are generally illustrated by lines drawn perpendicular to the wave crest in the direction of wave propagation known as wave rays. The angle of the wave along a ray follows Snell’s law. Wave diffraction is the bending of wave crests (changes in direction) due to along crest gradients in wave height. The point at which wave form becomes unstable and breaks is known as breaking point and the phenomenon is known as wave breaking. Wave breaking occurs when water particles at the crest travel much faster/farther than water particles in the trough. Wave breaking may be classified in four types as spilling, plunging, collapsing, and surging [8].

6. Differential Equations for Water Movement

Water movement is the key process in Water Resources Engineering. Most mathematical models describing water movement are consisting of differential equations deriving from the fundamental principles of the conservation of mass, energy and momentum. The continuity equation for one-dimensional unsteady gradually-varied flow, the Belanger equation, and unsteady flow equation are some example of differential Equations. The highest derivative is called the order of the differential equation and the power of its derivative of the highest order is called the degree of the differential equation. In general, the differential equations for water movement are second-order partial differential equations. The conservation laws in physics can be described by a PDE of the general form (Equation (78)):

$a\frac{{\partial }^{2}C}{\partial x^2}+b\frac{{\partial }^{2}C}{\partial x\partial y}+c\frac{{\partial }^{2}C}{\partial y^2}+d\frac{\partial C}{\partial x}+e\frac{\partial C}{\partial y}+fC+g=0$ (78)

where a, b, c, d, e, f and g may be constants or functions of x and y or of C. If a, b, c, d, e, f and g are constants or functions of x and y, the equation is linear. If $\left(b^2-4ac\right)>0$, the above equation is hyperbolic. For example wave equation is a hyperbolic equation. If $\left(b^2-4ac\right)=0$, the above equation is parabolic. For example diffusion equation is a parabolic equation. If $\left(b^2-4ac\right)<0$, the above equation is elliptic. The Laplace equation is an elliptic equation. Besides general form of PDE the conservation law can be seen as advection-diffusion process which is discussed in the following section.

7. Modeling Advection-Diffusion Processes

Advection-Diffusion process can be expressed by an equation known as advection-diffusion equation which comes from the conservation of contaminant mass. Advection refers to the transport of contaminants by the moving water where contaminants travel at the same velocity as the mass of water in which they are dissolved. Diffusion refers to the spreading of the contaminants in all direction due to the turbulence and non-uniform velocity distribution. Here, molecular diffusion is much smaller than that of turbulent diffusion. In addition to being advected and diffused, the amount of a contaminant dissolved in water may increase or decrease in time due to chemical reactions with other agents present in the environment. When such reactions are absent, the contaminant process is called a conservative process meaning that the total mass of the dissolved substance remains constant through the transport process. When such reaction processes are present, the contaminant process is called non-conservative indicating that the mass of the contaminant is growing or decaying at a certain rate.

$\frac{\partial C}{\partial t}=-u\frac{\partial C}{\partial x}+D_f\frac{{\partial }^{2}C}{\partial x^2}-rC$ (79)

Wave equation, diffusion equation and Laplace equation is also example of Advection-Diffusion process. Analytical solutions for advection-diffusion-reaction processes can be found only for simplified problems. For complicated boundary conditions or initial conditions, one needs numerical methods to get solution for the problem.

8. Transport of Particles in Water Movement

In addition to water movement, transports of sediment, salt and pollutants lead to differential equations. Advection-diffusion process is a common process in transports of sediment, salt and pollutants. Advection refers to the transport of contaminants by the moving water where contaminants travel at the same velocity as the mass of water in which they are dissolved. Diffusion refers to the spreading of the contaminants in all direction due to the turbulence and non-uniform velocity distribution. In an open channel flow, molecular diffusion is much lower compared to turbulent diffusion. The general differential form of conservation of contaminant mass equation can be expressed as:

$\frac{\partial C}{\partial t}+\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}+\frac{\partial q_z}{\partial z}=s$ (80)

where, C is the mass concentration per control volume. In addition, Fick’s law define the flux ( $q_x\mathrm{,}{q}_y$ and $q_z$ ) due to an advection and a diffusion process, that can be written as:

$q_x=uC-D_x\frac{\partial C}{\partial x}$ (81)

$q_y=vC-D_y\frac{\partial C}{\partial y}$ (82)

$q_z=wC-D_z\frac{\partial C}{\partial z}$ (83)

Apart from the general type of contaminants transporting particle, some other dissolved gas transport along with their reaction rate can also be part of hydraulics, which is very important for the supply and treatment of water.

9. Hydraulics Used in Water Supply and Treatment

The concentration of a gas (such as O₂) that can be present in solution is governed by: 1) the solubility of the gas as defined by Henry’s Law; 2) temperature; 3) presence of impurities (such as salinity); and 4) the partial pressure of the gas in the atmosphere. The saturation concentration of a gas dissolved in a liquid is

a function of the partial pressure of the gas. Henry’s law states that $P_g=H\frac{x_g}{P_T}$,

where, P_g and x_g is the mole fraction of the gas in air and in the liquid (water) respectively, P_T is the total pressure and H is the Henry’s constant that depends on the type of gas. For example, dissolved oxygen (DO) is measured in standard solution units such as mmol/L, mg/L, mL/L, or ppt. O₂ saturation is calculated as the percent of DO relative to a theoretical maximum concentration given the temperature, pressure, and salinity of the water. Well-aerated water will usually be 100% saturated. In general, lower temperatures, lower salinity, and higher atmospheric pressures lead to higher values of dissolved O₂. Higher barometric pressures lead to higher values of saturated DO. The correction for barometric

pressure may be written as $D{O}^{\prime }=D{O}_0\frac{P-u}{760-u}$, where, P is the barometric

pressure and u is the vapor pressure in mm (Hg). Microorganisms in water digest organic material as food in the presence of O₂ that lead to introduce a new term commonly known as the biochemical oxygen demand (BOD). The BOD on a water body depends on the microbial population dynamics which includes a lag phase, a constant growth phase, a stationary phase, and a decay phase. According to the popular mathematical model for the growth of BOD (i.e., the exponential model where the BOD grows asymptotically to the so-called ultimate BOD of the water), BOD can be written as $BO{D}_t=BO{D}_u\left(1-{\text{e}}^{-kt}\right)$, where k is the reaction rate constant (day⁻¹). The rate constant at T temperatures is given by $k_T=k_{20}{\left(1.047\right)}^{T-20}$. When a stream of untreated or partially treated sewage is discharged into a large body of water, the wastewater is diluted and their ultimate BOD (i.e., $L_o$ ), DO and temperature ( $T_o$ ) of the river-wastewater immediately after mixing can be given by $\left(Q_w{L}_w+Q_r{L}_r\right)/Q$, $\left(Q_{w}D{O}_w+Q_{r}D{O}_r\right)/Q$ and $\left(Q_w{T}_w+Q_r{T}_r\right)/Q$ respectively, where, $Q=Q_w+Q_r$, $Q_w$ = waste water flow rate and $Q_r$ = river flow rate. Initial oxygen deficit after mixing is given by $D_O=D{O}_{sat}-D{O}_0$. The Streeter-Phelps equations (see 84) describe the DO “sag curve” obtained as a result of simultaneous deoxygenation and reoxygenation that occurs when wastewater is discharged into receiving stream. According to the Streeter-Phelps model 5, the oxygen deficit after time t ( $t=0$ at the instant of mixing) is given by:

$D_t=\frac{k_d{L}_0}{k_r-k_d}\left({\text{e}}^{-k_{d}t}-{\text{e}}^{-k_{r}t}\right)+D_0{\text{e}}^{-k_{r}t}$ (84)

where, k_d and k_r are deoxygenation (day⁻¹) and reoxygenation rate (day⁻¹) respectively. The time at which the dissolved oxygen is critical is given by:

$t_c=\frac{1}{k_r-k_d}\ln \left(\frac{k_d{L}_o-k_r{D}_o+k_d{D}_o}{k_d{L}_o}\frac{k_r}{k_d}\right)$ (85)

At this time, the critical (maximum) O₂ deficit is given by:

$D_c=\frac{k_d{L}_o}{k_r}{\text{e}}^{-k_d{t}_c}$ (86)

The distance from the point of mixing to the critical location is given by $x_c=v{t}_c$, where, v = average stream velocity.

In the treatment of wastewater, operations in which some transformation is effected through the use of chemical reactions are termed chemical unit operations. The rate at which a chemical reaction progresses depends on the nature of the reaction and hydraulics used to describe that nature is called reactor hydraulics. Chemical reactions usually classified as zero order, first order, second order,

and saturation reactions. Zero-order reactions may be expressed as $\frac{\text{d}C}{\text{d}t}=-k$ that leading to $C_t=C_0-kt$. Similarly, for a first-order reaction follows $\frac{\text{d}C}{\text{d}t}=-kC$ that leading to $C_t=C_0{\text{e}}^{-kt}$ and for a second-order reaction follows $\frac{\text{d}C}{\text{d}t}=-k{C}^2$ that leading to $\frac{1}{C_t}=\frac{1}{C_0}+kt$. A reaction is classified as a saturation reaction if the rate of reaction saturates as the reaction progresses. This may be expressed as $\frac{\text{d}C}{\text{d}t}=\frac{kC}{a+C}$ that leading to $a\ln \frac{C_0}{C_t}+C_0-C_t=kt$. Required

hydraulic detention time for various orders of reaction also depicted in Figure 5, which is very useful for water treatment reactor. There are several further computational hydraulics applications that address water quality concerns [9].

10. Hydraulics of Sediment Transport

Sediment transport function is very crucial to predict complex river [10] hydrodynamics and morphodynamics. For decades, various methods have been used to establish sediment transport functions or formulas. They are often differ

Figure 5. (a) Dissolved oxygen sag curve (Streeter-Phelps); (b) Hydraulic detention time for various orders of reaction.

drastically from each other and from observations in the field. These formulas have been used to solve engineering and environmental problems. The sediment transport mechanics for cohesive and non-cohesive materials are different. Most of the sediment transport functions are proposed under the assumption of non-cohesive sediment. In this section, we will review of the basic concepts and approaches used in the derivation of incipient motion criteria and sediment transport functions. Due to the stochastic nature of sediment movement, incipient motion is critical for the study of sediment transport, channel degradation and stable channel design. In other words, it is difficult to define precisely at what flow condition a sediment particle will begin to move, hence it depends more or less on definition of incipient motion. Significant progress has been made on the study of incipient motion, both theoretically and experimentally. Out of all ranges, the emphasis here is more on Shear Stress Approach and Velocity Approach.

A sediment particle is at a state of incipient motion when one of the following conditions is satisfied: $F_L=W_S$, $F_D=F_R$ and $M_O=M_R$, where, overturning moment due to F_D and F_L, and resisting moment due to F_L and W_S (see Figure 6). One of the most prominent and commonly used incipient motion criteria is the Shields diagram (1936) based on shear stress approach, where, where the central premise is the shear stress $\tau$, the difference in density between sediment ( ${\rho }_S$ ) and water ( ${\rho }_w$ ), the diameter of the particle (d_s), the kinematic viscosity and the gravitational acceleration can be grouped into two dimensionless quantities that provides seminal plot to determine incipient motion (see Figure 6).

$d_s\frac{{\left(\frac{{\tau }_c}{{\rho }_w}\right)}^{1/2}}{\nu }=\frac{d{U}_{*}}{\nu }$ (87)

and:

$\frac{{\tau }_c}{d_s\left({\rho }_S-{\rho }_w\right)g}=\frac{{\tau }_c}{d_s\gamma \left[\left(\frac{{\rho }_S}{{\rho }_w}\right)-1\right]}$ (88)

Figure 6. Sediment transport modes and forces acting on a single spherical sediment particle along with Shields and HJulstrom diagram [11].

where, ${\tau }_C$ = critical shear stress at initial motion and $U_{\mathrm{<mo>\; *\; </mo>}}$ = shear velocity. Rouse (1939), White (1940), Brooks (1955), Liu (1958), Yang (1973), Govers (1987), Yalin and Karahan (1979) provides useful changes on the Shields diagram for practical use. While Lane (1953) developed stable channel design curves for trapezoidal channels with different typical side slopes, however, recently USBR (1987) synthesized stable channel design criteria based on the critical shear stress required to move sediment particles in channels under different flow and sediment conditions [11].

On the other hand, Fortier and Scobey (1926) carried out an extensive field survey of the maximum permissible mean velocity values in the canals of different materials. Hjulstrom (1935) carried out a detailed study of the movement of uniform materials at the bottom of the channels, which provides a convenient diagram (see Figure 6) based on the velocity approach. The relationship between sediment size and average flow velocity for erosion, transport and sedimentation can be illustrated by the Hjulstrom curve. Yang (1973) and Vanoni (1977) further improved the velocity approach that Govers (1987) and Talapatra and Ghosh (1983) tested experimentally [11]. Yang used theoretical fluid mechanics to calculate critical velocity by providing F_L, F_D, F_R, and W_S equations, which involve the fall velocity of the sediment particles. Sediment particle fall velocity is one of the crucial parameters used in most sediment transport functions or formulas. Different methods have been developed for the computation of sediment particle fall velocity. Rubey’s (1933) introduce fall velocity formula

${\omega }_S=F\sqrt{d_{s}g\left(G-1\right)}$, where, F is a function of particles diameter, specific gravity of sediment, kinematic viscosity of water and g, however for particle greater than 1mm, F can be considered as 0.79 [11]. Yang and Simoes (2002) use a shape

factor $SF=\frac{c}{\sqrt{ab}}$ for natural sand, where, a, b, c the length of the longest, the

intermediate, and the shortest mutually perpendicular axes of the particle, respectively. In water treatment, fall velocity commonly calculated based on the type of particle settling, that can be classified into four types: discrete particle settling, flocculent settling, hindered or zone settling, and compression settling. In addition, discrete particle fall velocity depends on the laminar, transition or turbulent region of flow that is defined in earlier section. Regime, regression, probabilistic, and deterministic approaches are the basic approaches used in the derivation of sediment transport functions or formulas. The concept of “regime” is similar to the concepts of “dynamic equilibrium” and “hydraulic geometry.” Various collections of regime equations have been formulated by various investigators, such as Blench (1969), Kennedy (1895) and Lacy (1929). Lacy’s (1929) regime equation describing the relationships among channel slope S, water discharge Q, and silt factor $f_S$ for sediment transport can be written as

$S=0.0005423\frac{f_S^{5/3}}{Q^{1/6}}$. Leopold and Maddock’s (1953) hydraulic geometry relationships commonly known as $W=a{Q}^b$, $y=c{Q}^f$ and $V=k{Q}^m$, where, W =

channel width, y = channel depth, V = average flow velocity, Q = water discharge and a, b, c, j, k, m are local constants. Yang et. al. (1981) applied the unit stream power theory for sediment transport, and the hydraulic geometry relationships between Q and S as $S=i{Q}^j$, where, i, j, constants. Shen and Hung (1972), Karim and Kennedy (1990) proposed sediment transport function adopting regression approach by considering flow velocity, sediment discharge, bed-form geometry, and friction factor. The regression approach may provide fairly accurate results due to the fact that the sediment transport is such a complex phenomenon and no single hydraulic parameter or combination of parameters can be found to describe sediment transport rate under all conditions. Einstein (1950) was a pioneer in sediment transport studies focused on a probabilistic approach, which is under the assumption of the beginning and ceasing of sediment motion can be expressed in terms of probability and the movement of bedload is a series of steps followed by rest periods. In addition to that, Einstein used the hiding correction factor and the lifting correction factor to get the best fitted theoretical result with the observed experimental data. Although Einstein’s bedload function is not common in engineering applications due to complex computational procedures, it has however been used as a theoretical foundation for the formulation of other transport functions. For example, Colby and Hembree (1955) used modified Einstein methods for the computation of total bed-material load. On the other hand, the deterministic approach is the existence of one-to-one relationship between independent and dependent variables. Commonly used independent variables are water discharge, average flow velocity, shear stress, and energy or water surface slope. More recently, the use of stream power and unit stream power have gained increasing attention as important parameters to calculate sediment concentration. Bagnold (1966) introduced the stream power concept based on classical physics, while, Engelund and Hansen (1972), and Ackers and White (1973) later used the concept as the theoretical basis for developing their sediment transport functions. Yang (1972) defines unit stream power as the rate of energy per unit weight of water available for transporting water and sediment in an open channel with reach length x and

total drop of y can be written as $\frac{\text{d}y}{\text{d}t}=\frac{\text{d}x}{\text{d}t}\frac{\text{d}y}{\text{d}x}=VS$, where, V = average flow

velocity, and S = energy or water surface slope. Besides that, Velikanov (1954) derived his transport function from the gravitational power theory and Pacheco-Ceballos (1989) derived a sediment transport function based on power balance between total power available and total power expenditure in a stream. This principle is also applicable to understanding the evolution of the natural channel network [12] [13] [14] [15] [16].

In this section, we comprehensively reviews basic approaches and theories used in the determination of noncohesive sediment concentration. In the following section we will discuss some basic finite differences numerical techniques to solve governing equations. Finite element and finite volume are another two other numerical technique which will be covered under subsequent sections.

11. Finite Difference Method

In Finite difference methods (FDM), partial derivatives are replaced with finite differences approximations, which results in a set of algebraic equations. These algebraic equations can be solved explicitly or implicitly. FDM uses a structured grid and thus requires grid transformation in multidimensional flow case. A typical finite difference grid can be shown in Figure 5. Finite difference method is popular because of its simplicity in formulation and implementation. Based on different ways of discretization, different numerical schemes are developed. Some of them are shown as bellows:

In a forward time difference, any partial time derivative of a function C can be expressed as Equation (89):

$\frac{\partial C}{\partial t}=\frac{C_j^{n+1}-C_j^n}{\Delta t}$ (89)

In a forward space difference, any partial space derivative of a function C can be expressed as Equation (90):

$\frac{\partial C}{\partial x}=\frac{C_{j+1}^n-C_j^n}{\Delta x}$ (90)

In a backward space difference, any partial space derivative of a function C can be expressed as Equation (91):

$\frac{\partial C}{\partial x}=\frac{C_j^n-C_{j-1}^n}{\Delta x}$ (91)

In a centered space difference, any partial space derivative of a function C can be expressed as Equation (92):

$\frac{\partial C}{\partial x}=\frac{C_{j+1}^n-C_{j-1}^n}{2\Delta x}$ (92)

where $\Delta t$ and $\Delta x$ is the time and space discretization respectively. When space discretizations are done on the current time level, then the scheme is known as an explicit scheme. For example Equations (90) to (92) are the examples of explicit space discretized scheme. On the other hand, when space discretizations are done on the next time level or using both the current and next time levels, the generated scheme is known as an implicit scheme. For example, if Equation (91) is discretized as follows:

$\frac{\partial C}{\partial x}=\frac{C_j^{n+1}-C_{j-1}^{n+1}}{\Delta x}$ (93)

or:

$\frac{\partial C}{\partial x}=\frac{1}{2}\left(\frac{C_j^n-C_{j-1}^n}{\Delta x}+\frac{C_j^{n+1}-C_{j-1}^{n+1}}{\Delta x}\right)$ (94)

then both will generate implicit schemes. An explicit scheme will not generate any systems of equation. Therefore, an explicit scheme can be solved directly. That means, unknowns at any nodal value can be calculated without solving other nodal values. On the contrary, an implicit scheme will generate a system of equations, and therefore, unknowns at any nodal value can not be solved independently. Apart from that another popular implicit scheme is Box Finite Difference Scheme, which is also known as four-point implicit scheme. According to this scheme, time derivative can be expressed as:

$\frac{\partial C}{\partial t}=\frac{1}{2}\left(\frac{C_{j+1}^{n+1}-C_{j+1}^n}{\Delta t}+\frac{C_j^{n+1}-C_j^n}{\Delta t}\right)$ (95)

and space derivative can be expressed as:

$\frac{\partial C}{\partial x}=\frac{1}{2}\left(\frac{C_{j+1}^{n+1}-C_j^{n+1}}{\Delta x}+\frac{C_{j+1}^n-C_j^n}{\Delta x}\right)$ (96)

and C can be expresses as:

$C=\frac{1}{4}\left(C_j^n+C_{j+1}^n+C_j^{n+1}+C_{j+1}^{n+1}\right)$ (97)

In order to improve accuracy, further schemes such as the Crank-Nickelson Scheme, Box Finite Difference Scheme, Lax-Wendroff Scheme, MacCormack Scheme were proposed to solve advection-diffusion processes. Prior to solving the advection-diffusion equation or any governing equation, the initial condition of the system has to be specified. Solution of the governing equations also require boundary conditions. The most common type of boundary conditions are Dirichlet, Neumann, and Cauchy boundary conditions. Dirichlet condition occurs when a portion of the boundary is at a prescribed concentration level (i.e., C

specified), Neumann condition occurs when a portion of the boundary has specified flow which is crossing the boundary curve normally (i.e., $\frac{\partial C}{\partial n}$ specified),

and Cauchy condition is a mixed boundary condition which represents some combination of the Direchlet and Neumann conditions on the boundary curve. It is important to mention that numerical method used to reformulate governing equations in an approximate manner where continuous operations replaced by discrete operations that causes truncation error. In addition, size of the discretized solution domain (i.e., cells) causes round off error (Figure 7). Model discretization and numerical solution should be consistent (i.e., correct), convergent (i.e., $\Delta x\to 0\to \epsilon \to 0$ ), and stable (i.e., error should not increase).

12. Finite Element Method

Another numerical technique used to obtain an approximate solution to a governing differential equation is the Finite Element Method. This is relatively new to the field of computational hydraulics and is supported by rigorous mathematical theory. FEM based on variational form of PDE, which derived using the method of weighted residuals. In this method, under a given boundary condition, the PDE (i.e., strong form) is used and the residuals are determined for the assumed solution due to the fact that the assumed solution is not exact. Then push the residual to zero at various points or intervals using different weighted residual method. There are three commonly used weighted residual methods (i.e., Collocation, subdomain, and Galerkins method) in FEM. Among them Galerkins method is popular, where the integral of residual equation R, multiplied by a weighting function, w, is pushed to be zero. The weighting functions are chosen to be the same form as each part of the approximate solution. If the residual, $R\left(u\right)$, for a given PDE, should equal 0 for the true solution. In other words, for approximate solution, non-zero residual measures accuracy $u_h\approx u$. Here, we require that the weighted residual equal zero:

Figure 7. (a) Schematic 1D space and time discretization; (b) Computational error as a function of $\Delta x$.

${\displaystyle {\int }_{\Omega }}wR\left(u\right)\text{d}x=\mathrm{0\; <mo>\; ;\; </mo>}\forall w\in W$ (98)

and $u_h$ is a linear combination of basic functions, can be expressed as

$u_h\left(x,t\right)={\displaystyle \underset{j}{\sum }}{u}_j\left(t\right){\varphi }_j\left(x\right)$ (99)

Eventually, transforms governing equation into a system of linear equations ( $Ax=b$ ), which can be solved using matrix operations. The FEM can be explain by three phases prepossessing, processing/solution and post-processing. In prepossessing phases, model geometry is discretized or subdivided into a finite number of smaller pieces, each piece is called an element, and all intersection points are called nodes. The process of dividing geometry is called meshing, which is one of the key steps of the FEM. After that, one element is considered at a time and assumes a shape function that depicts the physical behavior of the element. The stiffness matrix of each element is then developed and assembled them into a global matrix system for the entire geometry. Apply the required boundary condition and initial condition to the global matrix system. In solution phase, primary unknown is computed from the global system of linear equation. Finally, in the post-processing phase, other derived variables are determined on the basis of the nodal value of the primary unknowns.

For example, SLIM is a numerical model. The full form of SLIM is Second-generation Louvain-la-Neuve Ice-ocean Model. It is a hydrodynamic model based on finite element method (FEM). The special advantage of finite element method is that, user can use this technique for unstructured grids. One can refined computational grid arbitrarily in the areas of interest that focusing the computational power where it is needed. In this model no need of nested grids. A single model is able to resolve both the large-scale features, such as in the open sea and also small-scale phenomena in rivers, coasts and estuaries. SLIM consists of a 1D river model, a 2D depth averaged model and 3D model. For 1D river model consists of linear river where variable river width and cross-section, 2D model the domain is divided into triangular elements allowing accurate representation of complex geometry and 3D model uses triangular prismatic elements that are formed by extending the 2D mesh vertically.

13. Finite Volume Method

Besides computational hydraulics, majority of computational fluid dynamics (CFD) codes based on Finite Volume Method (FVM). This method is based on integral form of PDE, where the governing equation integrated over control volume. For example, the integral form of 1D advection-diffusion conservation law can be written as:

${\displaystyle {\int }_v}\frac{\partial C}{\partial t}\text{d}v=-{\displaystyle {\int }_v}u\frac{\partial C}{\partial x}\text{d}v+{\displaystyle {\int }_v}{D}_f\frac{\partial }{\partial x}\frac{\partial C}{\partial x}\text{d}v$ (100)

The integral form is based on a control volume v. This v volume can be cell-centered or vertex-centered as shown in Figure 8. Typically, the special model domain divides into discrete control volumes. Then an assumption is made on

Figure 8. (a) Cell-centered and (b) Vertex-centered finite volume mesh along with (c) Midpoint approximation rule.

how the value of integral changes within the control volume v, usually average value is taken at the center of the control volume v.

$C_i\left(t\right)=\frac{1}{\left|V_i\right|}{\displaystyle {\int }_{V_i}}C\left(x,t\right)\text{d}x$ (101)

The suitable time integration scheme is selected to solve the integral of Equation (101). Midpoint rule, Trapezoidal rule and Simpson’s rule are the most common approximations of integrals. In this method, all fluxes usually converted to surface integrals. This is an exact procedure due to the fact that the volumes are defined as polygons or polyhedra.

14. Concluding Remarks

Hydraulics is a prerequisite for those interested in a career in Water Resources Engineering. The fundamental principles of hydraulics, the computational aspects of hydraulic system analysis and design, and various engineering applications of the concepts are discussed in this manuscript. This manuscript can assist in developing an appropriate learning environment for comprehending and conducting pipe and open channel flow research.

Acknowledgements

The author wishes to express his gratitude to the editor and anonymous reviewers for their constructive suggestions and comments, which substantially improved the manuscript’s presentation and interpretation.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Çengel, Y.A., Turner, R.H., Cimbala, J.M. and Kanoglu, M. (2008) Fundamentals of Thermal-Fluid Sciences. McGraw-Hill, New York, 703.
[2]	Sarker, S. (2021) Hydraulics Lab Manual. 1-66. https://doi.org/10.31224/osf.io/mxcvw
[3]	Kumar, M.S.M. (2009) Computational Hydraulics. Lecture of Civil Engineering. Indian Institute of Science, Bengaluru, 1-404. https://nptel.ac.in/courses/105/108/105108125
[4]	Sarker, S. (2021) A Story on the Wave Spectral Properties of Water Hammer. 1-12. https://doi.org/10.31224/osf.io/nhuzq
[5]	Sarker, S. (2021) Separation of Flood Plain Flow: 1-D Momentum Equation Solver. 1-7. https://doi.org/10.31224/osf.io/sjcmv
[6]	Reza, A.A., Sarker, S. and Asha, S.A. (2014) An Application of 1-D Momentum Equation to Calculate Discharge in Tidal River: A Case Study on Kaliganga River. Technical Journal River Research Institute, 2, 77-86.
[7]	Sarker, S. (2021) Pipe Network Design and Analysis: An Example with WaterCAD. 1-10. https://doi.org/10.31224/osf.io/c3aky
[8]	USACE (2002) Coastal Engineering Manual (CEM). US Army Corps of Engineers.
[9]	Duan, Y.Z. and Liu, G.J. (2016) Water Resource Pricing Study Based on Water Quality Fuzzy Evaluation: A Case Study of Hefei City. Computational Water, Energy, and Environmental Engineering, 5, 99-111. https://doi.org/10.4236/cweee.2016.54010
[10]	Khan, I., Ahammad, M. and Sarker, S. (2014) A Study on River Bank Erosion of Jamuna River Using GIS and Remote Sensing Technology. International Journal of Engineering Development and Research, 2, 3365-3371. http://www.ijedr.org/papers/IJEDR1404002.pdf
[11]	Yang, C.T. (2003) Chapter 3: Noncohesive Sediment Transport. In: Sediment Transport: Theory and Practice, Krieger Publishing, Malabar, 1-111.
[12]	Sarker, S. (2021) Investigating Topologic and Geometric Properties of Synthetic and Natural River Networks under Changing Climate. UCF STARS.
[13]	Sarker, S., Veremyev, A., Boginski, V. and Singh, A. (2019) Critical Nodes in River Networks. Scientific Reports, 9, Article No. 11178. https://doi.org/10.1038/s41598-019-47292-4
[14]	Sarker, S., Veremyev, A., Boginski, V. and Singh, A. (2019) Spectral Properties of River Networks. AGUFM, EP51C-2107.
[15]	Sarker, S., Veremyev, A., Boginski, V. and Singh, A. (2018) On Critical Nodes in River Networks. American Geophysical Union Fall Meeting, EP33D-2446.
[16]	Sarker, S. and Singh, A. (2017) On the Topologic Properties of River Networks. American Geophysical Union Fall Meeting, IN33B-0128.