COMPUTATIONAL HYDRAULICS

SPRING 2014

INSTRUCTORS: PROF. DR. NECATİ AĞIRALİOĞLU and DOÇ. DR. MEHMET ÖZGER

CREDIT HOURS (3,0)

COURSE OUTLINE

1. Introduction to Basic concepts

2. Newton-Raphson Method
3. Least-Square Method
4. Numerical Interpolation
5. Numerical Integration
6. Method of Finite Difference
7. Basic forms of Finite Difference
8. Explicit and Implicit Methods
9. Stability, Consistency and Convergence
10. Initial and Boundary Conditions
11. Form and Occurrence of Some Partial Differential Equations
12. Numerical Solution of Parabolic Equations
13. Numerical Solution of Hyperbolic Equations
14. Numerical Solution of Elliptic Equations
15. Method of Characteristics
16. Applications in Civil Engineering

REFERENCES

1. AĞIRALİOĞLU, N., Sonlu Farklar Metodu, Unpublished Notes, 1985.

2. CHAPRA, S.C., CANALE, R.P., Numerical Methods for Engineers, McGraw Hill 1989.
3. GOLUB, G.H., ORTEGA, J.M., Scientific Computing and Differential Equation: An Introduction to
Numerical Methods, Academic Pres, 1993.
4. JENKINS, W.M., COULTHARD, JM., JESUS, E.C.de, BASIC Computing for Civil Engineers, Van Nostrand
Reinhold, 1983.
5. KOUTITAS, C.G., Elements of Computational Hydraulics, Pentech Pres, 1983.
6. ROACHE, P.J., Computational Fluid Dynamics, Hermosa Publishers, 1979.
7. FENTON, J., Computational Hydraulics, Wien, 2010.
8. Lick, W., Sediment and Contaminant Transport in Surface Waters, Taylor and Francis Group,2009.
9. KUMAR, M.S.M, Computational Hydraulics, India, 2010.

CONDICTION OF THE COURSE:

1. 70% attendance is strictly required.
2. Grading of the course will be as follows:
First midterm exam : 15
Second midterm exam : 15
Final exam : 30
Attendance : 10
Homework : 20
Quiz : 10
Total : 100

1
 CONTENTS

1. INTRODUCTION TO COMPUTATIONAL HYDRAULICS

1.1 WHAT IS COMPUTATIONAL HYDRAULICS?
1.2 BASIC PRINCIPLES OF COMPUTATIONAL HYDRAULICS
1.3 FORMS OF THE GOVERNING WATER-FLOW EQUATIONS
1.3.1 Integral (Control-Volume) Approach
1.3.2 Differential Equations
1.3.3 Numerical Integration of Differential Equations
1.4 THE MAIN DISCRETIZATION METHODS
1.5 TERMS OF MODELING HYDRAULIC MODELING
1.5.1 Introduction
1.5.2 General Terms in Modeling
1.5.3 Scaling and Design of Hydraulic Models
1.5.4 Types of Hydraulic Models
1.5.5 Types of Mathematical Models
1.5.6 Terms Referring to Numerical Modeling

2. NUMERICAL ANALYSIS
2.1 INTRODUCTION
2.2 NEWTON-RAPHSON METHOD
2.3 LEAST SQUARES METHOD
2.4 CURVE FITTING
2.5 NUMERICAL INTEGRATION
2.5.1 Numerical Integration By Trapezoidal Rule
2.5.2 Numerical Integration By Simpson’s Rule
2.6 SOLVED PROBLEMS
2.6.1 Calculation of the water at depth downstream of a spillway
2.6.2 Example of Regression
2.6.3 Example for Polynomial Fitting
2.6.4 Example for Numerical Integration
2.6.5 Example for Numerical Integration
2.6.6 Example for Numerical Integration by FORTRAN Programming
2.6.7 Example for Numerical Integration By BASIC Programming

3. FINITE DIFFERENCE METHODS

3.1 INTRODUCTION
3.2 DIFFERENCE EQUATION
3.3 PARTIAL DIFFERENTIAL EQUATIONS
3.4 DERIVATION OF SOME BASIC FINITE DIFFERENCE FORMS
2
 3.4.1 Taylor Series Expansions
3.4.2 Polynomial Fitting
3.4.3 Integral Method
3.5 TYPES OF FINITE DIFFERENCE SCHEMES
3.5.1 Basic Forms
3.5.2 Geometric Representation for Basic Finite Difference Forms
3.5.3 Explicit and Implicit Finite Difference Schemes
3.5.4 Comparison of Explicit And Implicit Schemes
3.6 ORDER OF SCHEMES
3.6.1 First Order Schemes
3.6.2 Second Order Schemes
3.7 PROCEDURE STEPS FOR FINITE DIFFERENCE SOLUTIONS
3.8 SOLVED PROBLEMS
3.8.1 Example of Finite Difference Approximation
3.8.2 An Example Of Finite Difference Solution
3.8.3 Coefficients For Derivatives

4. STABILITY, CONSISTENCY AND CONVERGENCE

4.1 INTRODUCTION
4.2 DEFINITIONS
4.3 STABILITY
4.4 CONSISTENCY AND CONVERGENCE
4.5 VON NEUMANN ANALYSIS
4.6 HIRT’S STABILITY ANALYSIS
4.7 SOLVED PROBLEM

5. GOVERNING FLOW EQUATIONS

5.1 INTRODUCTION
5.2 PRIMITIVE EQUATIONS
5.3 STREAM FUNCTION AND VORTICITY TRANSPORT EQUATIONS FOR PLANAR FLOWS
5.4 EQUATIONS FOR PLANAR FLOWS
5.5 CONSERVATION FORM
5.6 NORMALIZING SYSTEMS
5.8 ONE DIMENSIONAL MODEL TRANSPORT EQUATIONS

6. BOUNDARY AND INITIAL CONDITIONS

6.1 INTRODUCTION
6.2 DEFINITION OF STREAM FUNCTION AND EQUIPOTENTIAL FUNCTION
6.3 BOUNDARY CONDITIONS FOR PRIMITIVE EQUATIONS
6.4 DETERMINATION OF VELOCITY COMPONENTS FROM NON DIMENSIONAL STREAM FUNCTION
6.5 BOUNDARY CONDITIONS FOR OFF- STREAM FUNCTION EQUATIONS
3
 6.6 BOUNDARY CONDITIONS FOR VELOCITY AND STREAM FUNCTION EQUATIONS
6.7 BOUNDARY CONDITIONS FOR SHARP CORNERS
6.8 BOUNDARY CONDITIONS FOR IRREGULAR BOUNDARIES
6.9 INITIAL CONDITIONS
6.10 NUMBER OF BOUNDARY CONDITIONS
6.11 SOLVED PROBLEM

7. NUMERICAL SOLUTIONS OF ELLIPTIC EQUATIONS

7.1 INTRODUCTION
7.2. SOLUTION METHODS FOR POISSON’S EQUATION
7.2.1 Governing Equation
7.2.2 Boundary Conditions
7.2.3 Solution Methods
7.3 SOLUTION METHODS FOR LAPLACE’S EQUATION
7.3.1 Direct solution
7.3.2 Relaxation
7.3.2.1 Jacobi
7. 3.2.2 Gauss-Seidel
7.3.2.3 Red-Black ordering
7.3.2.4 Successive Over Relaxation (SOR)
7.4 SOLVED PROBLEMS

8. NUMERICAL SOLUTION OF PARABOLIC EQUATIONS

8.1 INTRODUCTION
8.2 SCHEMES FOR PARABOLIC EQUATIONS
8.3 COMPARISON OF SCHEMES FOR NUMERICAL SOLUTION OF PARABOLIC EQUATIONS
8.4 SOLVED PROBLEMS

