CE 470: Stream Riparian Systems and Fluvial Morphology
Name: Cassio Rampinelli
Date: October 02, 2019
Assignment: Homework 3
Note: This script was done in R programming language
______________________________________________________________________________
1.PURPOSE
This assignment is divided in two parts. First the units of four bed load transport formulas are presented (Meyer-Peter Muller, Einstein Brown, van Rijn and Chang). In the second part, three information are determined, as follows:
Flow frequency curve for the flows in the East Branch of the Ausable River;
Flow conditions (height, flow area, wetted perimeter, hydraulic radius and flow velocity) for flows between 10 cfs and the top of the channel in increments of 10 cfs for cross section at x=2,000;and
Sediment transport rate for each of the flows computed in (b) considering a sediment size, d50 = 10 mm.
The data derived from this task is important to determine the effective discharge for the study river reach, which will be focus of a next assignment. The basic data and information used as reference in this assignment can be accessed in Homework 1 and 2 reports.
2.METHODS
2.1 Bed Load - Sediment Transport Equations
2.1.1 Meyer-Peter Muller (1948)
Meyer-Peter Muller equation is given by the following expression:
\[ \Bigg[\frac{q_b(\gamma_s-\gamma)}{\gamma_s}\Bigg]^{2/3}\cdot\Bigg(\frac{\gamma}{g}\Bigg)^{1/3}\cdot \frac{0.25}{(\gamma_s-\gamma)d_m}=\frac{(k/k^{'})^{3/2}\cdot \gamma\cdot R\cdot S}{(\gamma_s-\gamma)d_m}-0.047 \]
In this equation \(q_b\) refers to the unit (or specific) transport rate per unit channel width. The usual units are \([(kg/s)/m]\), \([(N/s/)m]\), or \([m²/s]\), depending on the units adopted for the other variables. In this case, lets assume the following metric units.
\(R= hydraulic \ \ radius \ [m]\)
\(S= energy \ slope \ [dimensionless]\)
\(d_m= mean \ \ particle \ \ diameter \ \ of \ \ the \ \ bed \ \ mixture \ [m]\)
\(\gamma= specific \ \ weight \ \ of \ \ water \ \ [N/m^3]\)
\(\gamma_s= specific \ \ weight \ \ of \ \ sediment \ \ [N/m^3]\)
\(k= 1/n =reciprocal \ \ of \ \ the \ \ Manning \ \ n \ \ of \ \ channel \ \ bed\)
\(k^{'}=26/D_{60}^{1/6}\)
\(\frac{k} {k^{'}} = [dimensionless]\)
\(S = bed \ \ slope \ \ [dimensionless]\)
\(g = gravitational \ \ acceleration \ \ constant \ [m/s^{2}]\)
than,
\(q_b = bed \ \ load \ \ transport \ \ rate \ \ by \ \ weight \ \ per \ \ unit \ \ time \ \ and \ \ width \ \ [(N/m)/s]\)
2.1.2 Einstein Brown (1950)
Einstein Brown equation is given by the following expression:
\[ \Psi=\frac{(\gamma_s-\gamma)d}{\tau_b} \] \[ \Phi= 40\Bigg(\frac{1}{\Psi}\Bigg)^3 \ \ if \ \ \Psi \le 5.5 \]
\[ \Phi= \frac{1}{0.465}e^{-0.391\Psi} \ \ if \ \ \Psi > 5.5 \]
\[ \Phi= \frac{q_b}{\gamma_sF\Bigg[g\Bigg(\frac{\gamma_s-\gamma}{\gamma}\Bigg)d^3\Bigg]^{1/2}} \]
\[ F= \Bigg[\frac{2}{3}+\frac{36\cdot\nu^2}{g\cdot d^3 ((\gamma_s-\gamma)/\gamma)}\Bigg]^{1/2}-\Bigg[\frac{36\cdot\nu^2}{g\cdot d^3 ((\gamma_s-\gamma)/\gamma)}\Bigg]^{1/2} \]
Assuming
\(\tau_b= bed \ \ shear \ \ stress \ \ [N/m³]\)
and the other variables in S.I.
\(q_b= bed \ \ load \ \ discharge \ \ given \ \ as \ \ volume \ \ per \ \ unit \ \ time \ \ [m^3/s]\)
2.1.3 van Rijn (1984)
van Rijn equation is given by the following expression:
\[ q_{bi}=\frac{0.053}{D_{*i}^{0.3}}T_i^{2.1}p_{bi}[(s-1)gd_i^3]^{0.5} \ \ , \ \ if \ \ \tau_b \ge\tau_c \]
\[ q_{bi}=0\ \ , \ \ if \ \ \tau_b <\tau_c \]
\[ D^{*}=d\cdot\Bigg(\Bigg(s-1\Bigg)\cdot \frac{g}{\nu^2}\Bigg)^{1/3} \]
\[ T_i=\Bigg(\frac{\tau_b}{\tau_{c,i}}-1\Bigg) \]
Assuming international units for all variables as done in the previous equations and
\(\nu = kinematic \ \ viscosity \ \ [m^2/s]\)
It yelds. the unit (or specific) transport rate per unit channel width. \(q_b= unit \ \ bed \ \ load \ \ discharge \ \ per \ \ unit \\ channel \ \ given \ \ as \ \ mass \ \ per \ \ \ \ time \ \ \ \ per \ \ unit \ \ width \ \ of \ \ channel \ \ [(kg/s)/m]\)
2.1.4 Chang’s (1980)
Chang’s formula is given by the following expression:
\[ \Phi=6.62 \Bigg(\frac{1}{\Psi}-0.03\Bigg)^{5}\Psi^{3.9} \]
\[ \Phi=\frac{q_s}{[(s-1)gd^3]^{1/2}} \]
\[ \frac{1}{\Phi}=\frac{\tau_0}{g(\rho_s-\rho)d}=\frac{\gamma RS}{g(\rho_s-\rho)d}=\frac{ RS}{(s-1)d} \]
Assuming international units for all variables, It yelds.
\(q_b= unit \ \ bed \ \ load \ \ discharge \ \ per \ \ unit \\ channel \ \ given \ \ as \ \ mass \ \ per \ \ \ \ time \ \ \ \ per \ \ unit \ \ width \ \ of \ \ channel \ \ [(kg/s)/m]\)
2.2 Flow frequency curve
The flow frequency curve is a histogram showing the frequency of discharges among bins for different flow ranges. The code below uses the same data set employed in Homework 1 for East Branch of the Ausable River considering flow bins of 100 cfs. The R code with comments is presented next followed by the flow frequency curve.
#Load packae to handle time series
library(tseries)
#Load txt file and do manipulations/adjustments
dataSet=read.table(file="dataSet.txt",sep="\t")
dataSet=as.matrix(dataSet)
dataSet=dataSet[,-(1:2)]
dataSet=dataSet[,-3]
dataSet=dataSet[-(1:2),]
#Provide names for the two remaining columns
colnames(dataSet)=c("date","discharge_cfs")
#Save only the discharges from the dataSet table in a variable named discharges
discharges<-dataSet[,2]
#Convert the variables as numeric
discharges=as.numeric(discharges)
#Create a time series object to be used next
discharges.ts=ts(discharges,frequency=365.25,start=c(1924,9,5))
#Defining bins with 100 cfs widht
bins=seq(from=10, to=length(discharges.ts),by=100)
#Ploting the frequency with
hist(discharges.ts,breaks=bins,xlim=c(0,3000),xlab="discharges (cfs)",xaxt="none",main="Frequency Curve")
axis(side=1, at=seq(0,3000, 100))
2.3 Flow conditions at different stages at the cross-section x=2,000 ft
The same code presented at Homework 2 was employed to compute Flow Elevation (Ws), Flow Area (As), Hydraulic Radius (Rs), Wetted Perimeter (Ps), Energy Slope (S) and Flow Velocity(Vs) at cross section x= 2000 ft. For more details regarding the cross-section details as well as their stations one can consult Homework 2.
Flow conditions ranged from 10 to 200 cfs by increments of 10 cfs. For downstream boundary condition it was assumed normal depth, which means that for any flow at this section it is assumed that the energy slope is equivalent to bed slope (4/4000=0.001). To facilitate the implementation of the boundary condition for downstream (at Section 1), simulations were performed separated and a rating curve for Section 1 was established. The resultant rating curve for downstream boundary condition is presented below.
The code is presented below with comments. Some adjustments were necessary to set up the original code from Homework 2 to the described running conditions. At the end of the code the flow conditions Ws (Water Surface Elevation) , As (Flow Area), Ps (Wetted Perimeter), Rs (Hydraulic Radius), Vs (Flow Velocity) and Energy Slope (S) at Cross Section 2 (x= 2,000 ft) are presented for the set of flow raging from 10 to 200 cfs.
#Cross Sections Coordinates
x1<-c(0,3,10,18,23,31,38,41)
y1<-c(7,4,4,0,0,4,4,7)
x2<-c(1,5,15,20,28,30,35)
y2<-c(9,4,4,2,2,4,9)
x3<-c(0,4,12,17,25,30,38,38.1)
y3<-c(11,8,8,4,4,9,9,11)
#################
#Function to compute wetted area and wetted perimeter for anny geometry
###################
wettedArea.Perim<-function(x,y,y0,xleft,yleft,xright,yright)
{
##INPUTS
#x= vector with x coordinates
#y= vector with y coordinates
#y0=water surface elevation
#xleft=x left bank
#yleft=y left bank
#xright=x right bank
#yright=y right bank
##OUTPUT: a List with 8 elements
##[[1]]Left Bank Wetted Area
##[[2]]Channel Wetted Area
##[[3]]Right Bank Wetted Area
##[[4]]Total Wetted Area
##[[5]]Left Bank Wetted Perimeter
##[[6]]Channel Wetted Perimeter
##[[7]]Right Bank Wetted Perimeter
##[[8]]Total Wetted Perimeter
#########
x.left=xleft
y.left=yleft
x.right=xright
y.right=yright
if(y0<=y.right && y0<=y.left)
{
###############################
#LEFT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.left.side=x[1:ceiling(length(x)/2)]
#x left side of the main channel
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
#y left side of the bed
y.left.side=y[1:ceiling(length(y)/2)]
#y left side of the main channel
y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
###############################
#RIGHT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
#x left side of the main channel
x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
#y right side of the bed
y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
#y right side of the main channel
y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
#Computing increments
dxLeft=x.left-max(x.channel.left.side)
dyLeft=(y.left-min(y.channel.left.side))
dxRight=(x.right-min(x.channel.right.side))
dyRight=(y.right-min(y.channel.right.side))
x0left=max(x.channel.left.side)-abs(dxLeft/dyLeft)*(y0-min(y.channel.left.side))
x0right=min(x.channel.right.side)+abs(dxRight/dyRight)*(y0-min(y.channel.right.side))
#x0left=x.channel.left.side[length(x.channel.left.side)]-(y0- y.channel.left.side[length(y.channel.left.side)])*(x.channel.left.side[length(x.channel.left.side)]-x.left)/(y.left-y.channel.left.side[length(y.channel.left.side)])
#x0right=x.channel.right.side[min(which(x.channel.right.side<x.right))]+(y0-y.channel.right.side[1])*(x.right-x.channel.right.side[min(which(x.channel.right.side<x.right))])/(y.right-y.channel.right.side[1])
#Computing channel area
xx=yy=0
xx=c(x0left,x.channel.left.side[which(x.channel.left.side>x0left)],x.channel.right.side[which(x.channel.right.side<x0right)],x0right)
xx[2:length(xx)]=xx[length(xx):2]
yy=c(y0,y.channel.left.side[which(x.channel.left.side>x0left)],y.channel.right.side[which(x.channel.right.side<x0right)],y0)
yy[2:length(yy)]=yy[length(yy):2]
a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
a.channel=abs(a.channel)/2
#Computing perimeter
xxp=yyp=0
xxp=c(x0left,x.channel.left.side[which(x.channel.left.side>x0left)],x.channel.right.side[which(x.channel.right.side<x0right)],x0right)
yyp=c(y0,y.channel.left.side[which(x.channel.left.side>x0left)],y.channel.right.side[which(x.channel.right.side<x0right)],y0)
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.channel=sum(dists)
A.right=A.left=perim.left=perim.right=0
}
else{
if(y0>y.left && y0>y.right)
{
###############################
#LEFT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.left.side=x[1:ceiling(length(x)/2)]
#x left side of the main channel
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
#y left side of the bed
y.left.side=y[1:ceiling(length(y)/2)]
#y left side of the main channel
y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
###############################
#RIGHT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
#x left side of the main channel
x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
#y right side of the bed
y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
#y right side of the main channel
y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
######################
#Channel Area
######################
#Computing channel area
xx=yy=0
xx=c(x.channel.left.side[1],x.channel.left.side,x.channel.right.side,x.channel.right.side[length(x.channel.right.side)])
xx[2:length(xx)]=xx[length(xx):2]
yy=c(y0,y.channel.left.side,y.channel.right.side,y0)
yy[2:length(yy)]=yy[length(yy):2]
a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
a.channel=abs(a.channel)/2
#Computing perimeter
xxp=yyp=0
xxp=c(x.channel.left.side,x.channel.right.side)
yyp=c(y.channel.left.side,y.channel.right.side)
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.channel=sum(dists)
#Computing left increments
dxLeft=x.left.side[1]-x.left.side[2]
dyLeft=y.left.side[1]-y.left.side[2]
dxRight=x.right.side[length(x.right.side)]-x.right.side[length(x.right.side)-1]
dyRight=y.right.side[length(y.right.side)]-y.right.side[length(y.right.side)-1]
x0left=x.left.side[1]+abs(dxLeft/dyLeft)*(y.left.side[1]-y0)
x0right=x.right.side[length(x.right.side)]-abs(dxRight/dyRight)*(y.right.side[length(y.right.side)]-y0)
#Computing left area
xtemp.left=x.left.side[(which(x.left.side<=x.left))]
xtemp.left[1]=x0left
xtemp.left=c(xtemp.left,xtemp.left[length(xtemp.left)])
ytemp.left=y.left.side[(which(x.left.side<=x.left))]
ytemp.left[1]=y0
ytemp.left=c(ytemp.left,y0)
#Re-ording the x coordinates
xx=0
xx=xtemp.left
xx[2:length(xx)]=xx[length(xx):2]
#Re-ording the y coordinates
yy=0
yy=ytemp.left
yy[2:length(yy)]=yy[length(yy):2]
#Computing the area below the water surface level
A.left=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
A.left=abs(A.left)/2
#Computing perimeter
xxp=xtemp.left[1:(length(xtemp.left)-1)]
yyp=ytemp.left[1:(length(ytemp.left)-1)]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.left=sum(dists)
#Computing right area
xtemp.right=x.right.side[(which(x.right.side>=x.right))]
xtemp.right[length(xtemp.right)]=x0right
xtemp.right=c(xtemp.right[1],xtemp.right)
ytemp.right=y.right.side[(which(x.right.side>=x.right))]
ytemp.right[length(ytemp.right)]=y0
ytemp.right=c(y0,ytemp.right)
#Re-ording the x coordinates
xx=0
xx=xtemp.right
xx[2:length(xx)]=xx[length(xx):2]
#Re-ording the y coordinates
yy=0
yy=ytemp.right
yy[2:length(yy)]=yy[length(yy):2]
#Computing the area below the water surface level
A.right=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
A.right=abs(A.right)/2
#Computing perimeter
xxp=yyp=0
xxp=xtemp.right[2:(length(xtemp.right))]
yyp=ytemp.right[2:(length(ytemp.right))]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.right=sum(dists)
}
else
{
if(y0>y.left && y0<=y.right)
{
###############################
#LEFT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.left.side=x[1:ceiling(length(x)/2)]
#x left side of the main channel
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
#y left side of the bed
y.left.side=y[1:ceiling(length(y)/2)]
#y left side of the main channel
y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
###############################
#RIGHT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
#x left side of the main channel
x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
#y right side of the bed
y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
#y right side of the main channel
y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
#Computing increments
dxRight=(x.right-min(x.channel.right.side))
dyRight=(y.right-min(y.channel.right.side))
x0right=min(x.channel.right.side)+abs(dxRight/dyRight)*(y0-min(y.channel.right.side))
######################
#Channel Area
######################
#Computing channel area
xtemp=c(x.channel.left.side,x.channel.right.side[which(x.channel.right.side<x0right)])
xtemp=c(xtemp,x0right)
xtemp=c(xtemp[1],xtemp)
ytemp=c(y.channel.left.side,y.channel.right.side[which(x.channel.right.side<x0right)])
ytemp=c(ytemp,y0)
ytemp=c(y0,ytemp)
xx=yy=0
xx=xtemp
xx[2:length(xx)]=xx[length(xx):2]
yy=ytemp
yy[2:length(yy)]=yy[length(yy):2]
a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
a.channel=abs(a.channel)/2
#Computing channel perimeter
xxp=yyp=0
xxp=xtemp[2:length(xtemp)]
yyp=ytemp[2:length(xtemp)]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.channel=sum(dists)
######################
#Left Bank Area
######################
dxLeft=x.left.side[1]-x.left.side[2]
dyLeft=y.left.side[1]-y.left.side[2]
x0left=x.left.side[1]+abs(dxLeft/dyLeft)*(y.left.side[1]-y0)
#Computing left area
xtemp.left=x.left.side[1:max(which(x.left.side==x.left))]
xtemp.left[1]=x0left
xtemp.left=c(xtemp.left,xtemp.left[length(xtemp.left)])
ytemp.left=y.left.side[1:max(which(x.left.side==x.left))]
ytemp.left[1]=y0
ytemp.left=c(ytemp.left,y0)
#Re-ording the x coordinates
xx=0
xx=xtemp.left
xx[2:length(xx)]=xx[length(xx):2]
#Re-ording the y coordinates
yy=0
yy=ytemp.left
yy[2:length(yy)]=yy[length(yy):2]
#Computing the area below the water surface level
A.left=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
A.left=abs(A.left)/2
#Computing left bank perimeter
xxp=xtemp.left[1:(length(xtemp.left)-1)]
yyp=ytemp.left[1:(length(ytemp.left)-1)]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.left=sum(dists)
A.right=perim.right=0
}
else
{
if(y0<=y.left && y0>y.right)
{
###############################
#LEFT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.left.side=x[1:ceiling(length(x)/2)]
#x left side of the main channel
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
#y left side of the bed
y.left.side=y[1:ceiling(length(y)/2)]
#y left side of the main channel
y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
###############################
#RIGHT PART OF THE CROSS SECTION
###############################
#x left side of the bed
x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
#x left side of the main channel
x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
#y right side of the bed
y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
#y right side of the main channel
y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
#Computing increments
dxLeft=(x.left-max(x.left.side))
dyLeft=(y.left-min(y.channel.left.side))
x0left=x.left+abs(dxLeft/dyLeft)*(y0-min(y.channel.left.side))
######################
#Channel Area
######################
#Computing channel area
xtemp=c(x.channel.left.side,x.channel.right.side)
xtemp[1]=x0left
xtemp=c(xtemp,xtemp[length(xtemp)])
ytemp=c(y.channel.left.side,y.channel.right.side)
ytemp[1]=y0
ytemp=c(ytemp,y0)
xx=yy=0
xx=xtemp
xx[2:length(xx)]=xx[length(xx):2]
yy=ytemp
yy[2:length(yy)]=yy[length(yy):2]
a.channel=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
a.channel=abs(a.channel)/2
#Computing channel perimeter
xxp=yyp=0
xxp=xtemp[1:(length(xtemp)-1)]
yyp=ytemp[1:(length(xtemp)-1)]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.channel=sum(dists)
######################
#Right Bank Area
######################
dxRight=x.right.side[length(x.right.side)]-x.right.side[length(x.right.side)-1]
dyRight=y.right.side[length(y.right.side)]-y.right.side[length(y.right.side)-1]
x0right=x.right.side[length(x.right.side)]-abs(dxRight/dyRight)*(y.right.side[length(y.right.side)]-y0)
#Computing left area
xtemp.right=x.right.side[which(x.right.side>=x.right)]
xtemp.right=c(xtemp.right[1],xtemp.right)
ytemp.right=y.right.side[which(y.right.side>=y.right)]
ytemp.right[length(ytemp.right)]=y0
ytemp.right=c(y0,ytemp.right)
#Re-ording the x coordinates
xx=0
xx=xtemp.right
xx[2:length(xx)]=xx[length(xx):2]
#Re-ording the y coordinates
yy=0
yy=ytemp.right
yy[2:length(yy)]=yy[length(yy):2]
#Computing the area below the water surface level
A.right=sum((sum(xx[1:length(xx)-1]*yy[2:length(xx)])+xx[length(xx)]*yy[1])-sum(xx[2:length(xx)]*yy[1:length(xx)-1])-xx[1]*yy[length(xx)])
A.right=abs(A.right)/2
#Computing left bank perimeter
xxp=xtemp.right[2:(length(xtemp.right))]
yyp=ytemp.right[2:(length(ytemp.right))]
xsquared=(xxp[1:length(xxp)-1]-xxp[2:length(xxp)])^2
ysquared=(yyp[1:length(yyp)-1]-yyp[2:length(yyp)])^2
dists=sqrt(xsquared+ysquared)
perim.right=sum(dists)
A.left=perim.left=0
}
}
}
}
a.total=A.left+a.channel+A.right
perim.total=perim.left+perim.channel+perim.right
result=list(A.left,a.channel,A.right,a.total,perim.left,perim.channel,perim.right,perim.total)
return(result)
}
#################
#Function to compute the head loss
###################
deltaEcalc<-function(sfd,sfu,vd,vu,dx){
Sf1=sfd
Sf2=sfu
deltax=dx
v1=vd
v2=vu
Sfbar=mean(c(Sf1,Sf2))
veld=vd
velu=vu
if(veld>velu){C=0.1}
else C=0.3
#C=0.3
deltaEcalc=Sfbar*deltax+C*abs((velu^2)/(2*32.2)-(veld^2)/(2*32.2))
return(deltaEcalc)
}
###############################
#Standard step method function
###############################
findWsd<-function(ws1,v1,E1,Sf1,xup,yup,yd,posD,posU,nL,nC,nR,xleft,yleft,xright,yright,Q){
##INPUTS
#ws1= water surface elevation from downstream section
#v1= velocity from downstream section
#E1=Energy from downstream section
#Sf1= Slope of the Energy grade line at downstream section
#xup=x coordinates from the upstream section
#yup=y coordinates from the upstream section
#yd= y coordinates from the downstream section
#posD=Station of the downstream section
#posU=Station of the upstream section
#nL=Left bank Manning
#nC=Channel Manning
#nR-Right bank Manning
#xleft=x left bank
#yleft=y left bank
#xright=x right bank
#yright=y right bank
##OUTPUT: Water surface elevation from the downstream section
ws1=ws1
v1=v1
E1=E1
Sf1=Sf1
xup=xup
yup=yup
yd=yd
posD=posD
posU=posU
nL=nL
nC=nC
nR=nR
xleft=xleft
yleft=yleft
xright=xright
yright=yright
Q=Q
counter=1
wsu=(ws1-min(yd))+min(yup)
loop=0
while(counter==1){
#Section 2
#y02=min(y2)+(y0-min(y1))
ws2=wsu
restult=wettedArea.Perim(x=xup,y=yup,y0=ws2,xleft=xleft,yleft=yleft,xright=xright,yright=yright)
A2=restult[[4]]
v2=Q/A2
PT=restult[[8]]
Rh2=A2/PT
E2=ws2+v2^2/(2*32.2)
PL=restult[[5]]
PC=restult[[6]]
PR=restult[[7]]
n=((PL*(nL)^(1.5)+PC*(nC)^(1.5)+PR*(nR)^(1.5))/PT)^(2/3)
Sf2=((n*Q)/(1.49*A2*(Rh2^(2/3))))^2
#Checking
deltaE=E2-E1
deltaE.calc=deltaEcalc(Sf1,Sf2,v1,v2,(stations[posU]-stations[posD]))
ws2.calc=ws1+(v1^2)/(2*32.2)+deltaE.calc-((v2^2)/(2*32.2))
ws2.updated=(ws2.calc+ws2)/2
delta=abs(ws2.updated-ws2)
loop=loop+1
if(loop>=200){
counter=0
return(ws2.updated)
}
if(delta<0.001){
counter=0
return(ws2.updated)
}
else
counter=1
wsu=ws2.updated
}
}
#####################################
#####################################
##RUNNING THE CODE##################
####################################
#Simulated flows
Qs<-seq(from=10, to=200, by=10)
#Flow heights at Section x=2000 ft
Ws=0
#Flow Areas at Section x=2000 ft
As=0
#Wetted Permimeters at Section x=2000 ft
Ps=0
#Hydraulic Radius at Section x=2000 ft
Rs=0
#Velocity at Section x=2000 ft
Vs=0
#Energy Slope at Section x=2000 ft
S=0
#Boundary Condition
H0=0.3115*(Qs^0.499)
#Saving x an y matrix with the coordinates
xMatrix=qpcR:::cbind.na(x1,x2,x3)
yMatrix=qpcR:::cbind.na(y1,y2,y3)
#Stations of the sections
stations=c(0,2000,4000)
#Left Banks
xsleft<-c(x1[3],x2[3],x3[3])
ysleft<-c(y1[3],y2[3],y3[3])
#Right Banks
xsright<-c(x1[6],x2[6],x3[6])
ysright<-c(y1[6],y2[6],y3[6])
#Manning left bank
Nleft<-c(0.04,0.045,0.04)
#Manning channel
Nchannel<-c(0.025,0.024,0.025)
#Manning right bank
Nright<-c(0.035,0.024,0.035)
for (j in 1:length(Qs)){
#Water surface elevation at downstream section(ft)
wsd=ws1=H0[j]
#Vector with the water surface elevation of the sections
WSE=ws1
Q=Qs[j]
#Initiating Variables
A=0
v=0
P=0
Rh=0
E=0
Sf=0
nsec=3
#downstream x coordinates without NA
x1<-xMatrix[,1]
x1<-x1[!is.na(x1)]
#downstream y coordinates without NA
y1<-yMatrix[,1]
y1<-y1[!is.na(y1)]
#Section downstream
restult=wettedArea.Perim(x=x1,y=y1,y0=ws1,xleft=xsleft[1],yleft=ysleft[1],xright=xsright[1],yright=ysright[1])
A[1]=restult[[4]]
v[1]=Q/A[1]
P[1]=restult[[8]]
Rh[1]=A[1]/P[1]
E[1]=wsd+v[1]^2/(2*32.2)
PL=restult[[5]]
PC=restult[[6]]
PR=restult[[7]]
n=((PL*(Nleft[1])^(1.5)+PC*(Nchannel[1])^(1.5)+PR*(Nright[1])^(1.5))/P[1])^(2/3)
Sf[1]=((n*Q)/(1.49*A[1]*(Rh[1]^(2/3))))^2
######################################
#Loop to apply the standard step method along the cross sections
######################################
for(i in 2:nsec){
xd=yup=0
#downstream x coordinates without NA
xd<-xMatrix[,i-1]
xd<-xd[!is.na(xd)]
#downstream y coordinates without NA
yd<-yMatrix[,i-1]
yd<-yd[!is.na(yd)]
#upwnstream x coordinates without NA
xup<-xMatrix[,i]
xup<-xup[!is.na(xup)]
#upwnstream y coordinates without NA
yup<-yMatrix[,i]
yup<-yup[!is.na(yup)]
wsup=findWsd(ws1=WSE[i-1],v1=v[i-1],E1=E[i-1],Sf1=Sf[i-1],xup=xup,yup=yup,yd=yd,posD=(i-1),posU=i,nL=Nleft[i],nC=Nchannel[i],nR=Nright[i],xleft=xsleft[i],yleft=ysleft[i],xright=xsright[i],yright=ysright[i],Q)
restult=wettedArea.Perim(x=xup,y=yup,y0=wsup,xleft=xsleft[i],yleft=ysleft[i],xright=xsright[i],yright=ysright[i])
A[i]=restult[[4]]
v[i]=Q/A[i]
P[i]=restult[[8]]
Rh[i]=A[i]/P[i]
E[i]=wsd+v[i]^2/(2*32.2)
PL=restult[[5]]
PC=restult[[6]]
PR=restult[[7]]
n=((PL*(Nleft[i])^(1.5)+PC*(Nchannel[i])^(1.5)+PR*(Nright[i])^(1.5))/P[i])^(2/3)
Sf[i]=((n*Q)/(1.49*A[i]*(Rh[i]^(2/3))))^2
WSE[i]=wsup
}
Ws[j]=WSE[2]
As[j]=A[2]
Ps[j]=P[2]
Rs[j]=Rh[2]
Vs[j]=v[2]
S[j]=Sf[2]
}
#Simulated Flow Discharges (cfs)
Qs## [1] 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170
## [18] 180 190 200#Water Surface Elevations at Cross Section 2 (ft)
Ws## [1] 3.105144 3.195480 3.518984 3.788676 4.332319 4.501835 4.660426
## [8] 4.810400 4.953107 5.089204 5.219933 5.345726 5.467101 5.584379
## [15] 5.698064 5.808444 6.986889 6.877436 6.821552 6.800607#Flow Area at Cross Section 2 (ft^2)
As## [1] 10.97851 12.06489 16.18967 19.90830 31.40736 35.77253 39.90320
## [8] 43.85108 47.64526 51.29782 54.83775 58.27304 61.61467 64.86872
## [15] 68.04667 71.15454 105.70157 102.38757 100.70384 100.07425#Wetted Perimeter at Cross Section 2 (ft)
Ps## [1] 12.53860 12.90959 14.23816 15.34573 27.10914 27.56596 27.99333
## [8] 28.39749 28.78206 29.14882 29.50111 29.84011 30.16719 30.48324
## [15] 30.78960 31.08706 34.26277 33.96782 33.81722 33.76078#Hydraulic Radius at Cross Section 2 (ft)
Rs## [1] 0.8755767 0.9345677 1.1370621 1.2973184 1.1585524 1.2977070 1.4254536
## [8] 1.5441886 1.6553802 1.7598593 1.8588364 1.9528429 2.0424397 2.1280127
## [15] 2.2100539 2.2888797 3.0850266 3.0142523 2.9778866 2.9642164#Flow Velocity at Cross Section 2 (ft/s)
Vs## [1] 0.9108706 1.6577031 1.8530336 2.0092127 1.5919838 1.6772645 1.7542454
## [8] 1.8243564 1.8889604 1.9494005 2.0059176 2.0592713 2.1098871 2.1582051
## [15] 2.2043693 2.2486269 1.6083016 1.7580259 1.8867205 1.9985161#Energy Slope at Cross Section 2 (-)
S## [1] 0.0002569830 0.0007802776 0.0007506503 0.0007402408 0.0010136338
## [6] 0.0009691156 0.0009370301 0.0009123627 0.0008928579 0.0008776137
## [11] 0.0008649866 0.0008546216 0.0008460392 0.0008390052 0.0008330914
## [16] 0.0008281194 0.0002872319 0.0003537049 0.0004138615 0.00046714522.4 Sediment Transport Rate
For \(d_{50mm}=10 mm\) the Meyer-Peter Muller (1948) formula can be applied to comput the sediment transport rate given the formula was developed for sediments between 6.4 and 30 mm.
Since this formula requires \(d_{90mm}\), it was assumed uniform distribution for the sediments which means \(d_{90mm}=d_{50mm} =10mm\). A R function with the Meyer-Peter Muller formula was developed. As input, the hydraulic radius and the energy slope computed in the previous section was used.The hydraulic radius input unit is in (ft) and the function converts it to metric units. The function returns the sediment unit discharge (qs) per widht of the channe, in \((N/s)/m\). The input hydraulic radius is in (ft) and the function converts it to m. The function and the computed output are presented as follows. The output unit for the sediment transport rate is in \([(N/s)/m]\).
The R code is presented below followed by the sediment transport rates (qb) for the simulated flow conditions.
Note that for extremely low flow depth the formula can not compute sediment discharge appropriately, since the fraction before the constant 0.047 becomes less than it. In this case, it was assumed 0 for the sediment transport rate.
##############################
###Sediment Transport Rate
##############################
qsed<-function(Rh,S){
Rh=Rh
S=S
gamma=9810
gammas=24525
g=9.81
d=10/1000
n=0.024
d90=25/1000
k=1/n
kk=26/(d90)^(1/6)
R=Rh*0.3048 #converting hydraulic radius from ft to m
part1=gammas*sqrt((gammas/gamma-1)*g*d^3)
part2=((k/kk)^(3/2))*gamma*Rh*S
part3=(gammas-gamma)*d
if(part2/part3<0.047){
qsed=0
return(qsed)}
else{
qsed=part1*8*(part2/part3-0.047)^(3/2)
return(qsed)
}
}
qb=0
for(j in 1:length(Rs)){
qb[j]=qsed(Rs[j],S[j])
}
#Returns unit sediment transpor rate (N/s/m)
qb## [1] 0.00000000 0.00000000 0.00000000 0.23833814 1.54551837 2.23086867
## [7] 2.94523945 3.67231061 4.40703659 5.15394028 5.90182288 6.65361245
## [13] 7.40842397 8.16810815 8.92903479 9.69123543 0.01261662 0.79102868
## [19] 2.01477643 3.42654349DISCUSSION
The first part of this assignment evaluated the the units used for 4 sediment transport formulas. It was observed that the sediment transport rate (\(q_b\)) depends on the units of the parameters considered for each formula. In general \(q_b\) refers to the unit (or specific) transport rate per unit channel width. The usual units are \([(kg/s)/m]\), \([(N/s/)m]\), or \([m²/s]\). Assuming metric units, The Meyer-Peter Muller Formula returns the unit sediment transport rate in \([(N/s)/m]\). The Einstein Brown formula returns the sediment transport rate in \([m^3/s]\), while van Rijn and Cheng formulas return it in \([(kg/s)/m]\)
The flow frequency curve shows that the most frequency discharges are between 0 and 200 cfs, which corresponds to most of the simulated discharges (10 to 200 cfs).
For cross-section at x=2,000 ft, considering the simulated discharges, water surface elevation ranged between 3.11 ft to 6.8 ft. Considering the invert elevation for this cross section (invert =2ft), it corresponds to flow depths ranging from 1.1 to 4.8 ft. Considering the water surface elevations outputs and Cross Section 2, it can be verified that the bankfull discharge is between 40 and 50 cfs. The figure below shows the flow elevations for this reference discharges and the bank elevation at 4 ft.
plot(x2,y2,type="l",main="Cross-Sections", xlab="Station(ft)",ylab="Elevation(f)",xlim=c(0,40),ylim=c(0,12),lwd=3)
abline(h=Ws[4],col="darkblue")
abline(h=Ws[5],col="darkred",lty=2)
grid()
legend("top",lty=c(1,1,2),lwd=c(3,1,1),bty="n",cex=1.5,col=c("black","darkblue","darkred"),
legend=c("Section","WSE=3.79 ft and Q=40 cfs","WSE=4.33 ft and Q=50 cfs"))
Finally, with regards the sediment transport rate it could be noted that for flow discharges until 30 cfs it can be negligible. If we considered the sediment transport only until the bankfull discharge, considering that it could be assumed 50 cfs, the computed unit sediment transport rate at this discharge is 1.55 \((N/s)/m\). Adopting an average width of 11.5 ft (3.5 m) for the channel, it corresponds to a sediment rate of 5.43 \(N/s\) or 0.554 \(kg/s\) which is significant and must be taken in account in a channel restoration project.