Muskingum-Cunge Matlab

HW-Chapter 9 Naoki Mizukami
9-1. Detention basin with 700ft x 700ft with 4:1 side slope. Top of the berm is at elevation 9 ft.
Inflow hydrograph is a triangular shape with a peak discharge of 125 cfs at a time of 8 hr and
with a base time of 24 hr. The detention basin is initially empty. Route the inflow hydrograph
through the detention basin and determine the peak outflow and stage.
Stage-storage relationship is computed as follows

z z
V ( z )   A( z ) dz    700  4 z  2  dz 
 8 z  700 3  700 3
0 0
24
Stage-outflow relationship is given and summarized along with storage in the table below.
Table. Stage-storage-outflow relationships based on given outlet structures and detention basin geometry
Table. reservoir routing using storage-indication method

The peak outflow is 36.7 cfs at 19 hr. The maximum storage is 8 ft at 19 hr.
140
120
100
80
Q [cfs]
60
40
20
0
0 5 10 15 20 25
Time [hr]
Figure. Inflow and outflow hydrograph
10
6
Z [ft]
0
0 5 10 15 20 25
Time [hr]
Figure. Temporal change in stagee, Z
9-3. Rout the inflow hydrograph given using Muskingum method with Δt = 2 hr, Δx = 25000 ft,
θ = 3.6 hr, X = 0.2. Plot the inflow and outflow hydrographs and determine the percent reduction
in the inflow peak as well as the travel time of the peak
Figure. Inflow and outflow hydrograph

9-16. Using Muskingum-Cunge method to route the inflow hydrograph and river in example 9.6
for S0 = 0.001.
Figure. Inflow hydrograph
Average Q to compute normal depth, y0.

1T 500 * 7  ( 4500  500) * 6 / 2
Q  
T 0
Q (t ) dt 
7
 2500 cfs
Using average Q, find average velocity, flow area with manning’s equation
5 5
Kn A 3
1 Kn  y 0   b  m  y 0   3 1
Q  S 2    S 2
 
0 0
n 23 n 2
P b  2 y0 1  m 2 3
5
1.49  y 0  100  2 y 0   3
1
2500    0.001 2
 
2
0.025
100  2 y 0 1  2 2 3
Solving for y0, y0 = 4.65 ft,
2500
V  Q/A  = 4.91 ft/s
4.65  (100  2  4.65)
So wave celerity is calculated
5 5
ck  V   4.91  8.19
3 3
Δt = 0.5 hr is chosen so that 5 discretized points are in rising part of inflow hydrograph
1 Q  1 2500 
x   c k t     8.19  0.5  60  60    8659 ft
2  BS 0 c k  2  118 .6  0 .001  8.19 
So select Δx = 6000 ft
Compute Courant number, Cn as

t 0.5  60  60
C n  ck  8.19   2.46
x 6000
One more parameter, Muskingum weighting parameter, X is computed as
1 Q  1 2500 
X  1    1    0.29
2 BS 0 c k x  2  118 .6  0.001  8.19  6000 
Using parameters, X and Cn, M-C coefficients are computed as

0.5C n  X 0.5  2.46  0.29
C0    0.488
1  X  0.5C n 1  0.29  0.5  2.46
0.5C n  X 0.5  2.46  0.29
C1    0.785
1  X  0.5C n 1  0.29  0.5  2.46
1  0.5C n  X 1  0.5  2.46  0.29
C2    0.274
1  X  0.5C n 1  0.29  0.5  2.46
Outflow is computed using the Muskingum routing equations written by
Qi j 11  C 0  Qi j 1  C1  Qi j  C 2  Qi j 1
j
where the Qi refers to the flow rate at space j and time i. The result of routing is as follows
Table. Muskingum-Cunge routing with constant parameter.

Figure. Inflow and outflow hydrographs from Muskingum-Cunge routing.

9-18. Write the code for Muskingum-Cunge routing using the 4 point variable parameter method
and apply it to Example 9.6.
Matlab Code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Muskingum-Cunge with variable coefficient method
%Coefficients X, and ck (wave celerity) are dynamic in time and space due to
%varying Q, resulting in variable C0, C1, and C2.
%
% brief description of 4 point variable parameter method
%the 4 point method computes an average Q using inflow and outflow at the first step and next step and
%use the average Q to estimate M-C parameters. This computation is performed every time outflow at the
%next time step is computed. Since outflow at the next time step is unknown, three point Q are used for
%first guess for average Q. Iteration is necessary to make average Q converged at a certain error
%tolerance.
%
%
clear all
close all
%Reading inflow hydrograph