9. NUMERICAL SOLUTION OF HYPERBOLIC EQUATIONS

9.1 INTRODUCTION
9.2 SCHEMES FOR HYPERBOLIC EQUATIONS
9.3 SOLVED PROBLEMS

10. MATHEMATICAL MODELS IN OPEN CHANNEL FLOWS

10.1 INTRODUCTION
10.2 THE SAINT VENANT EQUATIONS
10.3 OPEN CHANNEL FLOW MODELS
10.3.1. Steady Flow
4
10.3.2 Quasi-Steady Approximation
10.3.3 Kinematic and Diffusion Wave Models
10.3.4 Inertial and Dynamic Wave Models
10.4 APPLICABILITY OF KINEMATIC AND DIFFUSION WAVE MODELS
10.5 UNSTEADY FLOW MODELS IN OPEN CHANNEL FLOW
10.5.1 Summary of Saint Venant Equations
10.5.2 Kinematic Wave Model
10.5.3 Numerical Solutions of the Kinematic Value
10.5.4 Nonlinear Kinematic Wave Scheme
10.6 SOLVED PROBLEMS

11. METHOD OF CHARACTERISTICS

11.1 INTRODUCTION
11.2 CHARACTERISTICS OF FIRST-ORDER PARTIAL DIFFERENTIAL EQUATIONS
11.3 SOLUTION OF COLOR EQUATION
11.4 KINEMATIC WAVE MODELING OF OVERLAND FLOW
11.4.1 Governing equations
11.4.2. Boundary and Initial Conditions
11.4.3 Analytical solution using characteristics method
11.4.4 Mathematical (analytical) Solutions
11.5 SOLVED PROBLEMS
11.4 .1 Example

REFERENCES

5
 CHAPTER 1

INTRODUCTION

1.1 WHAT IS COMPUTATIONAL HYDRAULICS?

Computational hydraulics (CH) or computational fluid dynamics (CFD) is the use of computers
and numerical methods to solve problems involving water and fluid flows.
CH has been successfully applied in many areas of water resources engineering and civil
engineering. Applications in civil engineering include: 1. Computation of normal depth, 2. Computation
of water surface profiles, 3. Contaminant transport in streams through an advection – dispersion
process, 4. Steady state groundwater flow system, 5. Flow in a pipe network, 6. Computation of
kinematic and dynamic equations, hydraulic structures such as weirs and spillways, sediment
transport.
Computational methods are narrower than numerical methods. In computational methods
the solution can be found by using only finite by using many numerical methods difference analysis.
The solution can be found by computer. In the numerical methods the solution can be found by hand,
calculator, or computer.
Engineering analysis methods can be classified as (1) Theoretical Methods, (2) Computational
Methods (Closer to Experimental Methods), (3) Experimental Methods.
Some remarks for this course:
1) Subject is almost new.
2) Its level is ideal for first year graduate students.
3) The course will be open for two semesters.
4) The second one will be related to special topics in computational methods.
5) Students would be enthusiastic to this course much more by applying computational
methods in their calculations.
This course is at an introductory level throughout as its aim is to facilitate the engineer, having
a separate knowledge and experience of engineering and computer programming.

1.2 BASIC PRINCIPLES OF COMPUTATIONAL HYDRAULICS

The approximation of a continuously-varying quantity in terms of values at a finite number of
points is called discretisation. The fundamental elements of any CH simulation are:
(1) The flow field is discretised; i.e. field variables (u, v, w, p, …) are approximated
by their values at a finite number of nodes.
(2) The equations of motion are discretised (approximated in terms of values at nodes):
control-volume or differential equations (continuous) algebraic equations
(discrete)
(3) The resulting system of algebraic equations is solved to give values at the nodes.

Continuous Discrete
curve approximation

Figure 1.1 Continuous curve and discrete approximation

6
 The main stages in a CH simulation are:
Pre-processing: (a) formulation of the problem (governing equations and boundary
conditions), (b) construction of a computational mesh (a set of control volumes).
Solving: (a) discretisation of the governing equations, (b) solution of the resulting algebraic
equations.
Post-processing: (a) visualization (graphs and plots of the solution); (b) analysis of results
(calculation of derived quantities: forces, flow rates, ... )

1.3 FORMS OF THE GOVERNING WATER-FLOW EQUATIONS

The equations governing fluid motion are based on fundamental physical principles:
· mass: change of mass = 0
· momentum: change of momentum = force × time
· energy: change of energy = work + heat
(In fluid flow these are usually expressed as rate equations; i.e. rate of change = …)
Additional equations may be applied for non-homogeneous fluids (e.g. containing dissolved
chemicals or imbedded particles).
When applied to a fluid continuum the same conservation principles may be expressed
mathematically as either: (a) integral (i.e. control-volume) or, (b) differential equations.

1.3.1 Integral (Control-Volume) Approach

This considers how the total amount of some physical quantity (mass, momentum, energy, …)
is changed within a specified region of space (control volume).
For an arbitrary control volume the balance of a physical quantity over an interval of time is
change = amount in - amount out + amount created

In hydraulics or in fluid mechanics this is usually expressed as rate by dividing with the time
interval (and transferring the net amount passing through the boundary to the LHS of the equation):

[RATE OF CHANGE (inside volume)] + [NET FLUX (out of boundary)] = [SOURCE(inside volume)]
(1.1)

The flux, or rate of transport across a surface, is composed of:

advection – movement with the fluid flow;
diffusion – net transport by random molecular or turbulent motion.

[RATE OF CHANGE (inside volume)] + [( ADVECTION + DIFFUSION)(through boundary of

volume)] = [ SOURCE (inside volume)] (1.2)

The important point is that this is a single, generic equation, irrespective of whether the
physical quantity concerned is mass, momentum, chemical content, etc. Thus, instead of
dealing with too many different equations we can consider the numerical solution of a generic
scalar-transport equation.
The finite-volume method, which is the subject of this course, is based on approximation of
these control-volume equations.

7
 V

Figure 2.2 Finite Volume Scheme

1.3.2 Differential Equations

In regions without shocks, interfaces or other discontinuities, the fluid-flow equations can also
be written in equivalent differential forms. These describe what is going on at a point rather than over
a whole control volume. Mathematically, they can be derived by making the control volumes infinitely
small.
The finite-difference method is based on the direct approximation of a differential form of the
governing equations.

1.3.3 Numerical Integration of Differential Equations

Numerical Integration of Ordinary Differential Equations: 1. Euler Method,
2. Modified Euler Method, 3. Runga - Kutta Method, 4. Predictor – Corrector Method.
Numerical Integration of Partial Differential Equations: 1. Characteristics Method, 2. Finite
Difference Method, 3. Finite Element Method, 4. Finite Volume Method, 5. Spectral Method, 6.
Boundary Element Method, etc.

1.4 SOLUTION METHODS OF PARTIAL DIFFERENTIAL EQUATIONS

1.4.1 Finite Difference Method
In computational hydraulics, finite difference, finite element, and finite volume methods are
widely used in CH and historically finite difference method is the oldest one of the three mentioned
above methods. Techniques in this method are published as early as 1910 by L. F. Richardson. Seminal
paper by Courant, Fredrichson and Lewy (1928) derived stability criteria for explicit time stepping.
First ever numerical solution: flow over a circular cylinder by Thom (1933). Scientific American article
by Harlow and Fromm (1965) clearly and publicly expresses the idea of “computer experiments” for
the first time and CFD is born. Its advantage its easiness to implement. The disadvantages are
restricted to simple grids and does not conserve momentum, energy, and mass on coarse grids.
Basic methodology in finite difference: The domain is discretized into a series of grid points. A
“structured” (i, j, k) mesh is required. The governing equations (in differential form) are discretized
(converted to algebraic form). First and second derivatives are approximated by truncated Taylor
series expansions. The resulting set of linear algebraic equations is solved either iteratively or
simultaneously.

