CE 573: Sediment Transport
HOME WORK 2: Chapter 3
Professor:Dr. Weiming Wu
By Cássio Rampinelli
February, 14th, 2019
OBS: This script was done in R programming language
__________________________________________________________________________________________________________________________
PROBLEM 3.2. Consider a 2-D uniform turbulent flow in a straight open channel with a channel bed slope of 0.0005 and a flow depth of 1.4 m, compare the production and dissipation rates of turbulence kinetic energy and analyze the kinetic energy budget. (hint: use Eqs. (3.7), (3.35) and (3.36) for the production rate)
Solution
After manipulating equations (3.35) and (3.36), the Production rate \(P_k\) (Eq. 3.7) can be re-written as follows:
\[ P_k=\frac{\tau_{0}}{\rho}\Bigg(1-\frac{y}{h}\Bigg)\cdot \Bigg(\frac{\sqrt{ \frac{ \tau_0}{\rho}}}{\kappa\cdot y}\Bigg) \]
Where
\[ R\simeq h =1.4 m \] Obs:assuming wide channel the hydraulic Radius (R) is approximately the flow depth
\[ \tau_0=\gamma\cdot R \cdot S_0=9810\cdot1.4\cdot0.0005=6.87 N/m² \] \[ \kappa=0.4 \] \[ \rho=1,000 kg/m³ \]
Then, the Production rate \(P_k\) turns to:
\[ P_k=\frac{6.87}{1000}\Bigg(1-\frac{y}{h}\Bigg)\cdot \Bigg(\frac{\sqrt{ \frac{6.87}{1000}}}{0.4\cdot y}\Bigg) \] \[ P_k=0.00687\Bigg(1-\frac{y}{1.4}\Bigg)\cdot \Bigg(\frac{1}{4.8259\cdot y}\Bigg) \]
The Dissipation rate \(\varepsilon\) is given by Eq. 3.20
\[ \frac{\varepsilon\cdot h}{u^3_*}=E_1\Bigg(\frac{y}{h}\Bigg)^{-1/2}\cdot exp\Bigg(-3\frac{y}{h}\Bigg) \]
Where
\[ u_*=\sqrt{ \frac{ \tau_0}{\rho}}=\sqrt{\frac{6.87}{1000}}=0.08289 m/s \] \[ E_1=9.8 \]
Then, the dissipation rate (\(\varepsilon\)) function can be re-arranged as:
\[ \varepsilon=\frac{u^3_*}{h}E_1\Bigg(\frac{y}{h}\Bigg)^{-1/2}\cdot exp\Bigg(-3\frac{y}{h}\Bigg) \] \[ \varepsilon=\frac{5.6952 \cdot 10^{-3} }{1.4}9.8\Bigg(\frac{y}{1.4}\Bigg)^{-1/2}\cdot exp\Bigg(-3\frac{y}{1.4}\Bigg) \]
\[ \varepsilon=3.987\cdot10^{-3}\Bigg(\frac{y}{1.4}\Bigg)^{-1/2}\cdot exp\Bigg(-3\frac{y}{1.4}\Bigg) \]
Finally, we can plot the dissipation rate (\(\varepsilon\)), the production rate (\(P_k\)) and the energy budget rate that is given by the production rate minus the dissipation rate (\(P_k-\varepsilon\)). Everything is plot as a function of the relation y/h, with y varying from 0 until 1.4 m (the water depth). The plot below summarizes that information.
The R code is shown bellow, followed by the plot.
#Clear R memory
rm(list=ls())
#Define the conveyance function and the input parameters
Pk<-function(yh)
{
#Flow depth
h=1.4
Pk=0.00687*(1-yh)*(1/(4.825937*yh*h))
return(Pk)
}
#Define the conveyance function and the input parameters
eps<-function(yh)
{
#Flow depth
eps=0.003986616*((yh)^{-1/2})*exp(-3*yh)
return(eps)
}
#Define the values of h in which Conveyance function will be applied
yh=seq(from=0.01, to=1, by=0.01)
yh2=seq(from=0.01, to=1, by=0.05)
yh3=seq(from=0.02, to=1, by=0.05)
#Applying the conveyance function in all values of the vector h and recorde the result in the vector k
PK=0
E=0
SUM=0
PK2=0
E2=0
SUM2=0
for(i in 1:(length(yh)) )
{
PK[i]=Pk(yh[i])
E[i]=eps(yh[i])
SUM[i]=PK[i]-E[i]
}
for(i in 1:(length(yh2)) )
{
PK2[i]=Pk(yh2[i])
}
for(i in 1:(length(yh3)) )
{
E2[i]=eps(yh3[i])
}#plottin the results
plot(yh,SUM,type='l',col="red",xlab="y/h ",ylab="Rate [J/s]",lwd=2)
lines(yh,E,col="black",lwd=1.5,lty=3,pch=15,cex=3)
points(yh3,E2,col="black",lwd=1.5,lty=3,pch=15,cex=1)
lines(yh,PK,col="blue",lwd=1.5,lty=2,pch=18,cex=3)
points(yh2,PK2,col="blue",lwd=1.5,lty=2,pch=19,cex=1)
#legend(c("Budget Energy Rate","Dissipation #Rate","Production Rate"),lty=c(1,3,2))
legend(0.5, 0.05, c("Budget Enery Rate", "Dissipation Rate", "Production Rate"), col = c("red", "black", "blue"),
text.col = "black", lty = c(1, 3, 2), pch = c(NA, 15, 19),
merge = TRUE)
#title(main="Energy Rate")
#Axis editing
#axis(1, at = seq(from=0,to=1,by=2), labels = c(0,.2,.4,.6,.8,1))
#rug(x = seq(from=0.1,to=0.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")PROBLEM 3.5. Consider a uniform turbulent flow in a wide, straight open channel. If the velocities at 0.2h and 0.8h are measured, can you estimate the unit discharge of this channel? (h is the flow depth. Hint: use the log law or power law).
Solution
Lets assume the Power Law, given by the following expression:
\[ \frac{u}{u_{max}}=\Bigg(\frac{y}{h}\Bigg)^{1/m} \]
Assuming two points, Point 1 (\(u_1\),\(y_1=0.2h\)) and Point 2(\(u_2, y_2=0.8h\)) in the velocity profile (\(u\)) along the the water depth (y), varying from 0 to h, we can solve a system with to points of the former equation to determine u_{max} and m.
\[ \frac{u_1}{u_{max}}=\Bigg(\frac{y_1}{h}\Bigg)^{1/m} \]
\[ \frac{1}{u_{max}}=\frac{1}{u_1}\Bigg(\frac{0.2h}{h}\Bigg)^{1/m} \] \[ \frac{1}{u_{max}}=\frac{0.2^{1/m}}{u_1}\ \ (Eq.1) \]
Repeating the samething for Point 2
\[ \frac{u_2}{u_{max}}=\Bigg(\frac{y_2}{h}\Bigg)^{1/m} \] \[ \frac{1}{u_{max}}=\frac{1}{u_2}\Bigg(\frac{0.8h}{h}\Bigg)^{1/m} \] \[ \frac{1}{u_{max}}=\frac{0.8^{1/m}}{u_2}\ \ (Eq.2) \]
Making (Eq.1)=(Eq.2)
\[ \frac{0.2^{1/m}}{u_1}=\frac{0.8^{1/m}}{u_2} \] \[ \frac{0.2^{1/m}}{0.8^{1/m}}=\frac{u_1}{u_2} \] \[ \Bigg(\frac{0.2}{0.8}\Bigg)^{1/m}=\frac{u_1}{u_2} \] \[ \Bigg(0.25\Bigg)^{1/m}=\frac{u_1}{u_2} \] \[ log\Bigg(0.25\Bigg)^{1/m}=log\Bigg(\frac{u_1}{u_2}\Bigg) \] \[ \frac{1}{m}log\Bigg(0.25\Bigg)=log\Bigg(\frac{u_1}{u_2}\Bigg) \] \[ m=\frac{log(0.25)}{\log(u_1/u_2)} \]
Plugging m in (Eq.2)
\[ \frac{1}{u_{max}}=\frac{0.8^{1/(\frac{log(0.25)}{\log(u_1/u_2)})}}{u_2} \]
\[ u_{max}=\frac{u_2}{0.8^{(\frac{\log(u_1/u_2)}{log(0.25)})}} \]
Finally, the averaged velocity (U) is given by the following equation:
\[ \frac{U}{u_{max}}=\frac{m}{m+1} \]
Plugging \(u_{max}\) and m in the previous equation
The averaged velocity (U) is given by the following equation:
\[ U=u_{max}\frac{m}{m+1} \] \[ U=\Bigg(\frac{u_2}{0.8^{(\frac{\log(u_1/u_2)}{log(0.25)})}}\Bigg)\frac{\frac{log(0.25)}{\log(u_1/u_2)}}{\frac{log(0.25)}{\log(u_1/u_2)}+1} \]
Finally, to estimate the unit discharge we should divide U by the channel width B.
PROBLEM 3.9. Consider the steady, uniform flow in a straight prismatic channel with a compound cross-section, as shown in Fig. 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.035, 0.025 and 0.35, respectively. The longitudinal slope of the channel is 0.001.
Figure 1: Reference Cross Section for Problem 3.5
(1) Find an adequate method to determine the flow conveyance of this channel.
Solution
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):
For: \(0\small\leq h \small\leq 5\)
\[ A(h)= 10h \]
For: \(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:
For: \(0 \small\leq h \small\leq 5\)
\[ R_H(h)=\frac{A(h)}{P_w} \] \[ R_H(h)=\frac{10h}{(10+2h)} \]
For: \(h > 5\)
\[ R_H(h)=\Bigg( \frac{[(h-5)\cdot10]}{(10+(h-5))} \Bigg)_{AL}+\Bigg( \frac{(h\cdot10)}{(10+2\cdot 5)} \Bigg)_{AM}+\Bigg( \frac{[(h-5)\cdot10]}{(10+(h-5))} \Bigg)_{AR} \]
Finally, replacing \(R(h)\) and \(A(h)\) in \(K(h)\)
For: \(0\small\leq h \small\leq 5\)
\[K(h)=\frac{10h }{n_{AM}}\Bigg(\frac{10h}{(10+2h)}\Bigg) R^{(2/3)}_H \]
For: \(h > 5\)
\[ K(h)=\Bigg[ \frac{[(h-5)\cdot10]}{n_{AL}}\Bigg( \frac{[(h-5)\cdot10]}{(10+(h-5))} \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-5))} \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-H)))^(2/3)+((h*bAM)/nAM)*((h*bAM)/(bAM+2*H))^(2/3)+(((h-H)*bAR)/nAR)*(((h-H)*bAR)/(bAR+(h-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.035
#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.35
#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()
(2) If the water depth in the main channel is 2.5 m, what is the flow discharge in this channel?
Solution
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=2.5)\) and substitute the answer in the rating curve.
Since we had defined the \(K(h)\) function, we call the function setting the parameters with the R code as follows:
conveyance(nAL=0.035,nAM=0.025,nAR=0.35,bAL=10,bAM=10,bAR=10,H=5,h=2.5)## [1] 1405.721
Then, \(K(h=2.5)\)=1405.72 m³/s, this brings to:
\[ Q = 1405.72\cdot 0.031623 \] \[ Q = 44.45 m³/s \]
(3) If the flow discharge is 400 m3/s, what is the water depth in this channel?
Solution
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=9.11 m)=400 m³/s.
Then, h=9.11 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 below.
#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.035
#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.35
#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-H)))^(2/3)+((h*bAM)/nAM)*((h*bAM)/(bAM+2*H))^(2/3)+(((h-H)*bAR)/nAR)*(((h-H)*bAR)/(bAR+(h-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] 9.110501PROBLEM 3.10. Consider a straight rectangular channel. The channel bed slope is 0.001, the average flow velocity is 1.5 m/s, and the flow depth is 0.5 m. Determine the bed shear stress when the channel width is 0.5 m and 5.0 m, using the following methods:
(2) Williams (1970) Method.
Solution
Computing the adjust factor of Williams.
\[ Y=\frac{1}{1+0.055h/B^2} \]
For B=0.5 m.
\[ Y=\frac{1}{1+0.055\cdot0.5/B^2} \] \[ Y=\frac{1}{1+0.055\cdot0.5/0.5^2} \]
\[ Y=0.901 \]
For B=5 m.
\[ Y=\frac{1}{1+0.055\cdot0.5/B^2} \]
\[ Y=\frac{1}{1+0.055\cdot0.5/5^2} \]
\[ Y=0.999 \]
The bed shear stress is given by.
\[ \tau_b=\gamma hSY \]
For Y=0.901 (B=0.5 m).
\[ \tau_b=9810 \cdot 0.5 \cdot 0.001 \cdot 0.901 \]
\[ \tau_b=4.419 \ N/m² \]
For Y=0.901 (B=5 m)
\[ \tau_b=9810 \cdot 0.5 \cdot 0.001 \cdot 0.999 \]
\[ \tau_b=4.900 \ N/m² \]
(3) Kinight et al. (1994) Method.
Solution
Knight et al. (1994) derived a relation for the integrated wall shear force, expressed in terms of a percentage of the total shear force carried by the channel (\(SF_w\))
\[ SF_w=C_{ef} \cdot exp\Bigg[4.605-3.23log\Bigg(\frac{P_b}{C_2P_w}+1\Bigg)\Bigg] \]
For B=0.5 m
\[ SF_w=1.0 \cdot exp\Bigg[4.605-3.23\cdot log\Bigg(\frac{0.5}{1.5\cdot 2\cdot0.5}+1\Bigg)\Bigg] \]
\[ SF_w=66.78 \ \% \]
From Eq.3.90
\[ P\tau=P_b\tau_b+P_w\tau_w \]
Dividing both sides of the equation by \(P\tau\)
\[ \frac{P\tau}{P\tau}=\frac{P_b\tau_b}{P\tau}+\frac{P_w\tau_w}{P\tau} \]
Assuming \(\frac{P_w\tau_w}{P\tau}=66.78 \%=0.6678\)
\[ 1=\frac{P_b\tau_b}{P\tau}+0.6678 \]
Assuming \(P=(0.5+2 \cdot 0.5)\) and \(\tau=\gamma \cdot R \cdot \ S=9810 \cdot (0.5 \cdot 0.5)/(0.5+2 \cdot 0.5) \cdot \ 0.001\)
\[ 1=\frac{P_b\tau_b}{(0.5+2\cdot0.5)\cdot \gamma \cdot R\cdot S}+0.6678 \]
\[ 1=\frac{P_b\tau_b}{(0.5+2\cdot0.5)\cdot 9810 \cdot (0.5 \cdot 0.5)/(0.5+2 \cdot 0.5) \cdot \ 0.001}+0.6678 \]
\[ 1=\frac{P_b\tau_b}{2.45}+0.6678 \]
\[ 1-0.6678=\frac{P_b\tau_b}{2.45} \] \[ P_b\tau_b=0.814 \]
\[ 0.5\tau_b=0.814 \]
\[ \tau_b=1.63\ N/m² \]
For B=5 m
\[ SF_w=1.0 \cdot exp\Bigg[4.605-3.23\cdot log\Bigg(\frac{5}{1.5\cdot 2\cdot0.5}+1\Bigg)\Bigg] \]
\[ SF_w=12.78 \ \% \]
From Eq.3.90
\[ P\tau=P_b\tau_b+P_w\tau_w \]
Dividing both sides of the equation by \(P\tau\)
\[ \frac{P\tau}{P\tau}=\frac{P_b\tau_b}{P\tau}+\frac{P_w\tau_w}{P\tau} \]
Assuming \(\frac{P_w\tau_w}{P\tau}=12.78 \%=0.1278\)
\[ 1=\frac{P_b\tau_b}{P\tau}+0.1278 \]
Assuming \(P=5+2 \cdot 0.5\) and \(\tau=\gamma \cdot R \cdot \ S=9810 \cdot (5 \cdot 0.5)/(5+2 \cdot 0.5) \cdot \ 0.001\)
\[ 1=\frac{P_b\tau_b}{(5+2\cdot0.5)\cdot \gamma \cdot R\cdot S}+0.1278 \]
\[ 1=\frac{P_b\tau_b}{(5+2\cdot0.5)\cdot 9810 \cdot (5 \cdot 0.5)/(5+2 \cdot 0.5) \cdot \ 0.001}+0.1278 \]
\[ 1-0.1278=\frac{P_b\tau_b}{21.39} \] \[ 21.39=5\tau_b \] \[ \tau_b=4.28 \ N/m² \]
(4) Guo and Julien (2005) method
Guo and Julien (2005) drived the following formula for the average bed shear stress:
\[ \frac{\tau_b}{\gamma \cdot h \cdot S}=\frac{4}{\pi}arctan\Bigg[exp\Bigg(-\frac{\pi h}{b}\Bigg)\Bigg]+\frac{\pi}{4}\frac{h}{b}exp\Bigg(-\frac{h}{b}\Bigg) \]
For B=0.5 m
\[ \frac{\tau_b}{9810 \cdot 0.5 \cdot 0.001}=\frac{4}{\pi}arctan\Bigg[exp\Bigg(-\frac{\pi 0.5}{0.5}\Bigg)\Bigg]+\frac{\pi}{4}\frac{0.5}{0.5}exp\Bigg(-\frac{0.5}{0.5}\Bigg) \]
\[ \frac{\tau_b}{9810 \cdot 0.5 \cdot 0.001}=0.3439 \] \[ \tau_b=1.67 \ N/m² \]
For B=5 m
\[ \frac{\tau_b}{9810 \cdot 0.5 \cdot 0.001}=\frac{4}{\pi}arctan\Bigg[exp\Bigg(-\frac{\pi 0.5}{5}\Bigg)\Bigg]+\frac{\pi}{4}\frac{0.5}{5}exp\Bigg(-\frac{0.5}{5}\Bigg) \]
\[ \frac{\tau_b}{9810 \cdot 0.5 \cdot 0.001}=0.8743 \] \[ \tau_b=4.29 \ N/m² \]