data1 = load('inhydro18.dat');
t = data1(:,1) % time [hr]
Q(:,1) = data1(:,2) %inflow [cfs]
dt = (t(2)-t(1))*60*60;% time interval [hr]
%River geometry
L = 18000; %river reach length [ft]
S0 = 0.0005; %bed slope [ft/ft]
b = 100; %bottom width [ft]
m = 2; %side slope
n = 0.025; %manning n
%some constant
Kh = 1.49; %manning equation
%define routing reach length, dx [ft]

dx = 9000;
%compute distance x where hydrographs are computed
x(1) = 0; %begining of river reach
for i = 2:L/dx+1
x(i) = x(i-1)+dx;
end
%initial condition at x
for j = 1:length(x)
Q(1,j) = Q(1,1);
end
%compute Muskingum-Cunge routing

%go thru all time steps and spatial step
for i = 2:length(t)
for j = 2:length(x)
%initial guess for Q using three points (inflow at 1st time step and 2nd time step
%and outflow at 1st time step.
Qave = (Q(i-1,j-1)+Q(i,j-1)+Q(i-1,j))/3;
%iteration till 4 point average Q are converged at 0.1 of error tolerance
while 1
%Define manning equation
func = @(y) Kh/n*(y*(b+2*y))^(5/3)/(b+2*y*sqrt(1+m^2))^(2/3)*sqrt(S0)-Qave;
%compute depth, y0, based on manning equation. use bisect method with 0.01 error tolerance
y0 = fbisect(func,0,10,0.01);
%compute flow area, A, velocity, V, and wave celerity, ck, top width, TOP, given depth y0
A = y0*(b+2*y0);
V= Qave/A;
ck = 5/3*V;
TOP = b+2*m*y0;
%Compute coefficients for Muskingum-Cunge equation

X = 0.5*(1-Qave/TOP/S0/ck/dx);
Cn = ck*dt/dx;
C0=(0.5*Cn-X)/(1-X+0.5*Cn);
C1=(0.5*Cn+X)/(1-X+0.5*Cn);
C2=(1-0.5*Cn-X)/(1-X+0.5*Cn);
%Compute outflow at next time step

Q(i,j)=C0*Q(i,j-1)+C1*Q(i-1,j-1)+C2*Q(i-1,j);
%Take average of 4 points

Qave1=(Q(i-1,j-1)+Q(i,j-1)+Q(i-1,j)+Q(i,j))/4;
%If difference between current 4 point average and previous one is less than 0.1
%leave iteration loop, otherwise use new average to compute new estimate of
%outflow at next time step
if abs(Qave-Qave1)< 0.1; break; end
Qave = Qave1;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of code
Results (applied this code to example 9.6)
Table. Muskingum-Cunge routing (inflow and routed hydrograph at x =9000 and x = 18000)
5000
infow
4500 x=9000 ft
x=18000 ft
4000
3500
3000
Q [cfs]
2500
2000
1500
1000
500
0
0 1 2 3 4 5 6 7 8
Time [hr]
Figure. Inflow and outflow hydrographs

Muskingum-Cunge Matlab

Uploaded by

Copyright:

Available Formats

Muskingum-Cunge Matlab

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Muskingum-Cunge Matlab

Uploaded by

Copyright:

Available Formats

HW-Chapter 9 Naoki Mizukami

Stage-storage relationship is computed as follows

Table. reservoir routing using storage-indication method

Figure. Inflow and outflow hydrograph

Figure. Inflow hydrograph

Average Q to compute normal depth, y0.

Compute Courant number, Cn as

Using parameters, X and Cn, M-C coefficients are computed as

Table. Muskingum-Cunge routing with constant parameter.

Figure. Inflow and outflow hydrographs from Muskingum-Cunge routing.

%Reading inflow hydrograph

%define routing reach length, dx [ft]

%compute Muskingum-Cunge routing

%Compute coefficients for Muskingum-Cunge equation

%Compute outflow at next time step

%Take average of 4 points

You might also like