1.4.2 Finite Element Method

Earliest use was by Courant (1943) for solving a torsion problem. Clough (1960) gave the
method its name. Method was refined greatly in the 60’s and 70’s, mostly for analysis of structural
mechanics problems. FEM analysis of fluid flow was developed in the mid- to late 70’s. Advantages:
highest accuracy on coarse grids. Excellent for diffusion dominated problems (viscous flow) and
viscous, free surface problems. Disadvantages: slow for large problems and not well suited for
turbulent flow.

1.4.3 Finite Volume Method

First well-documented use was by Evans and Harlow (1957) at Los Alamos and Gentry, Martin and
Daley (1966). It was attractive because while variables may not be continuously differentiable across
shocks and other discontinuities but mass, momentum and energy are always conserved. FVM
benefits from an advantage in memory use and speed for very large problems, higher speed flows,
turbulent flows, and source term dominated flows (like combustion). Late 70’s, early 80’s saw

8
development of body-fitted grids. By early 90’s, unstructured grid methods had appeared.
Advantages: basic FV control volume balance does not limit cell shape; mass, momentum, energy
conserved even on coarse grids; efficient, iterative solvers are well developed. Disadvantages: false
diffusion when simple numerics are used.
Finite volume: basic methodology: Integrate the differential equation over the control volume and
apply the divergence theorem. To evaluate derivative terms, values at the control volume faces are
needed and then it is required to make an assumption about how that values vary. The result is a set
of linear algebraic equations as one for each control volume solved iteratively or simultaneously.
Cells and nodes: Using finite volume method, the solution domain is subdivided into a finite number of
small control volumes (cells) by a grid. The grid defines the boundaries of the control volumes while
the computational node lies at the center of the control volume. The advantage of FVM is that the
integral conservation is satisfied exactly over the control volume.

1.5 TERMS OF HYDRAULIC MODELING

1.5.1 Introduction
ASCE (American Society of Civil Engineers) Task Committee on Glossary of Hydraulic Modeling Terms of
the Committee on Research of the Hydraulics Division studied on the hydraulic modeling terms. A
glossary is presented for the purposes of study by researchers with a view to generating a common
terminology among those involved with model in hydraulic phenomena. The glossary is divided into
five categories (Chesnutt et all. 1982).
1. General Terms in Modeling
2. Scaling and Design of Hydraulic Models
3. Types of Hydraulic Models
4. Types of Mathematical Models
5. Terms referring to numerical models
This glossary helps to generate a common terminology among those involved with modeling hydraulic
phenomena.

1.5.2 General Terms in Modeling

Air Model – Physical Models using air as fluid to model a hydraulic phenomena.
Adjustment – Variation of the parameters in a model to ensure a close reproduction by the
model of a set of prototype conditions.
Analog Model – Model using properties and behavior of the modeling materials that obey the
same mathematical laws but not the same physical laws as the prototype to produce analogous
results, e.g., an electrical analog model uses electrical potential to represent prototype velocity
potential and electrical current to represent prototype flow.
Calibration – Process of checking, adjusting, or standardizing operating characteristics of
instruments and model appurtenances on a physical model or coefficients in a mathematical model.
The process of evaluating the scale readings of an instrument in terms of the physical quantity to be
measured.
Conceptual Model – Simplification of prototype behavior used to demonstrate concepts.
Confirmation – Process in which a model of a specific design is built and tested to confirm that
a design is adequate, with respect to design conditions and that no major phenomenon has been
overlooked.
Distorted Model – Hydraulic Model in which horizontal and vertical scales are different.
Hybrid Model – Model combining at least two modeling techniques, e.g., physical and
numerical, in a closely coupled fashion.
Hydraulic Model – Physical model using water as a fluid.
Mathematical Model – Model using mathematical relationships to represent the prototype.
Model – A reproduction of the prototype, generally small-scale, but it may be larger or
geometrically distorted.

9
 One Dimensional Model – Model defined on one space coordinate, i.e., variables are averaged
over the other two directions, e.g., pipe flow models and numerical models of long wave propagation
in a narrow channel.
Physical Model – Model using the physical properties and behavior of modeling materials to
represent the prototype.
Prototype – The full-sized structure, system process, or phenomenon being modeled.
Quasi Steady-State Model – Model in which time-dependent variables are simulated by a
sequence of steady-state models.
Quasi Three-Dimensional Model – A combination of several two-dimensional models used to
simulate variations in three dimensions.
Simulation – Replication of the prototype using a model.
Steady-State Model – Model in which the variables being investigated are independent of
time.
Three-Dimensional Model – Model defined on three spaces coordinates, e.g., coastal models
and numerical models of explosions.
Two-Dimensional Model – Model defined on two space coordinates, i.e., variables are
averaged over the third direction (e.g., wave flume experiments, numerical models of storm surges).
Unsteady-State Model – Model in which the variables being investigated are time dependent.
Validation – Comparison of model results with a set of prototype data not used for verification.
Comparison includes: (1) Using a data set very similar to the verification data to determine the validity
of the model under conditions for which it was designed; (2) using a data set quite different from the
verification data to determine the validity of the model under conditions for which it was not designed
but could possibly be used; and (3) using post-construction prototype data to determine the validity of
the predictions based on model results.
Verification – Check of the behavior of an adjusted model against a set of prototype
conditions.

1.5.3 Scaling and Design of Hydraulic Models

The following symbols are used in this section: D = sediment grain size, L; g = gravitational
acceleration, LT-2; K = bulk modulus of elasticity, ML-1T-2; l = length, L; p = pressure, ML-1T-2; p0 =
reference pressure, ML-1T-2; pv = vapor pressure, ML-1T-2; V = velocity, LT-1 ; V* = shear velocity, LT-1; ν =
kinematic viscosity, L2T-1; ρ = density of water, ML-3; ρ’ = density difference between two fluids, ML-3; ρs
= density of sediment, ML-3; and σ = surface tension, MT-2.
Boundary Conditions – Conditions entered at the spatial boundaries of the model. (Slightly
different definition more appropriate to numerical models may be found in List 5.)
Boundary Effect – Consequence of dissimilarities between the model boundary conditions and
the conditions occurring in the prototype at the location of the model boundaries.
Cauchy Number = ρv2 / K – Dimensionless number related to the ratio of inertial force to elastic
force; see Mach Number.
Cavitation Number = (p-pv)/(ρv2)– Dimensionless number defined by the difference between
pressure and vapor pressure and vapor pressure divided by the dynamic pressure; Critical Cavitation
Number describes the onset of cavitation. Also called Cavitation Index.
Densimetric Froude Number = v/(g’h)1/2 g’=g(ρ1-ρ2)/ ρ– Froude Number for situations where
the density of one fluid is significant compared to the density of the adjoining fluid.
Dimensional Analysis – Derivation of dimensionless ratios, based on the fact that all functional
relationships (for both model and prototype) must be dimensionally homogeneous.
Dimensionless Number – Physically meaningful ratio of the parameters that is dimensionless.
These dimensionless ratios are useful in determining scaling laws since a particular dimensionless
number must be the same in model and prototype to achieve similarity. Examples are the common
force ratios, such as Froude and Reynolds numbers, and ratios which are of particular physical
significance.

10
 Distortion – Conscious departure from a scaling law often necessitated by a complex set of
