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 three bed load transport formulas are presented (Meyer-Peter Muller, Einstein Brown and van Rijn). 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 comput Flow Height (Hs), Flow Area (As), Hydraulic Radius (Rs), Wetted Perimeter (Ws) and Flow Velocity(Vs) at cross section x= 2000 ft. The employed code is presented below with comments. Some adjustments were necessary. After some simulations, the upstream boundary condition was set up to flows ranging from 400 cfs to 940, since those were the values that generated flow areas inside the main channel. For the downstream boundary condition the rating curve developed in Homework 1 was adopted. At the end of the code the vectors with each flow condition are presented.
#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)
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)
{
x.left.side=x[1:ceiling(length(x)/2)]
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
y.left.side=y[1:ceiling(length(y)/2)]
y.channel.left.side=y.left.side[(which(x.left.side>=x.left))]
x.right.side=x[(ceiling(length(x)/2)+1):length(y)]
x.channel.right.side=x.right.side[(which(x.right.side<=x.right))]
y.right.side=y[(ceiling(length(x)/2)+1):length(y)]
y.channel.right.side=y.right.side[(which(x.right.side<=x.right))]
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)
perim.channel
A.right=A.left=perim.left=perim.right=0
}
else{
if(y0>y.left && y0>y.right)
{
#Computing left area
x.left.side=x[1:ceiling(length(x)/2)]
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
x.bank.left.side=x.left.side[1:max(which(x.left.side==x.left))]
y.left.side=y[1:ceiling(length(y)/2)]
y.channel.left.side=y.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
y.bank.left.side=y.left.side[1:max(which(x.left.side==x.left))]
x0left= x.bank.left.side[min(which(y.bank.left.side<y0))]-(y0-y.bank.left.side[min(which(y.bank.left.side<y0))])*(x[min(which(y.bank.left.side<y0))]-x.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)])/(y.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)]-y.bank.left.side[min(which(y.bank.left.side<y0))])
xtemp.left=x.bank.left.side[which(x.bank.left.side>x0left)]
xtemp.left=c(x0left,xtemp.left,xtemp.left[length(xtemp.left)])
ytemp.left=y.bank.left.side[which(x.bank.left.side>x0left)]
ytemp.left=c(y0,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
x.right.side=x[(ceiling(length(x)/2)+1):length(x)]
x.channel.right.side=x.right.side[1:max(which(x.right.side==x.right))]
x.bank.right.side=x.right.side[max(which(x.right.side==x.right)):length(x.right.side)]
y.right.side=y[length(y):(ceiling(length(y)/2)+1)]
y.right.side=y.right.side[length(y.right.side):1]
y.channel.right.side=y.right.side[1:max(which(x.right.side==x.right))]
y.bank.right.side=y.right.side[length(y.right.side):max(which(x.right.side==x.right))]
y.bank.right.side=y.bank.right.side[length(y.bank.right.side):1]
x0right= x.bank.right.side[max(which(y.bank.right.side<=y0))]+(y0-y.bank.right.side[max(which(y.bank.right.side<=y0))])*(x.bank.right.side[min(which(y.bank.right.side>y0))]-x.bank.right.side[max(which(y.bank.right.side<=y0))])/(y.bank.right.side[min(which(y.bank.right.side>y0))]-y.bank.right.side[max(which(y.bank.right.side<=y0))])
xtemp.rigth=x.bank.right.side[which(x.bank.right.side<x0right)]
xtemp.rigth=c(xtemp.rigth[1],xtemp.rigth,x0right)
ytemp.right=y.bank.right.side[which(x.bank.right.side<x0right)]
ytemp.right=c(y0,ytemp.right,y0)
#Re-ording the x coordinates
xx=0
xx=xtemp.rigth
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.rigth[2:(length(xtemp.rigth))]
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)
perim.right
#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)
}
else
{
if(y0>y.left && y0<=y.right)
{
#Computing left area
x.left.side=x[1:ceiling(length(x)/2)]
x.channel.left.side=x.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
x.bank.left.side=x.left.side[1:max(which(x.left.side==x.left))]
y.left.side=y[1:ceiling(length(y)/2)]
y.channel.left.side=y.left.side[max(which(x.left.side==x.left)):length(x.left.side)]
y.bank.left.side=y.left.side[1:max(which(x.left.side==x.left))]
x0left= x.bank.left.side[min(which(y.bank.left.side<y0))]-(y0-y.bank.left.side[min(which(y.bank.left.side<y0))])*(x[min(which(y.bank.left.side<y0))]-x.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)])/(y.bank.left.side[match(min(y.bank.left.side[which(y.bank.left.side>y0)]),y.bank.left.side)]-y.bank.left.side[min(which(y.bank.left.side<y0))])
xtemp.left=x.bank.left.side[which(x.bank.left.side>x0left)]
xtemp.left=c(x0left,xtemp.left,xtemp.left[length(xtemp.left)])
ytemp.left=y.bank.left.side[which(x.bank.left.side>x0left)]
ytemp.left=c(y0,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)
A.right=perim.right=0
#Computing channel area
xx=yy=0
xx=c(x.channel.left.side[1],x.channel.left.side,x.channel.right.side)
xx[2:length(xx)]=xx[length(xx):2]
yy=c(y0,y.channel.left.side,y.channel.right.side)
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)
perim.channel
}
else
{
if(y0<y.left && y0>y.right)
{
#Computing right bank area
x.right.side=x[(ceiling(length(x)/2)+1):length(x)]
x.channel.right.side=x.right.side[1:max(which(x.right.side==x.right))]
x.bank.right.side=x.right.side[max(which(x.right.side==x.right)):length(x.right.side)]
y.right.side=y[length(y):(ceiling(length(y)/2)+1)]
y.right.side=y.right.side[length(y.right.side):1]
y.channel.right.side=y.right.side[1:max(which(x.right.side==x.right))]
y.bank.right.side=y.right.side[length(y.right.side):max(which(x.right.side==x.right))]
y.bank.right.side=y.bank.right.side[length(y.bank.right.side):1]
x0right= x.bank.right.side[max(which(y.bank.right.side<=y0))]+(y0-y.bank.right.side[max(which(y.bank.right.side<=y0))])*(x.bank.right.side[min(which(y.bank.right.side>y0))]-x.bank.right.side[max(which(y.bank.right.side<=y0))])/(y.bank.right.side[min(which(y.bank.right.side>y0))]-y.bank.right.side[max(which(y.bank.right.side<=y0))])
xtemp.rigth=x.bank.right.side[which(x.bank.right.side<x0right)]
xtemp.rigth=c(xtemp.rigth[1],xtemp.rigth,x0right)
ytemp.right=y.bank.right.side[which(x.bank.right.side<x0right)]
ytemp.right=c(y0,ytemp.right,y0)
#Re-ording the x coordinates
xx=0
xx=xtemp.rigth
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.rigth[2:(length(xtemp.rigth))]
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
#Computing channel area
xx=yy=0
xx=c(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(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)
perim.channel
}
}
}
}
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
P2=restult[[8]]
Rh2=A2/P2
E2=min(yup)+(ws2-min(yup))+v2^2/(2*32.2)
PL=restult[[5]]
PC=restult[[6]]
PR=restult[[7]]
n=((PL*(nL)^(1.5)+PC*(nL)^(1.5)+PR*(nR)^(1.5))/P2)^(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(delta<0.001){
counter=0
return(ws2.updated)
}
else
{
counter=1
if(ws2.updated>max(yup)){wsu=max(yup)-0.5}
else
wsu=ws2.updated
if(loop>100){
counter=0
return(ws2.updated=wsd+0.5)
}
}
}
}
#####################################
#####################################
##RUNNING THE CODE##################
####################################
#Simulated flows
Qs<-seq(from=400, to=940, by=10)
#Flow heights at Section x=2000 ft
Hs=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
#Boundary Condition
H0=0.302*(Qs^0.358)
#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
}
Hs[j]=WSE[2]
}
#Evaluated Discharges at cross section x=2000 (in cfs)
Qs## [1] 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560
## [18] 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730
## [35] 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900
## [52] 910 920 930 940#Computed Flow depth at cross section x=2000 (in ft)
Hs## [1] 3.079553 3.102457 3.125005 3.147212 3.169089 3.190649 3.211904
## [8] 3.232864 3.253539 3.273940 3.294076 3.313954 3.333584 3.352973
## [15] 3.372129 3.391058 3.409767 3.428263 3.446552 3.464640 3.482532
## [22] 3.500233 3.517749 3.535085 3.552245 3.569233 3.586055 3.602714
## [29] 3.619214 3.635559 3.651752 3.667798 3.683699 3.699459 3.715081
## [36] 3.730568 3.745923 3.761149 3.776248 3.791224 3.806078 3.820814
## [43] 3.835434 3.849939 3.864333 3.878617 3.892793 3.906864 3.920832
## [50] 3.934698 3.948465 3.962133 3.975706 3.989184 4.002569#Function to compute wetted permimeter and wetted area
Area<-function(H){
if((H-2)>4)H=4
b=8
B=b+(H-2)*2/2+(H-2)*5/2
Area=(B+b)*(H-2)/2
return(Area)
}
Permim<-function(H){
if((H-2)>4)H=4
Pe=8+sqrt(((H-2)*2/2)^2+(H-2)^2)+sqrt(((H-2)*5/2)^2+(H-2)^2)
return(Pe)
}
#Computing Flow and Perimeter
for(i in 1:length(Hs))
{
As[i]<-Area(Hs[i])
Ps[i]<-Permim(Hs[i])
}
# Computed Flow Areas at Section x=2000 ft (in ft^2)
As## [1] 10.67593 10.94663 11.21491 11.48086 11.74455 12.00607 12.26547
## [8] 12.52283 12.77820 13.03164 13.28321 13.53297 13.78096 14.02723
## [15] 14.27182 14.51479 14.75616 14.99600 15.23432 15.47117 15.70658
## [22] 15.94059 16.17323 16.40453 16.63452 16.86323 17.09069 17.31692
## [29] 17.54195 17.76581 17.98851 18.21009 18.43056 18.64995 18.86828
## [36] 19.08556 19.30182 19.51707 19.73134 19.94464 20.15699 20.36840
## [43] 20.57890 20.78849 20.99720 21.20503 21.41201 21.61814 21.82345
## [50] 22.02793 22.23162 22.43451 22.63662 22.83796 23.03855#Computed Wetted Permimeters at Section x=2000 (in ft)
Ps## [1] 12.43350 12.52757 12.62017 12.71136 12.80121 12.88975 12.97704
## [8] 13.06312 13.14803 13.23181 13.31451 13.39614 13.47676 13.55639
## [15] 13.63505 13.71279 13.78963 13.86559 13.94070 14.01498 14.08846
## [22] 14.16115 14.23309 14.30428 14.37475 14.44452 14.51360 14.58202
## [29] 14.64978 14.71691 14.78341 14.84930 14.91461 14.97933 15.04349
## [36] 15.10709 15.17015 15.23268 15.29469 15.35619 15.41720 15.47771
## [43] 15.53775 15.59732 15.65643 15.71510 15.77332 15.83110 15.88846
## [50] 15.94541 16.00195 16.05808 16.11382 16.16917 16.22414#Computed Hydraulic Radius at Section x=2000 (in ft)
Rs=As/Ps
Rs## [1] 0.8586423 0.8738031 0.8886497 0.9031964 0.9174566 0.9314430 0.9451671
## [8] 0.9586399 0.9718716 0.9848719 0.9976498 1.0102139 1.0225721 1.0347320
## [15] 1.0467008 1.0584852 1.0700917 1.0815263 1.0927947 1.1039024 1.1148547
## [22] 1.1256564 1.1363122 1.1468266 1.1572040 1.1674483 1.1775634 1.1875531
## [29] 1.1974209 1.2071702 1.2168043 1.2263262 1.2357391 1.2450457 1.2542488
## [36] 1.2633511 1.2723550 1.2812631 1.2900777 1.2988011 1.3074354 1.3159827
## [43] 1.3244451 1.3328244 1.3411227 1.3493417 1.3574832 1.3655489 1.3735403
## [50] 1.3814592 1.3893070 1.3970852 1.4047952 1.4124385 1.42001642.4 Sediment Transport Rate
For \(d_{50}=10 mm\) the Meyer-Peter Muller (1948) formula can be applied to comput the sediment transport rate. Given that this formula requires \(d_{90}\), it was assumed that \(d_90=25 mm\). The computed bed slope (S) was 0.001 for the simulated reach. Since all the flows were inside the main chain the Manning was assumed constant as 0.024. A R function with the Meyer-Peter Muller formula was developed. As input, the hydraulic radius was used. The function returns the sediment load unit discharge (qs) in \((kg/s)/m\). The input hydraulic radius is in (ft) and the function converts it to m. The function and the out put values are presented as follows. The output units are in \([(N/m)/s]\).
qsed<-function(Rh){
gamma=9810
gammas=24525
g=9.81
d=10/1000
n=0.024
d90=25/1000
S=4/4000
k=1/n
kk=26/(d90)^(1/6)
R=Rh*0.3048 #converting hydraulic radius from ft to m
part1=(((k/kk)^(3/2))*gamma*R*S)/(((gammas-gamma)*d))-0.047
#if(part1<0)(part1=0)
part2=part1*((gammas-gamma)*d)/0.25
part3=part1/(gamma/g)^(1/3)
qsediment=(((part3)^2))^(1/3)*gammas/(gammas-gamma)
return(qsediment)
}
qs=0
for(i in 1:length(Rs)){
qs[i]=qsed(Rs[i])
}
qs## [1] 0.03688648 0.03670063 0.03651818 0.03633897 0.03616286 0.03598972
## [7] 0.03581942 0.03565184 0.03548687 0.03532442 0.03516438 0.03500666
## [13] 0.03485118 0.03469785 0.03454661 0.03439737 0.03425006 0.03410463
## [19] 0.03396100 0.03381913 0.03367895 0.03354041 0.03340346 0.03326805
## [25] 0.03313414 0.03300168 0.03287062 0.03274094 0.03261258 0.03248552
## [31] 0.03235971 0.03223513 0.03211174 0.03198951 0.03186840 0.03174840
## [37] 0.03162947 0.03151159 0.03139472 0.03127885 0.03116395 0.03105000
## [43] 0.03093698 0.03082486 0.03071362 0.03060325 0.03049372 0.03038502
## [49] 0.03027713 0.03017003 0.03006370 0.02995812 0.02985329 0.02974918
## [55] 0.02964579DISCUSSION
In the first part this work evaluated 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/m)/s]\). 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. For cross-section at x=2,000 ft, flow elevation ranged between 3.08 ft to 4 ft, considering the discharge range from 400 cfs to 940 cfs. Then, it can be assumed that the bankfull discharge should be close to 940 cfs for this section. For the bankfull discharge the computed unit sediment bed load was 0.03 \([(N/m)/s]\). Considering the width of 15 ft for this section the estimated bed load discharge for this flow should be 0.135 N/s.