CE 573: Sediment Transport
HOME WORK 4: Chapter 5
Professor:Dr. Weiming Wu
By Cássio Rampinelli
March, 15th, 2019
OBS: This script was done in R programming language
__________________________________________________________________________________________________________________________
PROBLEM 5.3. Determine the sediment size at incipient motion when the shear velocity \(u_*\)=0.05 m/s. Use one of the modified Shields diagrams.
Solution
First we should compute the shear stress (\(\tau\)) employing the given shear velocity
\[ u_{*}=\sqrt{\frac{\tau}{\rho}} \] \[ \tau=u_*^2\rho \] \[ \tau=0.05^2 \cdot 1000 \] \[ \tau=2.5 \ N/m² \]
Second, we should compute the shields number (\(\Theta\))by:
\[ \Theta=\frac{\tau}{(\gamma_s-\gamma)d} \] \[ \Theta=\frac{2.5}{(25996.5-9810)d} \] \[ \Theta=\frac{2.5}{(16186.5)d} \]
Then, we can make the previous relation equals to some Shields Diagram approxiamtion expression, for instance, the one proposed by Paphitis (2001) that is given by:
\[ \Theta_c=\frac{0.273}{1+1.2D_*}+0.046(1-0.576e^{-0.02D_*}) \]
In which:
\[ D_*=d\Bigg[\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\frac{g}{\nu^2}\Bigg]^{1/3} \] \(\rho_s=2650 kg/m³\)
\(\rho=1000 kg/m³\)
\(g=9.81m/s²\)
\(\nu=10^{-6}\)
Setting the expressions equal to eachothers
\[ \frac{0.273}{1+1.2d\Bigg[\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\frac{g}{\nu^2}\Bigg]^{1/3}}+0.046\Bigg(1-0.576exp\Bigg[-0.02d\Bigg[\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\frac{g}{\nu^2}\Bigg]\Bigg)\Bigg]^{1/3}\Bigg)=\frac{2.5}{(16186.5)d} \]
Solving the expression above for d, using a solving function, brings to:
\[ d=0.00349\ m=3.49\ mm \]
PROBLEM 5.5. Consider a nonuniform sediment mixture consisting of three size classes with representative diameter 0.05, 0.5, and 5 mm and fraction of 0.3, 0.4 and 0.3, respectively. Determine the critical shear stress using the formula of Wu et al. (2000).
Solution
First we should compute \(p_{hk}\) and \(p_{ek}\) for each of the 3 classes regarding the three representative diameters given by:
\[ p_{hk}=\sum_{j=1}^Np_{bj}\frac{d_j}{d_k+d_j} \] \[ p_{ek}=\sum_{j=1}^Np_{bj}\frac{d_k}{d_k+d_j} \]
For k=1 (d=0.05 mm)
\[ p_{h1}=0.3\frac{0.05}{0.05+0.05}+0.4\frac{0.5}{0.05+0.5}+0.3\frac{5}{0.05+5} \] \[ p_{h1}=0.81 \]
\[ p_{e1}=0.3\frac{0.05}{0.05+0.05}+0.4\frac{0.05}{0.05+0.5}+0.3\frac{0.05}{0.05+5} \]
\[ p_{e1}=0.19 \]
For k=2 (d=0.5 mm)
\[ p_{h2}=0.3\frac{0.05}{0.5+0.05}+0.4\frac{0.5}{0.5+0.5}+0.3\frac{5}{0.5+5} \] \[ p_{h2}=0.5 \]
\[ p_{e2}=0.3\frac{0.5}{0.5+0.05}+0.4\frac{0.5}{0.5+0.5}+0.3\frac{0.5}{0.5+5} \]
\[ p_{e2}=0.5 \]
For k=3 (d=5 mm)
\[ p_{h3}=0.3\frac{0.05}{5+0.05}+0.4\frac{0.5}{5+0.5}+0.3\frac{5}{5+5} \] \[ p_{h3}=0.19 \]
\[ p_{e3}=0.3\frac{5}{5+0.05}+0.4\frac{5}{5+0.5}+0.3\frac{5}{5+5} \]
\[ p_{e3}=0.91 \]
Then, the one can calculate the critical shear stress for each representative diameter (or class k) using the Wu et al. (2000) formula proposed for the Shields sediment incipient motion:
\[ \frac{\tau_{ck}}{(\gamma_s-\gamma)d_k}=\Theta_c\Bigg(\frac{p_{ek}}{p_{hk}}\Bigg)^{-m} \]
Assuming \(\Theta_c=0.03\) and m=0.6, obteined from calibrated data, and \(\gamma_s=25,996.5 N/m³\), \(\gamma=9810 N/m³\) the critical shear stress fro each representative diamter can be calculated as follows:
For k=1 (d=0.05 mm)
\[ \frac{\tau_{c1}}{(25996.5-9810)0.00005}=0.03\Bigg(\frac{0.19}{0.81}\Bigg)^{-0.6} \]
\[ \tau_{c1}=0.058\ N/m² \]
For k=2 (d=0.5 mm)
\[ \frac{\tau_{c2}}{(25996.5-9810)0.0005}=0.03\Bigg(\frac{0.5}{0.5}\Bigg)^{-0.6} \]
\[ \tau_{c2}=0.24\ N/m² \]
For k=3 (d=5 mm)
\[ \frac{\tau_{c3}}{(25996.5-9810)0.005}=0.03\Bigg(\frac{0.81}{0.19}\Bigg)^{-0.6} \]
\[ \tau_{c2}=1.02\ N/m² \]
PROBLEM 5.6. Compare the original and revised Shields diagrams and the critical shear velocity Eq. (5.33) by graphically showing these relations in terms of \(\theta_c\) and \(D^*\)
We will use the following approximations for the revised shields diagram:
Soulsby (1997)
\[ \Theta_c=\frac{0.24}{D_*}+0.055(1-e^{-D_*/50}) \]
Paphitis (2001)
\[ \Theta_c=\frac{0.273}{1+1.2D_*}+0.046(1-0.576e^{-0.02D_*}) \]
Simões(2014) proposed the following equation (Eq. 5.33)
\[ A_c=0.215+\frac{6.79}{D_*^{1.7}}-0.075e^{-0.00262D_*} \]
Where
\[ A_c=\frac{u_{*c}}{\omega_s} \] \[ u_{*c}=A_c \cdot \omega_s \]
We will use the Cheng (1997a) expression for the fall velocity (\(\omega_s\))given by:
\[ \omega_s=\frac{\nu}{d}\Bigg(\sqrt{(25+1.2D_*^2)}-5\Bigg)^{1.5} \]
The relation between \(D_*\) and d is given by:
\[ d=\frac{D_*}{\Bigg[\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\frac{g}{\nu^2}\Bigg]^{1/3}} \]
The relation between \(\tau\) and the shear stress velocity is given by:
\[ u_*=\sqrt{\frac{\tau}{\rho}} \] \[ \tau=u_*^2\cdot\rho \]
Since Simões(2014) expression doesn’t provide the critical Shields number directly, one can make the following arrangement based on the previous relations.
\[ \Theta=\frac{\tau}{(\gamma_s-\gamma)d} \] \[ \Theta=\frac{u_*^2\cdot\rho}{(\gamma_s-\gamma)d} \] \[ \Theta=\frac{(A_c \cdot \omega_s)^2\cdot\rho}{(\gamma_s-\gamma)d} \]
The expression above is the Simões(2014) expression adjusted for the Shields number. This expression can be now compared with the Soulsby(1997), Wu and Wang(1999) and Paphitis (2001). The following plot presents the comparison between the mentioned equations. It can be seen that the Simões (2014) equation usually underestimates \(\Theta\) in relation to the Soulsby (1997) and Paphitis (2011) expresions.The R code is also presented.
paphitis<-function(D)
{
theta=0.273/(1+1.2*D)+0.046*(1-0.576*exp(-0.02*D))
return(theta)
}
soulsby<-function(D)
{
theta=0.24/D+0.055*(1-exp(-D/50))
return(theta)
}
simoes<-function(D)
{
A=0.215+6.79/(D^1.7)-0.075*exp(-0.00262*D)
omega=(1/1000000)/((D)/((2.65-1)*9.81/(1/1000000)^2)^(1/3))*(sqrt(25+1.2*D^2)-5)^1.5
ustar=(A*omega)
theta=ustar^2*1000/((25996.5-9810)*((D)/((2.65-1)*9.81/(1/1000000)^2)^(1/3)))
return(theta)
}
D=seq(from=0.01,to=10000,by=0.1)
ypaphitis=paphitis(D)
ysoulsby=soulsby(D)
ysimoes=simoes(D)
source("https://raw.githubusercontent.com/petrkeil/Blog/master/2016_07_05_Log_scales/loglogplot.r")
loglog.plot(xlab="D*", ylab=expression(paste(Theta)), ylim=c(0.0001, 0.5),xlim=c(1,10000))
lines(log(D),log(ysoulsby),lty=1,col="red",lwd=1.5)
lines(log(D),log(ypaphitis),lty=4,col="blue",lwd=1.5)
points(log(D),log(ysimoes),pch=15)
legend("top",c("Soulsby (1997)","Paphitis (2011)","Simões (2014)"),pch=c(NA,NA,15),col=c("red","blue","black"),lty=c(1,4,NA),lwd=c(1.5,1.5,NA))
PROBLEM 5.7. Determine the stability of the riprap of 10 cm in diameter placed on a bank with a 1V:2H side slope if the shear stress is 20 Pa.
First we compute \(D^*\)
\[ D_*=d\Bigg[\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\frac{g}{\nu^2}\Bigg]^{1/3} \] \[ D_*=2529.6 \]
For the calculated \(D^*\) we should use the last range from Eq.5.31, which brings to:
\[ \frac{\tau_{c}}{(\gamma_s-\gamma)d}=0.052 \]
\[ \frac{\tau_{c}}{(25996.5-9810)0.1}=0.052 \] \[ \tau_{c}=84.17\ N/m² \]
The critical shear stress limit (84.17 N/m²) is less than the informed value of 20 N/m². However we also should check the banks, as follows:
From equation 2.35, one can compute the repose angle (\(\phi_r\)) (where d is in mm):
\[ \phi_r=32.5+1.27d \] \[ \phi_r=32.5+1.27\cdot 10 \] \[ \phi_r=45.2 \ degrees \]
The bank angle (\(\beta\)) can be computed by the following equation, considering the given bank slope (where d is in mm):
\[ \beta=atan(1/2) \] \[ \beta=26.6\ degrees \]
The critical shear stress limit for the banks (\(\tau_{c\phi}\)) is given by:
\[ \tau_{c\phi}=\tau_c\sqrt{1-\frac{sin^2{(\beta)}}{\sin^2(\phi_r)}} \] \[ \tau_{c\phi}=84.17\sqrt{1-\frac{sin^2{(26.6)}}{sin^2(45.2)}} \] \[ \tau_{c\phi}=65.30 \ N/m² \]