prototype and laboratory conditions. The term is most commonly used for geometric distortion in
which the vertical and horizontal scales of a model are different.
Dynamic Similarity – Similarity of forces.
Euler Number = (p-p0)/ (ρV2) Dimensionless number related to the ratio of inertial force to
pressure force.
Froude Number = V/(gL)1/2 – Dimensionless Number related to ratio of inertial force to
gravitation force ; the length is often the hydraulic depth making this the ratio of velocity of flow to the
speed of propagation of a small disturbance. The numerical value of one distinguishes between
subcritical and super critical flow in open channels.
Geometric Similarity – Similarity in form
Kinematic Similarity – Similarity of motion.
Laboratory Effect – Consequences of necessary laboratory simplifications or physical
constraints on the model.
Mach number = V/(K/ρ)1/2 – Dimensionless number related to the ratio of inertial force to
elastic force ; also the ratio of velocity of flow to speed of sound in an infinite medium.
Mobility number (or shields’ parameter) = ρV2*/[(ρs-ρ)gD] – Dimensionless number related to
the ratio of shear force to submerged weight of a sediment particle.
Reynolds Number = Vl/ – Dimensionless ratio related to the ratio of inertial force to viscous
force; the length may include grain size, roughness, depth of flow, and pipe diameter, resulting in
different Reynolds Numbers for different purposes. The critical Reynolds Number describes the onset
of turbulence.
Scale (or scale ratio) – Ratio of variable in model to the corresponding variable in its prototype.
Scale Effect – Consequence of non-similarity between model and prototype, resulting from the
fact that not all pertinent dimensionless numbers are the same in model and prototype. In a “real”
model, economics dictates the use of certain materials, e.g., water as a fluid. This means that fluid
density and viscosity are not correctly scaled down from prototype to model and of necessity; some
dimensionless numbers which are not the same in model and prototype, result in scale effect.
Scaling laws – Conditions that must be satisfied to achieve desired similarity between model
and prototype.
Sediment Simulant – Material used in a model to simulate prototype sediment behavior.
Shields’ parameter – See Mobility Number.
Similarity (or similitude) – Correspondence between the behavior of a model and its prototype.
Similitude – See similarity.
Tilting. – Supplying additional slope to an open channel model, e.g., to achieve a model
Reynolds numbers that would be large enough for turbulent flow or cause adequate sediment motion
in the model.
Weber Number = V2lρ/σ – Dimensionless ratio related to the ratio of inertial force to surface
tension force. The quantity is useful in analyzing thin film flows and formation of droplets and bubbles.

1.5.4 Types of Hydraulic Models

Boundary Layer Model – Reynolds Number model in which the boundary layer is carefully
reproduced which conditions in the remainder of the body of water may not be accurately
represented.
Cavitation Model – Model designed to reproduce critical cavitation numbers accurately.
Coastal Model – Model of a coastal area. Often a movable bed model to reproduce coastal
sediment transport.
Densimetric Froude Number Model – Model of gravity dominated flow usually being
densimetric Froude number scaling. Buoyant jets and stratified flows are examples.
Fixed Bed Model – Model in which the bed and side materials are nonerodible.

11
 Froude Number Model (or Gravitational Model) – Model designed to emphasize similarity of
gravitational and inertial forces (Froude Number), while other forces, such as viscous (Reynolds
Number), may not be reproduced accurately. Open channel and coastal models are of this type.
Gravitational Model – See Froude Number Model
Ground Water Model – Model designed to study the movement of ground water (seepage
model)
Harbor Model – Model to determine wave agitation, tidal flushing, etc., in a harbor.
Ice Model – Model in which formation of ice, ice conditions, or ice forces are simulated.
Movable Bed Model – Model in which the bed or side materials, or both, are erodible and are
transported in a manner similar to the prototype.
Navigation Model – Fixed-bed model to study maneuverability of vessels under currents,
waves, wind, etc. (for design of navigable waterways)
Reynolds Number Model (or Viscosity Model) – Model designed to emphasize similarity of
viscous forces (Reynolds Number), while other forces, such as gravity (Froude Number), may not be as
accurately reproduced. Pipe-flow and wind-tunnel models are of this type.
Salinity Model – Densimetric Froude Number model in which salinity distributions are
simulated.
Semimovable Bed Model (Semigrid Model) – Model with only certain portions of the bed or
side materials erodible.
Ship Model – Froude Number model to determine power requirements and maneuverability of
ships (for ship design)
Structure Model – Model to determine flow patterns, forces, and responses to forces at or
near hydraulic structures.
Tidal Model – Froude number model in which tidal water level and current fluctuations are
reproduced.
Tracer Model – Fixed-bed model in which patterns of erosion and deposition are estimated
using tracers of bed material or lightweight sediment simulants.
Upper Region Model – Froude Number model in which similarity of the main body of water is
emphasized while the boundary layer properties may not be accurately reproduced.
Viscosity Model – See Reynolds Number Model
Wave Model – Froude number model in which gravity waves are reproduced to provide the
driving mechanism for currents, impact forces, wave agitation, sediment transport, etc.

1.5.5 Types of Mathematical Models

Analytical Model – Mathematical model in which the solution of the governing equations is
obtained by mathematical analysis, as opposed to numerical manipulation.
Deterministic Model – Mathematical model in which the behavior of every variable is
completely determined by the governing equations.
Dynamic Model – A mathematical model in which time is included as an independent variable.
Empirical Model – Representation of a real system by a mathematical description based on
experimental data rather than on general physical laws.
Heuristic Model – Representation of a real system by a mathematical description based on
reasoned, but unproven argument.
Interactive Model – Numerical model which allows interaction by the modeler during
computation.
Linear Model – Mathematical model based entirely on linear equations.
Nonlinear Model – Mathematical model using one or more nonlinear equations.
Numerical Model – Mathematical model in which the governing equations are not solved
analytically. Using discrete numerical values to represent the variables involved and using arithmetic
operations. The governing equations are solved approximately.
Probabilistic Model – Mathematical model in which the behavior of one or more of the
variables is either completely or partially described by equations of probability.
12
 Semi-Empirical Model – Representation of a real system by a mathematical description based
on general physical laws but containing coefficients determined from experimental data.
Simulation Model – Mathematical model which is used with actual or synthetic input data, or
both, to produce long-term time series or predictions.
Stochastic Model – See Probabilistic Model.
Theoretical Model – Representation of a real system by a mathematical description.

1.5.6 Terms Referring to Numerical Modeling

Accuracy – The ratio of the difference between the approximate solution obtained using a
numerical model and the exact solution of the governing equations, divided by the exact solution.
Algorithm – A set of numerical steps or routines to obtain a numerical output from a numerical
input.
Aliasing – Occurrence of an apparent shift in frequency of a periodic phenomenon. It arises as
the consequence of the choice of discrete space or time sampling points to represent a continuous
process. The choice may introduce a spurious periodic solution or mask a real periodic phenomenon.
Boundary Conditions – Definition or statement of conditions or phenomena occurring at the
boundaries of the model (Slightly different definition and it is more appropriate to physical models).
Characteristics Method – Method in which the governing partial differential equations of a
mathematical model are transformed into characteristic (ordinary differential) equations and solved
on a characteristic grid.
Consistency – Quality of a discrete representation behaving as the continuum equations and
solution as the discretization is made ever finer.
Convergence – State of tending to a unique solution. A given scheme is convergent if an
increasingly finer computational grid leads to a more accurate approximation of the unique solution.
Note that a numerical method may sometimes converge on a wrong solution.
Digitization – Representation of a continuous process by numerical (digital values).
Discretization – Procedure of representing a continuous variable by discrete values at specified
points in space or time, or both.
Discretization Error – Error introduced by the discrete representation of a continuous variable.
Dissipative Scheme – Computational scheme designed to dampen or eliminate numerical
production of energy.
Explicit Scheme – Scheme in which the governing equations of a numerical model are arranged
to update the dependent variables explicitly in terms of previously known values.
Filtering – Procedure of removing or smoothing undesired perturbations.
Finite Difference Method – Method of solving the governing equations of a mathematical
model are approximated over finite intervals in the dependent variable (written in finite difference
from).
Finite Element Method – Method of solving the governing equations of a numerical model by
dividing the spatial domain into elements in each of which the solution of the governing equations is
approximated by some continuous function.
Grid (computational) – Network of points covering the space or time-space domain of a
numerical model.
Implicit Scheme – Scheme in which the governing equations of a numerical modal are
arranged to obtain solutions for the dependent variable simultaneously at several grid points. The
computed values depend not only on known values but also on other unknown neighboring values at
the surrounding grid points.
Initial Conditions – Given values of dependent variables or relationship between dependent
and independent variables at the time of start-up of the computation.
Inverse scheme – Scheme for finding values of certain independent variables at given values of
the dependent variables and time (independent variable scheme) Usually involves the finding of
physical constants from a known (often experimental) solution.

