CE 572-01: Shallow Water Hydrodynamics
HOME WORK 1
Professor:Dr. Abul B.M.Baki
By Cássio Rampinelli
September, 23th, 2018
OBS: This script was done in R programming language
__________________________________________________________________________________________________________________________
QUESTION 1
Consider the steady, uniform flow in a straight prismatic channel with a compound cross-section, as shown in Figure 1. The side walls in the main channel and floodplains are vertical. The Manning`s n coefficients in the left floodplain, main channel and right floodplain are 0.04, 0.025 and 0.035, respectively. The longitudinal slope of the channel is 0.001.
Figure 1: Reference Cross Section for Question 1
(a) Find an adequate method to determine the flow conveyance of this channel.
Answer
First of all, we should define the concept of conveyance.
To start, let`s consider the Manning-Strickler Equation as follows: \[ Q = \frac{A}{n}R^{(2/3)}_HS^{(1/2)}_0 \] Where:
Q is the flow rate (m³/s);
n is the Manning`s Roughness coefficient (\(m^{-1/3}s\));
\(R_H\) is the hydraulic radius (m);
\(A\) is the flow are (m²);
\(S_0\) is the Energy Grade Line Slope
* Obs: in uniform flow we assume that \(S_0\) is equals to the bed slope*
In Manning-Strickler Equation both A and \(R_h\) are function of the water level inside the channel which is expressed by h. Then, we can restate the Manning`s Equation joining in one term, A and \(R_h\) that depend on h.
\[ Q = K(h)S^{(1/2)}_0 \] K is called the Conveyance Function that represents the ability of the cross section to convey flow
\[ K(h)=\frac{A(h)}{n}R_H(h)^{(2/3)}_H \]
with \(A(h)\) and \(R_H(h)\) determined by the cross section geometry.So, let`s develop these functions.
\(A(h)\) is the area of the cross section (Figure 1) and can be expressed as a piecewise function.As the Manning`s coefficient is different depending on the region of the cross section (see Figure 1), then, conveniently, we will express the functioin of the area considering the three partitions shown in Figure 1 (AL=Left Area; AM=Middle Area; AR=Right Area).:
Para: \(0\small\leq h \small\leq 5\) \[ A(h)= 10h \]
Para: \(h > 5\) \[ A(h)= ((h-5)\cdot10)_{AL} +(10\cdot h)_{AM} + ((h-5)\cdot10)_{AR} \]
\(R(h)\) is the flow area divided by the weted perimeter (\(P_m\)). Then, we, also, should use a peacewise function, as follows:
Para: \(0\small\leq h \small\leq 5\) \[ R_H(h)=\frac{A(h)}{P_w} \] \[ R_H(h)=\frac{10h}{(10+2h)} \] Para: \(h > 5\) \[ R_H(h)=\Bigg( \frac{[(h-5)\cdot10]}{(10+h)} \Bigg)_{AL}+\Bigg( \frac{(h\cdot10)}{(10+2\cdot 5)} \Bigg)_{AM}+\Bigg( \frac{[(h-5)\cdot10]}{(10+h)} \Bigg)_{AR} \]
Finally, replacing \(R(h)\) and \(A(h)\) in \(K(h)\)
Para: \(0\small\leq h \small\leq 5\) \[K(h)=\frac{10h }{n_{AM}}\Bigg(\frac{10h}{(10+2h)}\Bigg) R^{(2/3)}_H \] Para: \(h > 5\) \[ K(h)=\Bigg[ \frac{[(h-5)\cdot10]}{n_{AL}}\Bigg( \frac{[(h-5)\cdot10]}{(10+h)} \Bigg)^{(2/3)} \Bigg]_{AL}+\Bigg[ \frac{(h\cdot10)}{n_{AM}}\Bigg( \frac{(h\cdot10)}{20} \Bigg)^{(2/3)} \Bigg]_{AM}+\Bigg[ \frac{[(h-5)\cdot10]}{n_{AR}}\Bigg( \frac{[(h-5)\cdot10]}{(10+h)} \Bigg)^{(2/3)} \Bigg]_{AR} \] Then, we can plot the Conveyance and check how \(K(h)\) varies as a function of “h”.
To do so, a R code for Conveyance function was created and applyed to a set of “h”s values. The results were plotted next.
The function was created so that it can be applied to any geometry similar to that one presented in Figure 1. So, the geometric parameters of the section were reset as shown in Figure 2.
Figure 2: Reference Cross with generalized geometric parameters
The R code is shown bellow, followed by the plot.
#Clear R memory
rm(list=ls())
#Define the conveyance function and the input parameters
conveyance<-function(nAL,nAM,nAR,bAL,bAM,bAR,H,h)
{
#code of the piecewise function
if (h>H)
{
k=(((h-H)*bAL)/nAL)*(((h-H)*bAL)/(bAL+h))^(2/3)+((h*bAM)/nAM)*((h*bAM)/(bAM+2*H))^(2/3)+(((h-H)*bAR)/nAR)*(((h-H)*bAR)/(bAR+h))^(2/3)
}
else
k=((bAM*h)/nAM)*((bAM*h)/(bAM+2*h))^(2/3)
return(k)
}
#Define the values of h in which Conveyance function will be applied
h=seq(from=0.01, to=10, by=0.1)
#Define Question 1 n and geometric parameters of the cross section
#Manning coefficient for left part of the cross section
nAL=0.04
#Manning coefficient for the middle part of the cross section
nAM=0.025
#Manning coefficient for the right part of the cross section
nAR=0.035
#width of the left part of the cross section
bAL=10
#width of the middle part of the cross section
bAM=10
#width of the right part of the cross section
bAR=10
#height of the inferior part of the cross section
H=5
#initial value of the vector that will record the conveyance function vaules.
k=0
#Applying the conveyance function in all values of the vector h and recorde the result in the vector k
for(i in 1:length(h) )
{
k[i]=conveyance(nAL,nAM,nAR,bAL,bAM,bAR,H,h[i])
}Conveyance Function
#plottin the results
plot(h,k,type='l',col="red",xlab="h [m]",ylab="K(h) [m³/s]",axes=FALSE)
title(main="Conveyance Function")
#Axis editing
axis(1, at = seq(from=0,to=10,by=2), labels = c(0,2,4,6,8,10))
rug(x = seq(from=1,to=9,by=2), ticksize = -0.01, side = 1)
abline(v=seq(from=0,to=10,by=2), col="lightgray", lty="dotted")
axis(2, at = seq(from=0,to=14000,by=2000), labels = seq(from=0,to=14000,by=2000))
abline(h=seq(from=0,to=14000,by=2000), col="lightgray", lty="dotted")
Now, we can come back to the Manning-Strickler equation and rewrite the flow equation adding the slope of the channel given in the problem.
\[ Q = K(h)S^{(1/2)}_0 \] \[ Q = K(h)(0.001)^{(1/2)} \] \[ Q = K(h)0.031623 \] This is the conveyance function times a coefficient.
Finally, we can replot the graphs showing the relation between discharge and flow height and vice versa. That is the rating curve of the cross section
Rating Curve
par(mfrow=c(1,2))
#plottin the results for rating curve with Q in y axis and h in x axis
plot(h,k*sqrt(0.001),type='l',col="red",xlab="h [m]",ylab="Q(h) [m³/s]",axes=FALSE)
title(main="Q(h) X h")
#X axis editing
axis(1, at = seq(from=0,to=10,by=2), labels = c(0,2,4,6,8,10))
rug(x = seq(from=1,to=9,by=2), ticksize = -0.01, side = 1)
abline(v=seq(from=0,to=10,by=2), col="lightgray", lty="dotted")
#Y axis editing
axis(2, at = seq(from=0,to=450,by=50), labels = seq(from=0,to=450,by=50))
abline(h=seq(from=0,to=450,by=50), col="lightgray", lty="dotted")
box()
#plottin the results for rating curve with h in y axis and Q in x axis
plot(k*sqrt(0.001),h,type='l',col="red",xlab="Q(h) [m³/s]",ylab="h [m]",axes=FALSE)
title(main="h x Q(h)")
#Y axis editing
axis(1, at = seq(from=0,to=450,by=50), labels = seq(from=0,to=450,by=50))
abline(v=seq(from=0,to=450,by=50), col="lightgray", lty="dotted")
axis(2, at = seq(from=0,to=10,by=2), labels = c(0,2,4,6,8,10))
rug(x = seq(from=1,to=9,by=2), ticksize = -0.01, side = 1)
abline(h=seq(from=0,to=10,by=2), col="lightgray", lty="dotted")
box()
(b) If the water depth in the main channel is 7.5 m, what is the flow discharge in this channel?
Answer
To answer this question we can use the rating curve graph presented and take the value directly from there. Or, we can call the the Conveyance Function finding \(K(h=7.5)\) and substitute the answer in the rating curve.
Once we had defined the \(K(h)\) function, we call the function setting the parameters with the R code as follows:
conveyance(nAL=0.04,nAM=0.025,nAR=0.035,bAL=10,bAM=10,bAR=10,H=5,h=7.5)## [1] 8939.966Then, \(K(h=7.5)\)=8939.97 m³/s, this brings to:
\[ Q = 8939.97\cdot 0.031623 \] \[ Q = 282.71 m³/s \]
(c) If the flow discharge is 400 m³/s, what is the water depth in this channel?
Answer
This question follows the opposite process of the previous one. To solve that, we also can use the rating curve graph presented later and take the value directly from there. Or, we can recall the rating curve (Conveyance Function times the bed slope) and find for which “h”, Q=400 m³/s.
If we apply the opposite function for the rating curve it brings to: Q(h=8.78 m)=400 m³/s.
Then, h=8.78 m is the water depth of the channel related with the discharge of 400 m³/s.
The R code to calculate the inverse function is presented bellow.
#Clear R memory
rm(list=ls())
library(GoFKernel)## Loading required package: KernSmooth## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009 #Define the rating curve based on conveyance function
f<-function(h)
{
#Define Question 1 n and geometric parameters of the cross section
#Manning coefficient for left part of the cross section
nAL=0.04
#Manning coefficient for the middle part of the cross section
nAM=0.025
#Manning coefficient for the right part of the cross section
nAR=0.035
#width of the left part of the cross section
bAL=10
#width of the middle part of the cross section
bAM=10
#width of the right part of the cross section
bAR=10
#height of the inferior part of the cross section
H=5
#initialize conveyance
k=0
#Slope
S=0.001
#code of the piecewise function for conveyance
if (h>H)
{
k=(((h-H)*bAL)/nAL)*(((h-H)*bAL)/(bAL+h))^(2/3)+((h*bAM)/nAM)*((h*bAM)/(bAM+2*H))^(2/3)+(((h-H)*bAR)/nAR)*(((h-H)*bAR)/(bAR+h))^(2/3)
}
else k=((bAM*h)/nAM)*((bAM*h)/(bAM+2*h))^(2/3);q=k*sqrt(S)
return(q)
}
#Define the invers function of the flow
f.inv <- inverse(f,lower=0,upper=500)
#calculates which h produces Q=400 m³/s
f.inv(400)## [1] 8.781599
QUESTION 2
Consider the compound channel shown in the figure below. The coefficient of Manning for the main channel and its banks is 0.015 \(m^{-1/3}s\). The coefficient of Manning for the floodplains and their banks is 0.035 \(m^{-1/3}s\). The uniform flow depth in the channel is 4.0 m.
a) What is the flow discharge?
Answer
The flow discharge is 362.83 m³/s. (OBS: We assumed logitudinal slope \(S_{0}=0.001\))
To answer this question, a R code was developed to a generic section, with generic geometric parameters (Figure 4). A function was created to return the flow rate for any flow depth (h). Then, the given geometric parameters (Figure 3) were used on the created function to calculate the flow rate.
Figure 4: Generic Reference Cross Section for Question 2
The R code is presented below as the discharge result for the given parameters.
#Clear R memory
rm(list=ls())
#Define the function to calculate discharge given cross section geometric paramters
flowRate<-function(nAL,nAM,nAR,bAL,bAM,bAR,H,m1,m2,S0,h)
{
#code to caculate flow rate given a generic section
if (h>H)
{
x=m2*(h-H)
bx=m1*H
k=((((bAM+(bAM+2*bx))*h/2)+(bAM+2*bx)*(h-H))/nAM)*((((bAM+(bAM+2*bx))*h/2)+(bAM+2*bx)*(h-H))/(bAM+(sqrt(bx^2+H^2))*2))^(2/3)+((bAL*(h-H)+x*(h-H)/2)/nAL*((bAL*(h-H)+x*(h-H)/2)/(bAL+sqrt(x^2+(h-H)^2)))^(2/3))+((bAR*(h-H)+x*(h-H)/2)/nAR*((bAR*(h-H)+x*(h-H)/2)/(bAR+sqrt(x^2+(h-H)^2)))^(2/3))
flow=k*sqrt(S0)
}
else
{
bx=m1*h
k=(((bAM+(bAM+2*bx))*h/2)/nAM)*(((bAM+(bAM+2*bx))*h/2)/(bAM+(sqrt(bx^2+h^2))*2))^(2/3)
flow=k*sqrt(S0)
}
return(flow)
}
#Define Question 2 n and geometric parameters of the cross section
#Main channel bank slope coefficient
m1=2
#Floodplains main sope coefficient
m2=3
#Height of the main channel
H=2.5
#Manning coefficient for left part of the cross section
nAL=0.035
#Manning coefficient for the middle part of the cross section
nAM=0.015
#Manning coefficient for the right part of the cross section
nAR=0.035
#width of the left part of the cross section
bAL=10
#width of the middle part of the cross section
bAM=5
#width of the right part of the cross section
bAR=10
#initial value of the vector that will record the conveyance function vaules.
#Longitudinal slope
S0=0.001
#call the flowRate
flowRate(nAL,nAM,nAR,bAL,bAM,bAR,H,m1,m2,S0,h=4)## [1] 362.8313
b) Determine the energy coefficient. Is the use of the energy coefficient justified when computing gradually varied flow in this channel?
Answer
The energy coefficient is also known as Coriolis coefficient (\(\alpha\)) and is defined as:
\[ \alpha=\frac{\iint v³dA}{\overline{V}^{3}A} \] where
\(\alpha\) is the Coriolis coefficient
\(v\) is the velocity distribution in the channel (m/s)
\(\overline{V}\) is the mean flow velocity of the channel (m/s)
\(A\) is the total flow area (m²)
For this question, we will consider a discrete case where the total area will be divided into 3 parts. The velocity field variation will be simplified by 3 punctual velocities that represent each partition as shown in Figure 5
Hence, we can rewrite the Coriolis coefficient as:
\[ \alpha=\frac{\sum_{i=1}^{n} v_{i}^3\cdot A_i}{\overline{V}^{3}A} \] \[ \alpha=\frac{v_{1}^3\cdot A_1+v_{2}^3\cdot A_3+v_{3}^3\cdot A_3 }{\overline{V}^{3}(A_{1}+A_{2}+A_{3})} \]
Finally, to comput \(\alpha\) we need to calculate \(A_{1}\), \(A_{2}\), \(A_{3}\), \(v_{1}\), \(v_{2}\), \(v_{3}\), and \(\overline{V}\)
Given the data, we found \(\alpha\)=1.81 which represents a considerable magnitude justifying the use of the energy coefficient.
The R code to calculate each component of \(\alpha\) expression is presented bellow with the results.
#Define Question 2 n and geometric parameters of the cross section
#Main channel bank slope coefficient
m1=2
#Floodplains main sope coefficient
m2=3
#Height of the main channel
H=2.5
#Manning coefficient for left part of the cross section
nAL=0.035
#Manning coefficient for the middle part of the cross section
nAM=0.015
#Manning coefficient for the right part of the cross section
nAR=0.035
#width of the left part of the cross section
bAL=10
#width of the middle part of the cross section
bAM=5
#width of the right part of the cross section
bAR=10
#initial value of the vector that will record the conveyance function vaules.
#Longitudinal slope
S0=0.001
#flow depth
h=4
#Auxiliary geometric variables
x=m2*(h-H)
bx=m1*H
#Area 1
A1=(bAL*(h-H)+x*(h-H)/2)
A1## [1] 18.375#Area2
A2=(((bAM+(bAM+2*bx))*h/2)+(bAM+2*bx)*(h-H))
A2## [1] 62.5#Area3
A3=(bAR*(h-H)+x*(h-H)/2)
A3## [1] 18.375#Total Flow Area
A=A1+A2+A3
#discharge for Part 1
Q1=(((bAL*(h-H)+x*(h-H)/2)/nAL*((bAL*(h-H)+x*(h-H)/2)/(bAL+sqrt(x^2+(h-H)^2)))^(2/3)))*sqrt(S0)
#velocity for Part 1
v1=Q1/A1
v1## [1] 1.046369#discharge for Part 2
Q2=(((((bAM+(bAM+2*bx))*h/2)+(bAM+2*bx)*(h-H))/nAM)*((((bAM+(bAM+2*bx))*h/2)+(bAM+2*bx)*(h-H))/(bAM+(sqrt(bx^2+H^2))*2))^(2/3))*sqrt(S0)
#velocity for Part 2
v2=Q2/A2
v2## [1] 5.190035#discharge for Part 3
Q3=(((bAR*(h-H)+x*(h-H)/2)/nAR*((bAR*(h-H)+x*(h-H)/2)/(bAR+sqrt(x^2+(h-H)^2)))^(2/3)))*sqrt(S0)
#velocity for Part 3
v3=Q3/A3
v3## [1] 1.046369#channel mean velocity
V=(Q1+Q2+Q3)/(A1+A2+A3)
V## [1] 3.655731#energy coefficient
Alpha=((v1^3)*A1+(v2^3)*A2+(v3^3)*A3)/((A1+A2+A3)*(V^3))
Alpha## [1] 1.810613c) Compute the Froude number for the flow.
Answer
Considering the given data, \(F_r\)= 0.78.
The Froude number \(F_r\) is given by: \[ F_r=\frac{V}{\sqrt{gD_h}}=\frac{3.655731}{\sqrt{9.81\cdot2.255682}}=0.7771427 \] where
\(F_r\) is the Froude number
\(V\) is the mean flow velocity in the channel (m/s)
\(g\) is the gravity acceleration (9.81 m²/s)
\(D_h\) is the hydraulic depth (m)
\(D_h\) is defined as: \[ D_h=\frac{A}{T}=\frac{99.25}{44}=2.255682 \] where
\(A\) is the flow area (m²)
\(T\) Top width (m)
#top widht
Tt=bAM+bAL+bAR+2*x+2*bx
Tt## [1] 44#Hydraulic depth
Dh=A/Tt
Dh## [1] 2.255682#Froude number
Fr=V/(sqrt(9.81*Dh))
Fr## [1] 0.7771427d) Compute the bed shear stress on the channel bed
Answer
The shear stress (\(\tau\)) is defined as:
\[ \tau=\gamma\cdot R_h\cdot S_0 \] where
\(\gamma\) is the water specific weight (9810 N/m³)
\(R_h\) is the hydraulic radius (m)
\(S_0\) is the longitudinal sope
\[ \tau=9810\cdot 2.173334\cdot 0.01=21.32 N/m² \] The R code is presented below.
#wetted perimeter
wt=bAL+bAM+bAR+2*sqrt(x^2+(h-H)^2)+2*sqrt(bx^2+H^2)
wt## [1] 45.66717#Hydraulic Radius
Rh=A/wt
Rh## [1] 2.173334#Tau
Tau=9810*Rh*S0
Tau## [1] 21.3204#Froude number
Fr=V/(sqrt(9.81*Dh))
Fr## [1] 0.7771427QUESTION 3
A trapezoidal channel has a bottom width of b=5.0 m and side slopes of m=1.0. Calculate and plot the specific energy curves for discharges Q=2.0, 3.0, 4.0 and 5.0 m³/s.
Plot the specific energy curves in the same figure and assess the trend of critical depth varying with flow discharge.
Answer
To this question, we will consider Figure 6 bellow, as reference.
Figure 6: Cross section for this Question 3
The Specific Energy (\(H_s\)) as a function of h, for Q constant is given by: \[ H_s=h+\frac{Q²}{2gA²} \] where
\(h\) is the flow depth (m)
\(Q\) is the discharge (m³/s)
\(g\) is the gravity acceleration (9.81 m²/s)
\(A\) is the flow area (m²)
Considering the given cross section, we can express the flow area (A) as a function of h
\[ A(h)= \frac{[b+(b+2\cdot h\cdot m) \cdot h]}{2} \]
Rewriting the Specific Energy
\[ H_s=h+\frac{Q²}{2g\Bigg(\frac{[b+(b+2\cdot h\cdot m) \cdot h]}{2})\Bigg)^2} \]
Finally, we can set constant values to Q and varies h computing Hs to obtain the Specific Energy Curve. A R code was developed to create those curves considering a constant value for each given discharge.
The minimum values of each specific energy curve were taken as a reference and a linear regression model was adjusted resulting in the dashed black line that represents the trend of critical depth with respect to flow discharge. It is noticed in the plots that the critical depth is directly related to the flow discharge, hence, it tends to increase with the rising of the flow and vice versa.
The plots and the R code are presented below, including the critical depths for the given curves and the linear regression model adjusted.
#Clear R memory
rm(list=ls())
#defining Hs function for the first discharge
Hs1<-function(h)
{
Q=B=m=A=H=0
#Define the discharge
Q=2
#Define the cross section geometric parameters
#Base widht
B=5
#cross section side slope
m=1
A=(B+(B+2*h*m))*(h/2)
H=h+Q^2/(2*9.81*A^2)
return(H)
}
#defining Hs function for the second discharge
Hs2<-function(h)
{
Q=B=m=A=H=0
#Define the discharge
Q=3
#Define the cross section geometric parameters
#Base widht
B=5
#cross section side slope
m=1
A=(B+(B+2*h*m))*(h/2)
H=h+Q^2/(2*9.81*A^2)
return(H)
}
#defining Hs function for the third discharge
Hs3<-function(h)
{
Q=B=m=A=H=0
#Define the discharge
Q=4
#Define the cross section geometric parameters
#Base widht
B=5
#cross section side slope
m=1
A=(B+(B+2*h*m))*(h/2)
H=h+Q^2/(2*9.81*A^2)
return(H)
}
#defining Hs function for the fourth discharge
Hs4<-function(h)
{
Q=B=m=A=H=0
#Define the discharge
Q=5
#Define the cross section geometric parameters
#Base widht
B=5
#cross section side slope
m=1
A=(B+(B+2*h*m))*(h/2)
H=h+Q^2/(2*9.81*A^2)
return(H)
}
#initial variables to record the values for the specific energy curves
H1=0
H2=0
H3=0
H4=0
h=h2=h3=h4=0
h=seq(from=0.04,to=2,by=0.01)
h2=seq(from=0.06,to=2,by=0.01)
h3=seq(from=0.08,to=2,by=0.01)
h4=seq(from=0.1,to=2,by=0.01)
#Calculating the specific energy curves H x h for each discharge
for(i in 1:length(h))
{
H1[i]=Hs1(h[i])
}
for(i in 1:length(h2))
{
H=0
H2[i]=Hs2(h2[i])
}
for(i in 1:length(h3))
{
H=0
H3[i]=Hs3(h3[i])
}
for(i in 1:length(h4))
{
H=0
H4[i]=Hs4(h4[i])
}
#Critical Energy for Q1=2m³/s
min(H1)## [1] 0.3683484#Critical Depth for Q1=2m³/s
h[which.min(H1)]## [1] 0.25#Critical Energy for Q1=3m³/s
min(H2)## [1] 0.4782727#Critical Depth for Q1=3m³/s
h2[which.min(H2)]## [1] 0.33#Critical Energy for Q1=4m³/s
min(H3)## [1] 0.5745501#Critical Depth for Q1=4m³/s
h3[which.min(H3)]## [1] 0.39#Critical Energy for Q1=5m³/s
min(H4)## [1] 0.6618473#Critical Depth for Q1=5m³/s
h4[which.min(H4)]## [1] 0.45x=c(h[which.min(H1)],h2[which.min(H2)],h3[which.min(H3)],h4[which.min(H4)])
y=c(min(H1),min(H2),min(H3),min(H4))
dataValues=cbind(x,y)
dataValues=data.frame(dataValues)
linearMod<-lm(x~y,data=dataValues)
print(linearMod)##
## Call:
## lm(formula = x ~ y, data = dataValues)
##
## Coefficients:
## (Intercept) y
## 0.002797 0.676333criticalH<-function (h)
{
H=0.6778*h+0.0015
return(H)
}
Hc=0
hc=seq(from=0,to=5,by=0.1)
for(i in 1:length(hc))
{
Hc[i]=criticalH(hc[i])
}Specific Energy Plots
#ploting the results
plot(H2,h2,type="b",pch=8,col="black",lwd=1,xlab="Hs [m]",ylab="h [m]",ylim=c(0,2),xlim=c(0,5),axes=F,yaxs="i",xaxs="i", main="Specific Energy Curves")
lines(H1,h,type="b",lwd=1,pch=15,col="red",lty=1,cex=0.5)
lines(H3,h3,type="b",lwd=1,pch=13,col="blue",lty=1,cex=0.5)
lines(H4,h4,type="b",lwd=1,pch=17,col="darkgreen",lty=1,cex=0.5)
lines(hc,Hc,lwd=1,pch=0,col="black",lty=2,cex=0.5)
#x axis editing
axis(1, at = seq(from=0,to=5,by=1), labels = seq(from=0,to=5,by=1),lty=1)
#x grid
abline(v=seq(from=0,to=5,by=0.5), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=2,by=0.5), labels = seq(from=0,to=2,by=0.5))
#y grid
abline(h=seq(from=0,to=2,by=0.1), col="lightgray", lty="dotted")
#legend
legend(3, 1.4, legend=c("Q=2 m³/s","Q=3 m³/s", "Q=4 m³/s", "Q=5 m³/s","Critical Depth"),col=c("red","black","blue","darkgreen","black"), lty=c(1,1,1,1,2),pch=c(15,8,13,17,NA),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
QUESTION 4
A trapezoidal channel having a bottom width of b=10.0 m and side slopes of m=1.0, conveys a discharge of Q=40.0 m³/s. The channel has a constant bed slope of S0=0.001. The Manning friction coefficient for the channel is n=0.025m-1/3s.
Answer
Initially, we need to verify what is the normal depth \(h_0\) for the channel considering the given data. To do so, we apply the Manning equation considering the provided information. The normal depth obtained was \(h_0\)=1.99 m. Figure 7 and the R code employed to do the calculation are presented as follows.
Figure 7: Cross section for this Question 4
#Define Mannings Equation to a trapezoidal channel
maningEq<-function(h){
bAM=10
m1=1
S0=0.001
nAM=0.025
bx=m1*h
k=(((bAM+(bAM+2*bx))*h/2)/nAM)*(((bAM+(bAM+2*bx))*h/2)/(bAM+(sqrt(bx^2+h^2))*2))^(2/3)
Q=k*sqrt(S0)
return (Q)
}
#Define the invers function of the flow
f.inv <- inverse(maningEq,lower=0,upper=5)
#calculates which h produces Q=40 m³/s
f.inv(40)## [1] 1.989529
To answer the next questions is also helpful computes the critical depth. To do that, the specific energy curve will be derivated as was done in Question 3. The critical depth (\(h_c\)) obtained was \(h_c\)=1.13 m. The R code and the plot is presented below.
Hs<-function(h)
{
Q=B=m=A=H=0
#Define the discharge
Q=40
#Define the cross section geometric parameters
#Base widht
B=10
#cross section side slope
m=1
A=(B+(B+2*h*m))*(h/2)
H=h+Q^2/(2*9.81*A^2)
return(H)
}
h=0
h=seq(from=0.5,to=3,by=0.01)
H=0
for(i in 1:length(h))
{
H[i]=Hs(h[i])
}
#Critical Energy for Q1=40m³/s
min(H)## [1] 1.645554#Critical Depth for Q1=40m³/s
xc=c(0,min(H))
yc=c(h[which.min(H)],h[which.min(H)])
xc1=c(min(H),min(H))
yc1=c(0,h[which.min(H)])Specific Energy Plots
#ploting the results
plot(H,h,col="blue",lwd=2.5,xlab="Hs [m]",ylab="h [m]",ylim=c(0,3),xlim=c(1,3),axes=F,xaxs="i",yaxs="i", main="Specific Energy Curve for Q=40m³/s",lty=1,type="l")
lines(xc,yc,type="l",lwd=1,pch=15,col="red",lty=1,cex=0.5)
lines(xc1,yc1,type="l",lwd=1,pch=15,col="red",lty=1,cex=0.5)
#x axis editing
axis(1, at = seq(from=0,to=3,by=1), labels = seq(from=0,to=3,by=1),lty=1)
#x grid
abline(v=seq(from=0,to=3,by=0.5), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=3,by=0.5), labels = seq(from=0,to=3,by=0.5))
#y grid
abline(h=seq(from=0,to=3,by=0.1), col="lightgray", lty="dotted")
text(x=1.3, y = 1.3, labels = "hc=1.13 m")
text(x=1.8, y = 0.4, labels = "Hc=1.65 m")
box()
a) Determine the type of water surface profile to be expected if the channel terminates by a sudden drop of the bed.
Answer
Taking into account the verification done above, we can assume that the normal slope is less than the critical slope and that the flow depth is higher than the critical depth. In this case, we are in a Mild slope. Hence, if the channel terminates by a sudden drop of the bed, we expect that the water surface profile will decrease producing a profile M2, as shown in Figure 8.
Figure 8: Possible M2 profiles(Chow, Ven Te, 1959)
b) Determine the type of water surface profile to be expected if the outlet water level is 10 m.
Answer
Taking into account the same assumptions referred to Question 4(b), however, considering that the outlet water level is 10 m, hence, higher than the normal depth, we expect a water surface profile type M1, as shown in Figure 9.
Figure 9: Possible M1 profiles(Chow, Ven Te, 1959)
c) Make a MatLab code to calculate the water surface profiles for a reach length of 20000 m under the above two conditions.
Answer
To compute the water surface profile we will apply the Direct Step Method (or Step-by-Step Method). According to Chow (1959), this method was first suggested by the Polish engineer Charnomskil in 1914 and then by Husted in 1924.
This idea of this method is breaking up the water profile into different segments and then finding the length of each one of them by accounting for the energy grade line slope at the midpoint of each segment. To do so, we need to know the water depth at the beginning and at the end of the reach we want to compute the water surface profile (usually that information are given or can be calculated by the data available). To start, we calculate the maximum water level difference between the first and the last cross-sections of the reach. Then, the difference is divided by the number of segments which will be used to discretize the reach, this will result in the incremental value of water depth (\(\small \Delta\)y). Starting from the beginning /or from the end of the reach, this incremental value can be added to or subtracted from the knowning water depth in the boundary section. Finally, to be able to draw the water profile we also need the length increment (\(\small \Delta\)x) associated with the water level variation so that we can have the coordinates (x,y) of the water surface profile. For the derivation of the describe the method, consider Figure 10 below, with two cross-sections of a channel reach.
For the derivation of the described the method, consider Figure 10 below, with two cross-sections of a channel reach.
Figure 10: A channel reach for the derivation of the Direct Step Method
By the principle of the energy conservation, we can say that the total energy (H) in section 1 is equivalent to the total energy in section 2. After some algebra we can obtain how long (\(\small \Delta\)x) should be.
\[ H_1=H_2 \]
\[ S_0 \cdot\Delta x +y_1+\alpha_1\frac{v_1^2}{2g}=S_f \cdot\Delta x +y_2+\alpha_2\frac{v_2^2}{2g} \]
Where
\(S_0\) is the bed slope
\(S_f\) is the energy grade line slope
\(\Delta\)x is the length increment (m)
\(y_1\) and \(y_2\) are the flow depth in each section (m)
\(v_1\) and \(v_2\) are the flow velocity in each section (m/s)
\(\alpha_1\) and \(\alpha_2\) are the Coriolis coefficients
g is the gravity acceleration (m/s²)
Assuming: \[ E_1= y_1+\alpha_1\frac{v_1^2}{2g}; E_2=y_2+\alpha_2\frac{v_2^2}{2g} \] Then we can rewrite the former equation as: \[ S_0 \cdot\Delta x +E_1=S_f \cdot\Delta x +E_2 \]
Which brings to: \[ \Delta x =\frac{E_2-E_1}{S_0-S_f}=\frac{\Delta E}{S_0-S_f} \]
Considering that the energy grade line can be curved between the two cross-section, in this method, as a simplification, we will assume that \(S_f\) is the mean of the grade line slope between the two cross sections \(\overline{S_f}\)
\[ \Delta x =\frac{\Delta E}{S_0-\overline{S_f}} \]
All the terms of the right-hand side of the equation we can calculate from the given data. It is also useful considering the following:
1-For subcritical flow: proceed the calculation by upstream;
2-For supercritical flow: proceed the calculation by downstream;
3-Consider the \(\Delta\)y increment less or equal to 0.5 m per segment/step
4-The geometry of the channel should be prismatic and the same for all cross sections
Considering the general presentation of the Direct Standard Method we can start calculating the water surface profiles for the two conditions (a) and (b). The calculations will be presented in R language code instead of Matlab.
For the first condition, we are interested in calculating the water surface profile in a channel that terminates by a sudden drop of the bed. Figure 11 presents an overview of the problem.
Figure 11: Overview of the the first condition
Before starting, it is necessary to recall the normal depth (\(y_0\)=1.99) and the critical depth (\(y_c\)=1.13) that were previously calculated. Hence, y*=0.79 m. Following is presented the R code to the direct step method and the plots of the water profile. Two plots are presented. One with a large overview of the 20 km of the reach, and the second one with a zoom in the first 800 m of the reach. Visually, is possible to infer that the most significant change of the water profile occurs in the first 500 m. Due to this point, the differences between the normal depth and water surface level are less than 12 cm.
#Clear R memory
rm(list=ls())
###########
###Setting discretization parameters
############
yBegin=1.13 #Define the water depth for the intial downstream section (critical depth)
yFinal=1.99#Define the water depth for the final upstream section (normal depth)
nOfSteps=65080 #Number of steps for calculation
deltaY=(yFinal-yBegin)/nOfSteps #Define y increment
###########
###Setting geometric parameters and input data
###########
b=10#botton channel widht
m=1# side channel slope
S0=0.001# longitudinal channel slope
Q=40# flow discharge
n=0.025 #manning roughness
yc=1.13 #critical depth
yn=1.99 #normal depth
###########
###Initializing state variables
###########
y=0#record flow depth at botton reference
y[1]=yBegin
Tt=0#record top width
Tt[1]=b+2*m*y[1]
P=0#record the wetted perimeter
P[1]=b+2*sqrt((y[1])^2+(m*y[1])^2)
A=0#record flow area
A[1]=(b+Tt[1])*y[1]/2
Rh=0#record hydraulic radius
Rh[1]=A[1]/P[1]
v=0#record flow velocity
v[1]=Q/A[1]
E=0#record especific energy
E[1]=y[1]+(v[1])^2/(2*9.81)
deltaE=0#record delta energy
Sf=0#record energy slope
##Obs: this variable will be initialized after function definitions
AvgSf=0#record average slope
deltaX=0#comput delta x
xBed=0#record x bed coordinate
yBed=0#record y bed elevation
wselv=0#record water surface elevation
ycelev=0#record critical elevation
ynelev=0#record normal elevation
##########
####SlopeEnergyFunction
###########
SEF<-function(Qd,nd,Ad,Rhd)
{
SE=((Qd*nd)/(Ad*(Rhd)^(2/3)))^2
return(SE)
}
##########
####deltaXFunction
###########
deltaXFunction<-function(deltaEd,S0d,Sfmed)
{
deltaXd=(deltaEd)/(S0d-Sfmed)
return(deltaXd)
}
####Initializing values
Sf[1]=SEF(Q,n,A[1],Rh[1])
#Applying the Direct Step Method
for(i in 1:nOfSteps)
{
#updates all variables for the first step calculations
y[i+1]=y[i]+deltaY
Tt[i+1]=b+2*m*y[i+1]
P[i+1]=b+2*sqrt((y[i+1])^2+(m*y[i+1])^2)
A[i+1]=(b+Tt[i+1])*y[i+1]/2
Rh[i+1]=A[i+1]/P[i+1]
v[i+1]=Q/A[i+1]
E[i+1]=y[i+1]+(v[i+1])^2/(2*9.81)
deltaE[i]=E[i+1]-E[i]
Sf[i+1]=SEF(Q,n,A[i+1],Rh[i+1])
AvgSf[i]=mean(c(Sf[i],Sf[i+1]))
deltaX[i]=deltaXFunction(deltaE[i],S0,AvgSf[i])
}
#####
###Initialize variables for plot routine
#####
xBed[1]=0#record x bed coordinate
xBed[2]=4*yc #experimental constant
yBed[1]=0#record y bed elevation
yBed[2]=xBed[2]*S0#update position of the bed for the second section
wselv[1]=0.7*yc#record water surface elevation
wselv[2]=yc+yBed[2]
ycelev[1]=yc#record critical elevation
ycelev[2]=yc+yBed[2]
ynelev[1]=yn#record normal elevation
ynelev[2]=yn+yBed[2]
########
###Plot routine
#######
for(i in 1:nOfSteps)
{
xBed[i+2]=xBed[i+1]-deltaX[i]
yBed[i+2]=xBed[i+2]*S0
wselv[i+2]=y[i+1]+yBed[i+2]
ycelev[i+2]=yc+yBed[i+2]
ynelev[i+2]=yn+yBed[i+2]
}Water Surface Profile Plots
par(mfrow=c(1,1))
#ploting the results
plot(xBed,yBed,col="black",lwd=3,xlab="Distance [m]",ylab="Elevation [m]",ylim=c(-2,25),xlim=c(-1000,25000),axes=F,xaxs="i",yaxs="i", main="Water Surface Elevation",lty=1,type="l")
lines(xBed,wselv,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed,ycelev,type="l",lwd=1,col="red",lty=2)
lines(xBed,ynelev,type="l",lwd=1,col="darkgreen",lty=3)
lines(c(0,0),c(0,-2), col="black", lty=1,lwd=3)
#x axis editing
axis(1, at = seq(from=0,to=22000,by=2000), labels = seq(from=0,to=22000,by=2000),lty=1)
#x grid
abline(v=seq(from=0,to=22000,by=1000), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=25,by=5), labels = seq(from=0,to=25,by=5))
#y grid
abline(h=seq(from=0,to=25,by=1), col="lightgray", lty="dotted")
#legend
legend(13000, 7.5, legend=c("Bed Elevation","Water Surface Elevation", "Normal Elevation", "Critical Elevation"),col=c("black","blue","darkgreen","red"), lty=c(1,1,3,2),lwd=c(3,2,1,1),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
#ploting the results
plot(xBed,yBed,col="black",lwd=3,xlab="Distance [m]",ylab="Elevation [m]",ylim=c(-0.8,3),xlim=c(-20,820),axes=F,xaxs="i",yaxs="i", main="Water Surface Elevation-Zoom in the first 800 m",lty=1,type="l")
lines(xBed,wselv,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed,ycelev,type="l",lwd=1,col="red",lty=2)
lines(xBed,ynelev,type="l",lwd=1,col="darkgreen",lty=3)
lines(c(0,0),c(0,-0.8), col="black", lty=1,lwd=3)
#x axis editing
axis(1, at = seq(from=0,to=820,by=100), labels = seq(from=0,to=820,by=100),lty=1)
#x grid
abline(v=seq(from=0,to=820,by=50), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=3,by=0.5), labels = seq(from=0,to=3,by=0.5))
#y grid
abline(h=seq(from=0,to=3,by=0.1), col="lightgray", lty="dotted")
#legend
legend(520, 0.35, legend=c("Bed Elevation","Water Surface Elevation", "Normal Elevation", "Critical Elevation"),col=c("black","blue","darkgreen","red"), lty=c(1,1,3,2),lwd=c(3,2,1,1),cex=0.9,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
For the second condition, we have at the starting point, normal depth (\(y_0\)=1.99) and at the end of the reach flow depth (\(y_final\)=10 m). As shown in the plots at around 9500 the backward effect is insignificant.
#Clear R memory
rm(list=ls())
###########
###Setting discretization parameters
############
yBegin=10 #Define the water depth for the intial downstream section (critical depth)
yFinal=1.99#Define the water depth for the final upstream section (normal depth)
nOfSteps=100000 #Number of steps for calculation
deltaY=(yFinal-yBegin)/nOfSteps #Define y increment
###########
###Setting geometric parameters and input data
###########
b=10#botton channel widht
m=1# side channel slope
S0=0.001# longitudinal channel slope
Q=40# flow discharge
n=0.025 #manning roughness
yc=1.13 #critical depth
yn=1.99 #normal depth
###########
###Initializing state variables
###########
y=0#record flow depth at botton reference
y[1]=yBegin
Tt=0#record top width
Tt[1]=b+2*m*y[1]
P=0#record the wetted perimeter
P[1]=b+2*sqrt((y[1])^2+(m*y[1])^2)
A=0#record flow area
A[1]=(b+Tt[1])*y[1]/2
Rh=0#record hydraulic radius
Rh[1]=A[1]/P[1]
v=0#record flow velocity
v[1]=Q/A[1]
E=0#record especific energy
E[1]=y[1]+(v[1])^2/(2*9.81)
deltaE=0#record delta energy
Sf=0#record energy slope
##Obs: this variable will be initialized after function definitions
AvgSf=0#record average slope
deltaX=0#comput delta x
xBed=0#record x bed coordinate
yBed=0#record y bed elevation
wselv=0#record water surface elevation
ycelev=0#record critical elevation
ynelev=0#record normal elevation
##########
####SlopeEnergyFunction
###########
SEF<-function(Qd,nd,Ad,Rhd)
{
SE=((Qd*nd)/(Ad*(Rhd)^(2/3)))^2
return(SE)
}
##########
####deltaXFunction
###########
deltaXFunction<-function(deltaEd,S0d,Sfmed)
{
deltaXd=(deltaEd)/(S0d-Sfmed)
return(deltaXd)
}
####Initializing values
Sf[1]=SEF(Q,n,A[1],Rh[1])
#Applying the Direct Step Method
for(i in 1:nOfSteps)
{
#updates all variables for the first step calculations
y[i+1]=y[i]+deltaY
Tt[i+1]=b+2*m*y[i+1]
P[i+1]=b+2*sqrt((y[i+1])^2+(m*y[i+1])^2)
A[i+1]=(b+Tt[i+1])*y[i+1]/2
Rh[i+1]=A[i+1]/P[i+1]
v[i+1]=Q/A[i+1]
E[i+1]=y[i+1]+(v[i+1])^2/(2*9.81)
deltaE[i]=E[i+1]-E[i]
Sf[i+1]=SEF(Q,n,A[i+1],Rh[i+1])
AvgSf[i]=mean(c(Sf[i],Sf[i+1]))
deltaX[i]=deltaXFunction(deltaE[i],S0,AvgSf[i])
}
#####
###Initialize variables for plot routine
#####
xBed[1]=0#record x bed coordinate
#xBed[2]=4*yc #experimental constant
yBed[1]=0#record y bed elevation
#yBed[2]=xBed[2]*S0#update position of the bed for the second section
wselv[1]=yBegin#record water surface elevation
#wselv[2]=yc+yBed[2]
ycelev[1]=yc#record critical elevation
#ycelev[2]=yc+yBed[2]
ynelev[1]=yn#record normal elevation
#ynelev[2]=yn+yBed[2]
########
###Plot routine
#######
for(i in 1:nOfSteps)
{
xBed[i+1]=xBed[i]-deltaX[i]
yBed[i+1]=xBed[i+1]*S0
wselv[i+1]=y[i+1]+yBed[i+1]
ycelev[i+1]=yc+yBed[i+1]
ynelev[i+1]=yn+yBed[i+1]
}
i=1
xBed2=xBed[length(xBed)]
yBed2=yBed[length(yBed)]
ycelev2=ycelev[length(ycelev)]
ynelev2=ynelev[length(ynelev)]
wselv2=wselv[length(wselv)]
while(xBed2[i]<=20000)
{
xBed2[i+1]=xBed2[i]+1000
yBed2[i+1]=xBed2[i+1]*S0
ycelev2[i+1]=yc+yBed2[i+1]
ynelev2[i+1]=yn+yBed2[i+1]
wselv2[i+1]=ynelev2[i+1]
i=i+1
}Water Surface Profile Plots
par(mfrow=c(1,1))
#ploting the results
plot(xBed,yBed,col="black",lwd=3,xlab="Distance [m]",ylab="Elevation [m]",ylim=c(-2,25),xlim=c(-100,22000),axes=F,xaxs="i",yaxs="i", main="Water Surface Elevation",lty=1,type="l")
lines(xBed,wselv,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed,ycelev,type="l",lwd=1,col="red",lty=2)
lines(xBed,ynelev,type="l",lwd=1,col="darkgreen",lty=3)
lines(xBed2,wselv2,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed2,ycelev2,type="l",lwd=1,col="red",lty=2)
lines(xBed2,ynelev2,type="l",lwd=1,col="darkgreen",lty=3)
lines(xBed2,yBed2,type="l",lwd=3,col="black",lty=1)
lines(c(0,0),c(0,-2), col="black", lty=1,lwd=3)
#x axis editing
axis(1, at = seq(from=0,to=22000,by=2000), labels = seq(from=0,to=22000,by=2000),lty=1)
#x grid
abline(v=seq(from=0,to=22000,by=1000), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=25,by=5), labels = seq(from=0,to=25,by=5))
#y grid
abline(h=seq(from=0,to=25,by=1), col="lightgray", lty="dotted")
#legend
legend(12000, 10, legend=c("Bed Elevation","Water Surface Elevation", "Normal Elevation", "Critical Elevation"),col=c("black","blue","darkgreen","red"), lty=c(1,1,3,2),lwd=c(3,2,1,1),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
#ploting the results
plot(xBed,yBed,col="black",lwd=3,xlab="Distance [m]",ylab="Elevation [m]",ylim=c(-2,15),xlim=c(-100,10000),axes=F,xaxs="i",yaxs="i", main="Water Surface Elevation-Zoom in the first 10000 m",lty=1,type="l")
lines(xBed,wselv,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed,ycelev,type="l",lwd=1,col="red",lty=2)
lines(xBed,ynelev,type="l",lwd=1,col="darkgreen",lty=3)
lines(xBed2,wselv2,type="l",lwd=2.0,col="blue",lty=1)
lines(xBed2,ycelev2,type="l",lwd=1,col="red",lty=2)
lines(xBed2,ynelev2,type="l",lwd=1,col="darkgreen",lty=3)
lines(xBed2,yBed2,type="l",lwd=3,col="black",lty=1)
lines(c(0,0),c(0,-2), col="black", lty=1,lwd=3)
#x axis editing
axis(1, at = seq(from=0,to=22000,by=2000), labels = seq(from=0,to=22000,by=2000),lty=1)
#x grid
abline(v=seq(from=0,to=22000,by=1000), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=25,by=5), labels = seq(from=0,to=25,by=5))
#y grid
abline(h=seq(from=0,to=25,by=1), col="lightgray", lty="dotted")
#legend
legend(6000, 5, legend=c("Bed Elevation","Water Surface Elevation", "Normal Elevation", "Critical Elevation"),col=c("black","blue","darkgreen","red"), lty=c(1,1,3,2),lwd=c(3,2,1,1),cex=0.9,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
QUESTION 5
In a long rectangular channel, stagnant water at a depth of h1=5.0 m is stored behind the gate. The channel has a smooth bed of negligible slope. The gate is suddenly and completely opened by an accident.
Determine the water surface profiles at different instants, t=50, 100, 150 and 200 sec after the accident, in the following two cases.
a) The channel downstream the gate is dry.
Answer
First of all, let`s consider the noncoservation form Saint Venant Equations, as presented below:
\[ V\frac{\partial y}{\partial x}+y\frac{\partial V}{\partial x}+\frac{\partial y}{\partial t} = 0\ (Continuity\ Equation) \]
\[ \frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}+g\frac{\partial y}{\partial x} -g(S_0 -S_f) = 0\ (Momentum Equation) \]
For a frictionless and horizontal channel of prismatic and rectangular cross-section the Momentum equation can be rewritten without the gravity and friction force terms.
\[ \frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}+g\frac{\partial y}{\partial x} = 0\ \]
Where
\(V\) is the flow velocity
\(y\) is the flow depth
\(x\) is the distance in regards to some reference
\(g\) is the gravity acceleration
\(S_0\) is the longitudinal bed slope
\(S_f\) is the energy grade line slope
Ritter addressed that equations to solve an ideal dry-bed dam break case for a rectangular and horizontal channel. That is, actually, the same case of this question, with the difference that instead of a dam we have a gate, that is suddenly removed. Ritter found an analytical solution to the problem resulting in the following analytical equations that allow us determining the flow depth in any position and in any time, starting from the gate vertical axis (or dam axes) going to the downstream direction. Figure 12 and 13 give an overview of the problem.
\[ y(x,t)=\frac{1}{9g}\Bigg(2\sqrt {gh_0}-\frac{x}{t}\Bigg)^{2} \]
\[ V(x,t)=\frac{2}{3}\Bigg(\sqrt {gh_0}+\frac{x}{t}\Bigg) \]
Orgaz & Chanson (2017) presented an interesting review of Ritter study, including some original Ritter manuscripts at the end of the paper.
The R code and the plots for the water surface profiales at instants t=50, 100, 150 and 200 are presented as follows. An annimated plot is also presented.
#Clear R memory
rm(list=ls())
###########
###Setting geometric parameters and input data
###########
h0=5#reservoir height
#distance discretization
x=seq(from=-2000, to=3000,by=1)
#time discretization
t=seq(from=1,to=200,by=1)
#matrix to record the profiles along x for each t
matrixResult=matrix(0,nrow=length(t),ncol=length(x))
##########
####Ritter function for water surface profile
###########
yFunction<-function(t,x)
{
Y=(1/(9*9.81))*(2*sqrt(9.81*h0)-(x/t))^(2)
return(Y)
}
########
###Computation routine
#######
for(i in 1:length(t))
{
for(j in 1:length(x))
{
matrixResult[i,j]=yFunction(t[i],x[j])
}
}Water Surface Profile Plots
par(mfrow=c(1,1))
#ploting the results
plot(x[600:4590],matrixResult[200,600:4590],col="darkgreen",lwd=1.5,xlab="Distance [m]",ylab="Elevation [m]",axes=T,xaxs="i",yaxs="i", main="Water Surface Profiles for different time steps",lty=1,type="l",ylim=c(0,5),pch=8)
lines(x[948:3880],matrixResult[150,948:3880],type="l",lwd=1.5,col="blue",lty=2,pch=15)
lines(x[1297:3300],matrixResult[100,1297:3300],type="l",lwd=1.5,col="red",lty=3,pch=16)
lines(x[1651:2660],matrixResult[50,1651:2660],type="l",lwd=1.5,col="purple",lty=4,pch=17)
#x axis editing
axis(1, at = seq(from=-2000,to=3000,by=1000), labels = seq(from=-2000,to=3000,by=1000),lty=1)
#x grid
abline(v=seq(from=-2000,to=3000,by=500), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=5,by=1), labels = seq(from=0,to=5,by=1))
#y grid
abline(h=seq(from=0,to=5,by=0.5), col="lightgray", lty="dotted")
#legend
legend(1300, 4.5, legend=c("t= 200 s","t=150 s", "t=100 s", "t=50 s"),col=c("darkgreen","blue","red","purple"), lty=c(1,2,3,4),lwd=c(1.5,1.5,1.5,1.5),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
Figure 14: Water Surface Profile along the time
b) The channel downstream of the gate is filled with stagnant water at a depth of h0=0.5 m
Answer
In this case, before the sudden removal of the gate consider that the downstream water depth is \(h_0\). If the negative wave in the previous dry bed case is valid, it needs to intersect with the downstream water surface at C, how presented in the Figure 15 below. From the previous analysis, the negative wave should have a velocity at C, but, the downstream water is still in height \(h_0\). The, this connection is not physically right. (Baki/Wu, classes notes)
Figure 15: Discontinuity in water profile (Baki/Wu, classes notes)
Then, Figure 16 below is the valid water surface profile (Baki/Wu classes notes). There is a discontinuity at C, where a positive wave occurs. Upstream of point C, the water depth is \(h_2\).
Figure 16: Discontinuity in water profile (Baki/Wu, classes notes)
After some algebra, from the below two equations \(h_2\) can be determine by trial and error, then \(c_tc\).
\[ c_{tc}=\sqrt{gh_0}\sqrt{\frac{h_2}{2h_0}\Bigg(1+\frac{h_2}{h_0}\Bigg)} \] \[ c_{tc}=\frac{(-2\sqrt{gh_2}+2\sqrt{gh_1})h_2}{h_2-h_0} \]
A R code was developed to solve the system, converting the problem in finding the root of a equation. The code and the plots are presented as follows.
#Clear R memory
rm(list=ls())
###########
###Setting geometric parameters and input data
###########
h0=0.5#reservoir height
h1=5
h2=1
#distance discretization
x=seq(from=-2000, to=3000,by=1)
#time discretization
t=seq(from=1,to=200,by=1)
#matrix to record the profiles along x for each t
matrixResult=matrix(0,nrow=length(t),ncol=length(x))
##########
####Ritter function for water surface profile
###########
yFunction<-function(t,x)
{
Y=(1/(9*9.81))*(2*sqrt(9.81*h1)-(x/t))^(2)
return(Y)
}
h2Function<-function(h2)
{
Y=((sqrt(9.81*h0))*sqrt((h2)/(2*h0)*(1+(h2)/(h0))))-(((-2*sqrt(9.81*h2)+2*sqrt(9.81*h1))*h2)/(h2-h0))
return(Y)
}
#find the root
h2=uniroot(h2Function,interval=c(1,5))
h2=h2$root
#wave celerety
ctc=(((-2*sqrt(9.81*h2)+2*sqrt(9.81*h1))*h2)/(h2-h0))
########
###Computation routine
#######
for(i in 1:length(t))
{
for(j in 1:length(x))
{
matrixResult[i,j]=yFunction(t[i],x[j])
}
}Water Surface Profile Plots
#plotting results for t=200
plot(x[600:which.min(abs(matrixResult[200,] - h2))],matrixResult[200,600:which.min(abs(matrixResult[200,] - h2))],col="darkgreen",lwd=1.5,xlab="Distance [m]",ylab="Elevation [m]",axes=T,xaxs="i",yaxs="i", main="Water Surface Profiles for different time steps",lty=1,type="l",ylim=c(0,5),pch=8,xlim=c(-2000,3000))
#which.min(abs(matrixResult[200,] - h2))
xx=x[which.min(abs(matrixResult[200,] - h2))]
xx=c(xx,(ctc*200),(ctc*200),3500)
yy=matrixResult[200,which.min(abs(matrixResult[200,] - h2))]
yy=c(yy,yy,0.5,0.5)
lines(xx,yy,col="darkgreen",lwd=1.5,lty=1)
#plotting results for t=150
lines(x[948:which.min(abs(matrixResult[150,] - h2))],matrixResult[150,948:which.min(abs(matrixResult[150,] - h2))],type="l",lwd=1.5,col="blue",lty=2,pch=15)
xx=x[which.min(abs(matrixResult[150,] - h2))]
xx=c(xx,(ctc*150),(ctc*150),3500)
yy=matrixResult[150,which.min(abs(matrixResult[150,] - h2))]
yy=c(yy,yy,0.5,0.5)
lines(xx,yy,lwd=1.5,col="blue",lty=2,pch=15)
#plotting results for t=100
lines(x[1297:(which.min(abs(matrixResult[100,1297:3000] - h2))+1297)],matrixResult[100,1297:(which.min(abs(matrixResult[100,1297:3000] - h2))+1297)],type="l",lwd=1.5,col="red",lty=3,pch=16)
xx=x[(which.min(abs(matrixResult[100,1297:3000] - h2))+1297)]
xx=c(xx,(ctc*100),(ctc*100),3500)
yy=matrixResult[100,(which.min(abs(matrixResult[100,1297:3000] - h2))+1297)]
yy=c(yy,yy,0.5,0.5)
lines(xx,yy,lwd=1.5,col="red",lty=3,pch=16)
#plotting results for t=50
lines(x[1651:which.min(abs(matrixResult[50,] - h2))],matrixResult[50,1651:which.min(abs(matrixResult[50,] - h2))],type="l",lwd=1.5,col="purple",lty=4,pch=17)
xx=yy=0
xx=x[which.min(abs(matrixResult[50,] - h2))]
xx=c(xx,(ctc*50),(ctc*50),3500)
yy=matrixResult[50,which.min(abs(matrixResult[50,] - h2))]
yy=c(yy,yy,0.5,0.5)
lines(xx,yy,lwd=1.5,col="purple",lty=4,pch=17)
#x axis editing
axis(1, at = seq(from=-2000,to=3000,by=1000), labels = seq(from=-2000,to=3000,by=1000),lty=1)
#x grid
abline(v=seq(from=-2000,to=3000,by=500), col="lightgray", lty="dotted")
#y axis editing
axis(2, at = seq(from=0,to=5,by=1), labels = seq(from=0,to=5,by=1))
#y grid
abline(h=seq(from=0,to=5,by=0.5), col="lightgray", lty="dotted")
#legend
legend(1300, 4.5, legend=c("t= 200 s","t=150 s", "t=100 s", "t=50 s"),col=c("darkgreen","blue","red","purple"), lty=c(1,2,3,4),lwd=c(1.5,1.5,1.5,1.5),cex=1,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")
box()
Figure 17: Water Surface Profile along the time