13
 Markov Chain – Sequential probabilistic process in which the probability of the state of the
variety under analysis is dependent upon its state at the previous steps in sequence. An Nth-order
Markov Chain allows dependence on up to N previous steps.
Monte Carlo Method – Technique of stochastic sampling or selection of random numbers to
generate synthetic data.
Numerical Diffusivity – Second order term introduced as a result of discretization of the
governing differential equations using finite differences, which are of the same form as diffusion
terms.
Numerical Dispersion – Numerical smearing or smoothing of a solution resulting from
numerical diffusivity.
Roundoff Error – Cumulative error introduced by rounding of the results from individual
arithmetic operations because only a finite number of digits can be retained after each operation.
Schematization – Representation of a continuum by discrete elements, e.g., dividing a real
river into reaches with constant parameters.
Scheme (numerical or computational) – Systematic program of action for solving the governing
equations of a mathematical model.
Stability (numerical or computational) – Ability of a scheme to control the propagation or
growth of small perturbations introduced in the calculations. A scheme is unstable if it allows the
growth of error so that it eventually obliterates the true solution.
Truncation Error – The error introduced by replacing the differentials of a differential equation
by finite differences using truncated Taylor series expansions.

14
 CHAPTER 2

NUMERICAL ANALYSIS

2.1 INTRODUCTION
Numerical analysis is the area of mathematics and computer science that creates, analyzes,
and implements algorithms for solving continuous mathematical problems numerically. Such problems
originate generally from real-world applications of algebra, geometry and calculus, and they involve
variables which vary continuously; these problems occur throughout the natural sciences, social
sciences, engineering, medicine, and business. During the past half-century, the growth in power and
availability of digital computers has led to an increasing use of realistic mathematical models in science
and engineering, and numerical analysis of increasing sophistication has been needed to solve these
more detailed mathematical models of the world. The formal academic area of numerical analysis
varies from quite theoretical mathematical studies to computer science issues.
With the tremendous surge of importance of using computers to carry out numerical
procedures in solving mathematical models of the scientific world, an area known as scientific
computing or computational science has taken shape during the 1980s and 1990s.
This area looks at the use of numerical analysis from a computer science perspective. It is
concerned with using the most powerful tools of numerical analysis,
computer graphics, symbolic mathematical computations, and graphical user interfaces to make it
easier for a user to set up, solve, and interpret complicated mathematical models of the real world.
Numerical algorithms are almost as old as human civilization. The Rhind Papyrus (1650 BC) of
ancient Egypt describes a root finding method for solving a simple equation. Archimedes of Syracuse
(287- 212 BC) created much new mathematics, including the “method of exhaustion” for calculating
lengths, areas, and volumes of geometric figures. When used as a method to find approximations, it is
in much the spirit of modern numerical integration; and it was an important precursor to the
development of the calculus by Isaac Newton and Gottfried Leibnitz.
A major impetus for developing numerical procedures was the invention of the calculus by
Newton and Leibnitz, as this led to accurate mathematical models for physical reality, first in the
physical sciences and eventually in the other sciences, engineering, medicine, and business. These
mathematical models cannot usually be solved explicitly, and numerical methods to obtain
approximate solutions are needed.
Another important aspect of the development of numerical methods was the creation of
logarithms by Napier (1614) and others, giving a much simpler manner of carrying out the arithmetic
operations of multiplication, division, and exponentiation.
Newton created a number of numerical methods for solving a variety of problems, and his
name is attached today to generalizations of his original ideas. Of special note is his work on root
finding and polynomial interpolation. Following Newton, many mathematics’ giants of the 18th and
19th centuries made major contributions to the numerical solution of mathematical problems.
Foremost among these are Leonhard uler (1707-1783), Joseph-Louis Lagrange (1736-1813), and Karl
Friedrich Gauss (1777-1855).

2.2 NEWTON-RAPHSON METHOD

If you've ever tried to find a root of a complicated function algebraically, you may have had
some difficulty. Using some basic concepts of calculus, we have ways of numerically evaluating roots of
complicated functions. Commonly, we use the Newton-Raphson method. This iterative process follows
a set guideline to approximate one root, considering the function, its derivative, and an initial x-value.

15
 You may remember from algebra that a root of a function is a zero of the function. This means
that at the "root", the function equals zero. As shown in Figure 2.1, the derivative of the function can
be written as

f(x)

f(xi)

f(xi)-0

0 xi+1 xi x

Figure 2.1

df / dx  f ' xi    f xi   0 /xi  xi1  (2.1)

Rearranging this equation one can obtain

x i 1  x i  f x i  / f ' x i  (2.2)

We can find these roots of a simple function such as: f(x) =x3 – 4x +1 simply by setting the
function to zero:

f(x) =x3 – 4x +1 = 0 (2.3)

The Newton-Raphson method uses an iterative process to approach one root of a function.
The specific root that the process locates depends on the initial, arbitrarily chosen x-value.

f ( xn )
xn 1  xn 
f '( xn ) (2.4)

Here, xn is the current known x-value, f(xn) represents the value of the function at xn, and f'(xn)
is the derivative (slope) at xn. xn+1 represents the next x-value that you are trying to find. Essentially,
f'(x), the derivative represents f(x)/dx (dx = delta-x). Therefore, the term f(x)/f'(x) represents a value of
dx.

16
 f ( x) f ( x)
  x
f '( x) f ( x) / x (2.5)

The more iteration that are run, the closer dx will be to zero (0). To see how this works, we will
perform the Newton-Raphson method on the function that we investigated earlier, f(x) =x3 – 4x +1.
Below are listed the values that we need to know in order to complete the process.

f(x) =x3 – 4x +1 = 0 (2.3)

f′(x) = 3x2 - 4 (2.6)

f′′(x) = 6x (2.7)

Theoretically, we could execute an infinite number of iterations to find a perfect

representation for the root of our function. However, this is a numerical method that we use to lessen
the burden of finding the root, so we do not want to do this. Therefore we will assume that the
process has worked accurately when our delta-x becomes less than 0.1. This value of precision should
be specific to each situation. A much more, or a much less, precise value may be appropriate when
using the Newton-Raphson method in class. The table below shows the execution of the process.

Table 2.1 Execution of the Process

n Xn F(x n) f'(x n) X n+1 dx
0 5,000 106,000 71,000 3,507 -
1 3,507 30,106 32,898 2,592 0,915
2 2,592 8,045 16,154 2,094 0,498
3 2,094 1,805 9,153 1,897 0,197
4 1,897 0,237 6,793 1,862 0,035
5 1,862 0,007 6,400 1,861 0,001

A graphical representation can also be very helpful. Below, you see the same function f(x) =x3 –
4x +1 (shown in blue). The process here is the same as above. In the first iteration, the red line is
tangent to the curve at x0. The slope of the tangent is the derivative at the point of tangency, and for
the first iteration is equal to 71. Dividing the value of the function at the initial x (f(5) = 106) by the
slope of the tangent (71), we find that the delta-x is equal to 1,493. Subtracting this from five (5) we
find that the new x-value is equal to 3.507. Another way of considering this is to find the root of this
tangent line. The new x-value (xn+1) will be equal to the root of the tangent to the function at the
current x-value (xn). As seen from the table, the solution is found as 1,861.

17
 Figure 2.2 Newton-Raphson Method: First iteration Initial x=5: New x= 3,507

2.3 LEAST SQUARES METHOD

Least squares method sometimes is called numerical approximation and interpolation or
regression analysis.
Assume that for a function f(x) a series of xi values and the corresponding yi=fi(i=1,..... n) are
given. The interpolation is defined on approximate computation of the f value for x≠xi. Such series of
points (fi,xi) may arise from laboratory experiments or may be given in table describing the use of an
instrument or a nomograph.
When incontrollable errors are included in the fi values (on with laboratory experiments), then
the curve passing from all points on the f-x plane may not be smooth and the cloud of points is usually
approximated by a curve passing among the point in a weighted way.

f(x) 
 
=y
  

 
 Cloud of points




x
Figure 2.3 Scatter diagram

This is usually done by the least squares method, achieving a minimum of the variance
between the smooth curve and the experimental points.

18
 This analysis is useful in finding the “best fit” of a set of experimental observations to a mathematical
relationship. If the best straight line is sought, then it is more accurate to use linear regression analysis.
The independent variable is x and this is assumed to be error free. The dependent variable is y and this
is the result of experimental observation and thus contains errors.
The general straight line is:

y= mx + c (2.7)

where m is the slope of the line and c is the point at which the line crosses the y-axis, or known as the
“y-intercept”.
Suppose we have n values of xi and corresponding observation yi

. We use the principle of “Least-square” to find the best straight line to represent the relationship
between y and x.
The principles state that the sum of the squares of the “errors” or “deviations” should be a minimum,
i.e, ∑(yi=y)2 should be a minimum. Now

(2.8)

The errors may be negative or positive so the square of errors are always positive.

y × × × × × ×
× × × × ×
× × × ×
× × × × × ×
× × × × ×
× × × x

Figure 2.4 Scatter diagram and straight line

The minimization in Eq. 2.8 is achieved by making the derivatives zero,

(2.9)

(2.10)

Differentiating the right-hand side of equation partially with respect to m and c we obtained:

(2.11)

19
 (2.12)

(2.13)

And

(2.14)

(2.15)
In which the summation are carried out from 1 to n . Solving Equations 2.13 and 2.15
(2.16)

(2.17)
+
(2.18)

(2.19)

(2.20)

The line y = mx + c is the line of regression of y on x as seen from Figure 2.5.

y=mx+c

c  tan=m

Figure 2.5 Interception point and slope of the straight line

2.4 POLYNOMIAL CURVE FITTING

There are many instances in civil engineering where it is useful to fit a mathematical
relationship, in the form of a polynomial, to experimental data. The general form of the polynomial is
taken to be;

y = a0 + a1x + a2x2 + a3x3+ ... + anxn (2.21)

The largest power of x denotes the order of the polynomial. There if n=2 in Equation (2.21), we have a
quadratic relationship between y and x;

y = a1+ a1x + a2x2 (2.22)

20
If it is assumed that the intervals x, at which observed values of y are available, are ”exact” then any
errors are confined to the observations y. If any observed values yi is known, the the corresponding
”residual” ∆i is;

∆i = yi - (a0 + a1xi + a2xi2 + ... + anxin) (2.23)

The fitting of a polynomial involves choosing values for the coefficients a0>a1....an which together
minimize the sum of the squares of all the residuals as given by equations such an Eq.(2.23). If there
are p observed values of y and if E is the sum of the squares of residuals then,

E= {yi - ( a0 + a1xi + a2xi2 + ...+ anxin)}2

= ∆i2 (2.24)

The proceeding is to differentiate Eq.(2.24) with respect to the coefficients a0, a1, ... , an is turn
and equate to zero. The resulting equations are called the “normal“ equations and their solution gives
the “best” values of the coefficient in the polynomial of chosen order n. Suppose N=2, then we have a
quadratic relationship in Eq. (2.22) and,

∆i= yi - (a0 + a1xi + a2xi2) (2.25)

E= ∆i2 = (yi - a0 + a1xi + a2xi2)2 (2.26)

Differentiating Eq.(2.26) with respect to a0, a1 and a2,

= -2yi + 2a0 + 2a1xi + 2a2xi2 = 0 (2.27)

= -2xiyi + 2a0xi + 2a1xi2 + 2a2xi3 = 0 (2.28)

= -2xi2yi + 2a0xi2 + 2a1xi3 + 2a2xi4 = 0 (2.29)

or,
│ pa0 + a1∑xi + a2∑xi2 = ∑yi
│
│a0∑xi + a1∑xi2 + a2∑xi3 = ∑xiyi
│ (2.30)
│a0∑xi2 + a1∑xi3 + a2∑xi4= ∑xi2yi

In each case the summation is taken over the p observation.

Equations (2.30) are the normal equation corresponding to a second-order polynomial.
The normal equations in the general case of a polynomial of order n are;

21
 *└ ┘= └ ┘ (2.31)

2.5 NUMERICAL INTEGRATION

2.5.1 Numerical Integration By Trapezoidal Rule
Numerical integration is defined as approximate numerical evaluation of the integral

(2.32)

Several computational schemes for numerical integration exist and can be derived from
Newton’s formula. The function error introduced can be computed from the integration of a
polynomial of degree n.
Suppose we have a function

(2.33)

Which can be represented as in Fig. 2.6 and for which we have a sequence of numerical values of y,
h1,h2,h3... at equal intervals along the x-axis.

Figure 2.6 Scheme for trapezoidal rule

The area between the curve and the x-axis is obtained approximately by adding the areas of
the trapezoid formed by drawing chords between the ordinate, i.e.

(2.34)

Thus the function (curve) is replaced by straight lines between each successive pair of
ordinates. This is the “trapezoidal” rule.

2.5.2 Numerical Integration by Simpson’s Rule

A more accurate evaluation of area may generally be obtained using Simpson’s rule. In this
rule, the function is replaced by a parabolic segment between successive groups of three ordinates as
shown in Fig.2.7. The shaded area in Fig. 2.8 is, there fore

22
 Figure 2.7 Area of parabolic segment in Simpson’s rule approximation

Figure 2.8 Scheme for Simpson’s rule

and total area between the curve and the x axis is

(2.35)

Substituting for

(2.37)

Hence in the general case Fig. 2.7, with n coordinates (n-1, intervals) the total area is

(2.38)

(2.39)

23
2.6 SOLVED PROBLEMS
2.6.1 Calculation the water depth downstream of a spillway
Problem 1: A spillway has a span width of 20 m and the maximum depth of overflow of 10 m, as
shown in the following figure. The discharge per span is calculated as 470 m 3/s. The difference
between the bottom elevations of upstream and downstream of spillway is 55 m. (a) Find the head of
water over spillway assuming the discharge coefficient of the spillway as C = 2,10. (b) Find the water
depth at the downstream of the spillway.

Figure: The spillway

Using the Bernoulli equation to determine flowrate over a dam assumes that the velocity
upstream of the dam is negligible (V1=0) and that the nappe is exposed to atmospheric pressure above
and below. Experiments have shown that the Bernoulli equation alone does not adequately predict
the flow, so empirical constants have been determined which allow better agreement between
equations and real flows. To obtain better accuracy than the Bernoulli equation alone provides, use
our weir calculations.

Discharge formula for a spillway : Q = Cxbxh01,5 →470 = 2,1x20xh01,5 , h0= 5,00 m.

The main equation used in describing open channel flow is the Bernoulli equation. This equation can
be derived by integrating the expression resulting from the application of Newton's second law to
open channel flow. Common form of Bernoulli's equation, valid at any arbitrary point along a
streamline where gravity is constant, is:

V12/(2g) + P1/γ + z1 = V22/(2g) + P2/γ + z2 = H = Constant

where:
V is the water flow speed at a point on a streamline,
g is the acceleration due to gravity,
z is the elevation of the point above a reference plane, with the positive z-direction pointing upward –
so in the direction opposite to the gravitational acceleration,
P is the pressure at the chosen point, and
γ is the specific weight the water at all points in the fluid.
H = total flow energy.
For any open channel flow (P1 = P2), and flow over a spillway (V1= 0,00; z = h), Rearranging the above
equation, yields,

H = V2/2g + h2 = (Q2/(b2h2)22g) + h2 = q2/(2gh22) + h2

where
H = total flow energy
Q = discharge over spillway
q = discharge per unit with over spillway
V2 = flow velocity at the downstream of the spillway
h2 = flow depth at the downstream of the spillway
24
 b2 = span width of the spillway
g = gravitational acceleration.
55 = (470/20)2/[2x9,81xh22] +h2
The roots of this equation can be found by using Newton-Raphson method as below:
h23 - 55 h22 + 28,147 = 0
f(h) = h23 - 55 h22 + 28,147
f´(h) = 3h22 – 110h2

Calculation table using Newton-Raphson method.

n Xn f(xn) f'(x n) X n+1 dx
0 4 -787,853 -392 1,990171 -
1 1,990171 -181,813 -207,036 1,112001 0,87817
2 1,112001 -38,4879 -118,61 0,78751 0,32449
3 0,78751 -5,47408 -84,7656 0,722931 0,064579
4 0,722931 -0,21979 -77,9545 0,720112 0,002819
5 0,720112 -0,00042 -77,6566 0,720106 5,41E-06
h2 = 0,72 m. This is the result calculated for the downstream water depth of the spillway.

Problem 2: The following equation is given. Please find the roots of this equation. To approximate the
solution of the following equation using the Newton-Raphson method, start the process by taking h=1

Solution :

For Xi=1 f(x) = -27×1+1.236=-24.764 and f `(x) = -51

Xi+1 =

Xi = 0.51

f(x) = 0.13265 - 7.0222 + 1.236 = -5.65

f `(x) = 0.7803 - 27.54 = -26.7597

25
Xi+1 =

Xi = 0.28

f(x) = 0.022 - 2.1168 + 1.236 = 1.1495

f `(x) = 0.2352 - 15.12 =-15,8726

Xi+1 =

Xi = 0.22

Xi = 0.215

2.6.2 Regression Method

Problem 1: The following table gives the results of 12 laboratory observations yi at equal intervals of
time t

t(x) 0 10 20 30 40 50 60 70 80 90 100 110 Sec

yi 3.00 3.50 3.60 3.90 3.90 4.10 4.90 5.20 5.30 5.40 6.00 6.50 cm
a) Find the line of regression of y on x
b) For t=200 find yi (extrapolation)
c) Find yi for t=85 (interpolation)

Solution: (a) For n=12

Therefore,

(b) For

(c) For

26
 FORTRAN:
LEAST SQUARE CURVE FITTING
DIMENSION X(100),F(100)
READ(5.5) N,(X(I),F(I),I=1,N)
5 FORMAT (I4/(10F 7.0))
SUM1=0.
SUM2=0.
SUM3=0.
SUM4=0.
DO 7 I=1,N
SUM1=SUM1+X(I)
SUM2=SUM2+F(I)
SUM3=SUM3+X(I)××2
SUM4=SUM4+X(I)*F(I)
7 CONTINUE
A=(SUM 4*N-SUM1*SUM2)/(N*SUM3-SUM1××2)
B=(SUM3*SUM2-SUM1*SUM4)/(N*SUM3-SUM1××2)
WRITE(6.8) A,B
8 FORMAT (2F (0.4))
STOP
END

Problem 2: Find the relationship between volume and water level in the reservoir for the given data.

Water level, h Reservoir volume, V

(m) (106 m3)
0 0
5 2,25
10 8,05
15 17,75
20 30,5
25 46
30 65,7
35 89,5

The characteristic relationship can be written in the form of

V = a*hb

Take the logarithmic transformation of both sides of the equation

log(V) = log(a) + b*log(h)

log(V) = y, log(a) = c, log(h) = x

Then, equation turns out a simple linear equation as:

y = c + b*x

The regression equation coefficients can be found by using following formulations.

27
 𝑛 ∑ 𝑥𝑖 𝑦𝑖 − ∑ 𝑥𝑖 ∑ 𝑦𝑖 ∑ 𝑥𝑖2 ∑ 𝑦𝑖 − ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑥𝑖
𝑏= 𝑐=
𝑛 ∑ 𝑥𝑖2 − (∑ 𝑥𝑖 )2 𝑛 ∑ 𝑥𝑖2 − (∑ 𝑥𝑖 )2

h 5 10 15 20 25 30 35 Σ
V 2250000 8.05*106 17.75*106 30.5*106 46*106 65.7*106 89.5*106 259.75*106
x=logh 0,699 1 1,176 1,301 1,398 1,477 1,544 8,595
y=logV 6,352 6,906 7,249 7,484 7,663 7,818 7,952 51,424
x*y 4,440 6,906 8,526 9,737 10,712 11,547 12,278 64,146
x^2 0,488 1 1,383 1,693 1,954 2,182 2,384 11,085

∑ 𝑥𝑖 = 8,595 ∑ 𝑦𝑖 = 51.424 , ∑ 𝑥𝑖 𝑦𝑖 = 64.146 , ∑ 𝑥𝑖2 = 11.085

2
(∑ 𝑥𝑖 ) = 73.878 , ∑ 𝑥𝑖 ∑ 𝑦𝑖 = 441.997,

∑ 𝑥𝑖2 ∑ 𝑦𝑖 = 3799.102

7 ∗ 64.146 − 441.997
𝑏= = 1.892
7 ∗ 11.085 − 73.878
3799.102 − 64.146 ∗ 8.595
𝑐= = 5.023
7 ∗ 11.085 − 73.878
b = 1.892

c = 5.023 = Log (a)  a = 105506.487

The relationship between h and V is obtained as:

𝑽 = 𝟏𝟎𝟓𝟓𝟎𝟔. 𝟒𝟖𝟕𝒉𝟏.𝟖𝟗𝟐

2.6.3 Polynomial Curve Fitting

Problem 1: Fit a second degree polynomial to the following data

X 0 1 2 3 4 5
Y 5 15 26 27 30 35

y = a0 + a1x + a2x2
p = 6, n = 2;
∑x i = 15, ∑xi2=55, ∑xi3=225, ∑xi4=979,
∑yi=138, ∑xiyi=443, ∑xi2yi=1717

The normal equations are, therefore

6 15 55 𝑎0 138
[15 55 225] [𝑎1 ] = [ 443 ]
55 225 979 𝑎2 1717

From which
28
 a0 = 5.61 6a0 + 15a1 + 55a2 = 138
a1 =10.7 15a0 + 55a1+ 225a2 = 443
a2 = -1.02 55a0+ 225a1 + 979a2 =1717

Hence the polynomial is,

y=5.61+10.7x-1.02x2

2.6.4 Numerical Integration

Problem 1: Find analytical and numerical solution for the integral of the following equation.

(a) Analytical solution

6
3

= (360+144+6) - (45+36+3) = 510 – 84 = 426

(b) Numerical solution by hand (Simpson’s rule)

N = 7 number of ordinates
n-1= number of intervals

I=a/3 [sum of end ordinates + 4(sum of even ordinates)+2(sum of other odd ordinates)]
x1=3.0 , x2=3.5 x3=4.0 y1=4.5 x5=5 x6=5.5 x7= 6.0
f1(3,0)=5*3*3+8*3+1=70
f2(3,5)=5*3,5*3,5+8*3,5+1=90,5
f3(4,0)=5*4*4+8*4+1=113
f4(4,5)=5*4,5*4,5+8*4,5+1=138,25
f5(5,0)=5*5*5+8*5+1=166
f6(5,5)=5*5,5*5,5+8*5,5+1=196,25
f7(6,0)=5*6*6+8*6+1=229.

I=426,166.

Analytical solution is 426. The difference between analytical and numerical solution is 0,116. The
numerical solution is obtained by hand. If we increase the number ordinates and use a small computer
program, the difference between analytical and numerical solution should be less than 0,166.

Problem 2: Find the emptying time of a reservoir that is required to lower the water level from h1=50
m to h2=40 m. The reservoir volume-elevation relationship is given as V = 8000 h2,7 , where V is the
volume in m3 and H is the water level elevation in m. The outlet works is located at the elevation of h’=

29
10 m, its diameter is 2,5 m and discharge coefficient, c, is 0,63 (assumed constant during in the
variation of reservoir elevation. The inflow through the reservoir is constant as I = Qi =10 m3/s.

Solution: The continuity equation is written as

𝑑𝑉
= 𝑂 − 𝐼 = 𝑄𝑜 − 𝑄𝑖
𝑑𝑡
in which, V is volume, t is time, O is outflow, and I is inflow
𝑡2 ℎ2
𝑑𝑉
∫ 𝑑𝑡 = ∫
𝑄𝑜 − 𝑄𝑖
𝑡1 ℎ1

𝑉 = 𝑎ℎ𝑏
𝑑𝑉 = 𝑎 ∗ 𝑏 ∗ ℎ𝑏−1 𝑑ℎ
𝑉 = 8000ℎ2,7
𝑑𝑉 = 21600ℎ0,7 𝑑ℎ
𝐷2
𝑄𝑜 = 𝐶𝑜 ∗ 𝜋 ∗ ∗ √2𝑔(ℎ − ℎ′ )
4
2,52
𝑄𝑜 = 0,63 ∗ 𝜋 ∗ ∗ √2𝑔(ℎ − 10)
4

𝑡 ℎ 𝑑𝑉 ℎ 𝑎∗𝑏∗ℎ 𝑏−1
Find 𝑇 = ∫𝑡 2 𝑑𝑡 = ∫ℎ 2 = ∫ℎ 2 2 0.5
1 𝑄 1 𝑐 − 𝑄𝑔 1 𝐶 𝑜 ∗𝜋∗𝑟 ∗(2𝑔∗ℎ) − 𝑄𝑖

𝑡2 50
21600 ∗ ℎ0.7 𝑑ℎ
𝑇 = ∫ 𝑑𝑡 = ∫
𝑡1 40 13.698 ∗ (ℎ − 10)0.5 − 25

To solve this equation we will use Simpson’s formula. The number of ordinates, n, is taken as 5, so we
have 4 intervals.

ℎ2 − ℎ1 10
𝑎= = = 2.5
𝑛 − 1 4

21600 ∗ ℎ0.7
𝑓(ℎ) =
13.698 ∗ (ℎ − 10)0.5 − 25

ℎ1 = 40.0 , ℎ2 = 42.5 , ℎ3 = 45.0 , ℎ4 = 47.5 , ℎ5 = 50.0

i hi f(hi)
1 40 5710,7
2 42,5 5614,5
3 45 5536,3
4 47,5 5472,1
5 50 5418,9
30
 2.5
𝑇= ( 5710,7 + 5418,9 + 4(5614,5 + 5472,1) + 2(5536,3))
3
T = 55456,9 seconds = 15,4 hours

2.6.5 Numerical Integration by FORTRAN Programming

Problem: Compute the volume of water collected during one year in a reservoir at the upstream of a
dam built on a river where discharges (measured every month) are given in the following table by
means of trapezoidal rule. What must be the storage capacity of the reservoir to meet the need for a
continuous water supply if 42 m3/s?
For the computation of the reservoir capacity the curve of accumulated volume with time is evaluated
numerically.
The integration is performed with h=30 days, for all months, and the results of the interpretation are
printed continuously.

Month 1 2 3 4 5 6 7 8 9 10 11 12
Discharge
x106m3/day 8 6 5 5 4 1 1 1 5 6 7 7

The FORTRAN code is:

NUMERICAL INTEGRATION
DIMENSION F(100)
READ(5,5)N,H,(F(I),I=1,N)
5 FORMAT (I4,F7.0)/(10F7.10))
SUM=0
NN=N-1
DO 10 I=1, NN
SUM=SUM+(F(I)+F(I+1))/2*/t
WRITE(6,15) SUM
15 FORMAT (E14.6)
10 CONTINUE
STOP
END

2.6.6 Numerical Integration by BASIC Programming

Problem: Write a Basic code to calculate the area under a curve usind trapezoidal and Simpson’s rules

PROGRAM (BASIC)
10 PRINT “AREA UNDER CURVE BY TRAPEZOIDAL AND SIMPSON’S RULES”
20 PRINT “TYPE IN NUMBER OF ORDINATES”
30 INPUT N
40 IF (N-1)/2=INT(N/2) THEN 60
50 PRINT “N MUST BE ODD;THAT NUMBER WAS EVEN”:GO TO 20
60 DIM H(N)
70 PRINT “TYPE IN INTERNAL”
80 INPUT L
90 PRINT “TYPE IN ORDINATES; EACH FOLLOWED BY ‘RETURN’”
100 FOR I=1 TO N: INPUT H(I)
110 FOR I=2 TO N-1 :H(I) =2*H(I): NEXT
120 FOR I= 1 TO N: A=A+H(I):NEXT

31
130 A=A*L/2
140 PRINT: PRINT “ AREA BY TRAPEZOIDAL RULE=”;A
150 A=0
160 FOR I=2 TO N-1 STEP 2:H(I)=2*H(I): NEXT170 FOR I=1 TO N:A=A+H(I): NEXT
180 A=A*L/3
190 PRINT : PRINT “AREA BY SIMPSON’S RULE=”; A
200 END

2.6.7 Numerical Integration by MATLAB

Problem: Write a m-file to compute the integral of f(x)=sin2(x)+x2 between the boundries [0  ] by
numerical integration for 100 segments using trapezoidal rule and Simpson’s rules .
𝐻
Hint: 𝐼 = (𝑓(1) + 2 ∗ ∑𝑛𝑖=2 𝑓(𝑖) + 𝑓(𝑛 + 1)) n=number of segments, H=interval width
2

a=0;
b=pi;
n=100;
h=(b-a)/n;
x=linspace(a,b,n+1);
f=@(x)sin(x).^2+x.^2;
fpm=feval(f,x);
Itrap = (fpm(1) + 2*sum(fpm(2:n)) + fpm(n+1)) * (h/2)
Isimp = (fpm(1) + 4*sum(fpm(2:2:n)) + 2*sum(fpm(3:2:n-1)) + fpm(n+1)) * (h/3)

32