CE 573: Sediment Transport
HOME WORK 3: Chapter 4
Professor:Dr. Weiming Wu
By Cássio Rampinelli
February, 25th, 2019
OBS: This script was done in R programming language
__________________________________________________________________________________________________________________________
PROBLEM 4.2. Describe the relationship of the drag coefficient and sediment size. How do the particle shape and angularity affect the drag coefficient?
Solution
Several formulas have been developed to consider particle size as the representative property of particle geometry. Others have added particle shape and roundness. If we consider a spherical particle and only the effect of the size of the particle (represented by its diameter), for instance, biger particles tend to have higher settling velocities, which means lower drag coefficient. McNown et al. (1951) studied the effect of particle shape on the settling velocity for particles with different shapes. They introduced a correction factor K for the drag force that is equal to the ratio of the settling velocity of a sphere with the same volume and weight as the particle to the settling velocity of the particle. For Reynolds numbers less than 0.1, the results demonstrated that non-spherical particles experience higher drag and settle slower than the equivalent spheres in most of the cases, except for ellipsoid with the primary (settling) axis nearly twice the other two axes (\(l_1/\sqrt(l_2l_3)\) around 2).
It is observed that by virtue of a decrease in sphericity (low shape factors) and increase in angularity the settling velocities are usually lower, which means that the drag coefficient is higher. In addition, flow separation which increases drag is more likely to occur for non-spherical particles which may ultimately influence the settling path. Thus, natural particles tend to have a lower settling velocity compared to perfectly spherical particles.
Figure 1: Drag coefficient as function of Reynolds number and particle shape (Wu and Wang, 2006)
Figure 2: Comparison of Wu and Wang (2006) formula with the method of U.S. Interagency Committee
PROBLEM 4.9. Calculate the settling velocities of natural sediment particles with size of 0.01 mm and 1 mm at water temperature of 20C, using the formulas of Rubey (1933), Zhang (1961), Wu and Wang (2006), Ahrens (2000), and Jimenez and Madsen (2003).
Solution
* Rubey (1933), for d= 1 mm;
and assuming, g=9.81 \(m/s^2\), F=0.79, \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ w_s=F\sqrt{\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\cdot gd} \] \[ w_s=0.79\sqrt{\Bigg(\frac{2650}{1000}-1\Bigg)\cdot 9.81\cdot 0.001} \] \[ w_s=0.10\ m/s \]
* Rubey (1933), for d= 0.01 mm;
and assuming, g=9.81 \(m/s^2\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\), \(\nu=10^{-6}\ m^2/s\) F is computed by:
\[ F=\Bigg[\frac{2}{3}+\frac{36\nu^2}{gd^3(\rho_s/\rho-1)}\Bigg]^{1/2}-\Bigg[\frac{36\nu^2}{gd^3(\rho_s/\rho-1)}\Bigg]^{1/2} \]
\[ F=\Bigg[\frac{2}{3}+\frac{36\cdot 10^{-6^{2}}}{9.81 \cdot (0.00001)^3(2650/1000-1)}\Bigg]^{1/2}-\Bigg[\frac{36\cdot 10^{-6^{2}}}{9.81 \cdot (0.00001)^3(2650/1000-1)}\Bigg]^{1/2} \] \[ F=0.00707 \]
Then,
\[ w_s=0.00707\sqrt{\Bigg(\frac{2650}{1000}-1\Bigg)\cdot 9.81\cdot 0.00001} \]
\[ w_s=8.99\cdot10^{-5} \ m/s \]
* Zhang (1961), for d= 1 mm;
and assuming, g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ w_s=\sqrt{\Bigg(13.95\cdot\frac{\nu}{d}\Bigg)^2+1.09\Bigg(\frac{\rho_s}{\rho}-1\Bigg)gd}-13.95\frac{\nu}{d} \]
\[ w_s=\sqrt{\Bigg(13.95\cdot\frac{10^{-6}}{0.001}\Bigg)^2+1.09\Bigg(\frac{2650}{1000}-1\Bigg)9.81 \cdot (0.001)}-13.95\frac{10^{-6}}{0.001} \]
\[ w_s=0.12 \ m/s \]
* Zhang (1961), for d= 0.01 mm;
and assuming, g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ w_s=\sqrt{\Bigg(13.95\cdot\frac{10^{-6}}{0.00001}\Bigg)^2+1.09\Bigg(\frac{2650}{1000}-1\Bigg)9.81 \cdot (0.00001)}-13.95\frac{10^{-6}}{0.00001} \] \[ w_s=6.32\cdot10^{-5}\ m/s \]
* Wu and Wang (2006), for d= 1 mm;
and assuming,\(S_p=0.7\), g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ M=53.5e^{-0.65S_p};N=5.65e^{-2.5S_p}; n=0.7+0.9S_p \] \[ M=53.5e^{-0.65\cdot 0.7};N=5.65e^{-2.5 \cdot 0.7}; n=0.7+0.9\cdot 0.7 \] \[ M=33.94;N=0.98; n=1.33 \]
\[ D_*^3=\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\cdot g\cdot d^3/\nu^2 \]
\[ D_*^3=\Bigg(\frac{2650}{1000}-1\Bigg)\cdot 9.81\cdot (0.001)^3/10^{-6^{2}} \] \[ D_*^3=16186.5 \]
\[ w_s=\frac{M\nu}{Nd}\Bigg[\sqrt{\frac{1}{4}+\Bigg(\frac{4N}{3M^2}D_*^3\Bigg)^{1/n}}-\frac{1}{2}\Bigg]^n \] \[ w_s=\frac{33.94\cdot10^{-6}}{0.98\cdot 0.001}\Bigg[\sqrt{\frac{1}{4}+\Bigg(\frac{4\cdot0.98}{3\cdot33.94^2}16186.5\Bigg)^{1/1.33}}-\frac{1}{2}\Bigg]^{1.33} \] \[ w_s=0.12 \ m/s \]
* Wu and Wang (2006), for d= 0.01 mm;
and assuming,\(S_p=0.7\), g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ M=53.5e^{-0.65S_p};N=5.65e^{-2.5S_p}; n=0.7+0.9S_p \] \[ M=53.5e^{-0.65\cdot 0.7};N=5.65e^{-2.5 \cdot 0.7}; n=0.7+0.9\cdot 0.7 \] \[ M=33.94;N=0.23; n=1.33 \]
\[ D_*^3=\Bigg(\frac{\rho_s}{\rho}-1\Bigg)\cdot g\cdot d^3/\nu^2 \]
\[ D_*^3=\Bigg(\frac{2650}{1000}-1\Bigg)\cdot 9.81\cdot (0.00001)^3/10^{-6^{2}} \]
\[ D_*^3=0.0161865 \]
\[ w_s=\frac{M\nu}{Nd}\Bigg[\sqrt{\frac{1}{4}+\Bigg(\frac{4N}{3M^2}D_*^3\Bigg)^{1/n}}-\frac{1}{2}\Bigg]^n \]
\[ w_s=\frac{33.94\cdot10^{-6}}{0.23\cdot 0.00001}\Bigg[\sqrt{\frac{1}{4}+\Bigg(\frac{4\cdot0.23}{3\cdot 33.94^2}0.0161865\Bigg)^{1/1.33}}-\frac{1}{2}\Bigg]^{1.33} \]
\[ w_s=6.36\cdot10^{-5} \ m/s \]
* Ahrens (2000), for d= 1 mm;
and assuming,g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ w_s=C_1\cdot \Delta\cdot g\cdot d^2/\nu^2+C_t \sqrt{\Delta\cdot g \cdot d} \]
Computing \(\Delta\)
\[ \Delta=\frac{\rho_s}{\rho}-1 \] \[ \Delta=\frac{2650}{1000}-1 \] \[ \Delta=1.65 \]
Computing \(C_1\)
\[ C_1=0.055\cdot tanh\Bigg[12\cdot A^{-0.59}\cdot exp(-0.0004\cdot A)\Bigg] \]
In which A is given by:
\[ A=\Delta\cdot g \cdot d^3/\nu^2 \] \[ A=1.65\cdot 9.81 \cdot (0.001)^3/10^{-6^{2}} \] \[ A=16186.5 \]
Substituting A in \(C_1\) expression:
\[ C_1=0.055\cdot tanh\Bigg[12\cdot16186.5^{-0.59}\cdot exp(-0.0004\cdot 16186.5)\Bigg] \] \[ C_1=3.343957\cdot 10^{-6} \]
Computing \(C_t\)
\[ C_t=1.06\cdot tanh\Bigg[0.016\cdot A^{0.50}\cdot exp(-120/A)\Bigg] \]
\[ C_t=1.023381 \]
Substituting \(C_1\) and \(C_t\) in Ahrens (2000) formula
\[ w_s=3.343957\cdot 10^{-6}\cdot 1.65\cdot 9.81\cdot (0.001)^2/10^{-6}+1.023381 \sqrt{1.65\cdot 9.81 \cdot 0.001} \]
\[ w_s=0.13 \ m/s \]
* Ahrens (2000), for d= 0.01 mm;
and assuming,g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
Updating A :
\[ A=\Delta\cdot g \cdot d^3/\nu^2 \]
\[ A=1.65\cdot 9.81 \cdot (0.00001)^3/10^{-6^{2}} \]
\[ A=0.0161865 \]
Substituting A in \(C_1\) expression:
\[ C_1=0.055\cdot tanh\Bigg[12\cdot 0.0161865^{-0.59}\cdot exp(-0.0004\cdot 0.0161865)\Bigg] \] \[ C_1=0.055 \]
Substituting A in \(C_t\) expression
\[ C_t=1.06\cdot tanh\Bigg[0.016 \cdot 0.0161865^{0.50}\cdot exp(-120/(0.0161865))\Bigg] \]
\[ C_t=0 \]
Substituting \(C_1\) and \(C_t\) in Ahrens (2000) formula
\[ w_s=0.055 \cdot 1.65\cdot 9.81\cdot (0.00001)^2/10^{-6}+1.025491 \sqrt{1.5\cdot 9.81 \cdot 0.00001} \]
\[ w_s=8.90\cdot10^{-5} \ m/s \]
* Jimenez and Madsen (2003), for d= 1 mm;
and assuming,\(S_p=0.7\),g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
\[ \frac{w_s}{\sqrt{\Delta gd}}=\Bigg(C+\frac{B}{S_*}\Bigg)^{-1} \]
Computing \(S_*\))
\[ S_*=\frac{d}{4\nu}\sqrt{\Delta gd} \] \[ S_*=\frac{0.001}{4\cdot10^{-6}}\sqrt{1.65\cdot 9.81\cdot 0.001} \] \[ S_*=31.80654 \]
For Natural sediments (\(S_p=0.7\)) and assuming C=0.954 and B=5.121
\[ \frac{w_s}{\sqrt{1.65\cdot 9.81\cdot 0.001}}=\Bigg(0.954+\frac{5.121}{31.80654}\Bigg)^{-1} \] \[ w_s=0.11 \ m/s \]
* Jimenez and Madsen (2003), for d= 0.01 mm;
and assuming,\(S_p=0.7\),g=9.81 \(m/s^2\), \(\nu=10^{-6}\ m^2/s\), \(\rho=1000 kg/m^3\), \(\rho_s=2650 kg/m^3\)
Computing \(S_*\)
\[ S_*=\frac{d}{4\nu}\sqrt{\Delta gd} \] \[ S_*=\frac{0.00001}{4\cdot10^{-6}}\sqrt{1.65\cdot 9.81\cdot 0.00001} \] \[ S_*=0.03180654 \]
For Natural sediments (\(S_p=0.7\)) and assuming C=0.954 and B=5.121
\[ \frac{w_s}{\sqrt{1.65\cdot 9.81\cdot 0.00001}}=\Bigg(0.954+\frac{5.121}{0.03180654}\Bigg)^{-1} \]
\[ w_s=7.85\cdot 10^{-5} \ m/s \]
PROBLEM 4.10. Draw a graph similar to Fig. 4.15 to illustrate the relation of the standard fall diameter with the nominal diameter and shape factor for naturally worn quartz particles. Use the procedure described in Appendix 4.3.
Solution
Following the guidelines described in Appendix 4.3:
a) Assume a nominal diameter and a Corey shape factor, and calculate the settling velocity using Eq. (4.34), as presented below:
\[ w_s=\frac{M\nu}{Nd}\Bigg[\sqrt{\frac{1}{4}+\Bigg(\frac{4N}{3M^2}D_*^3\Bigg)^{1/n}}-\frac{1}{2}\Bigg]^n \]
b) Submit the calculated settling velocity into Eq. (4.11) to calculate the diameter of the equivalent sphere, which is the fall diameter of the sediment particle.
\[ \frac{w_s}{[g\cdot \nu \cdot (\rho_s-\rho)/\rho]^{1/3}}=\Bigg[\Bigg(\frac{18}{D_*^2}\Bigg)^{0.898(1+0.936D_*)/(1+D_*)}+\Bigg(\frac{0.317}{D_*}\Bigg)^{0.449}\Bigg]^-1.114 \]
A R code was developed to process the described steps (a) and (b), considering the shape factors (\(S_p=1; 0.9; 0.7; 0.5 and 0.3\)). For each shape factors several diameters ranging from 0.05 mm until 10 mm were considered. The code and the final plot are presented below.
obs: The sieve diameter was consider 0.9 of the nominal diameter
#Clear R memory
rm(list=ls())
#Library for inverse function
library(GoFKernel)## Loading required package: KernSmooth## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009#Eq.4.34
ws<-function(dn,sp)
{
M=53.5*exp(-0.65*sp)
N=5.65*exp(-2.5*sp)
n=0.7+0.9*sp
rhos=2650
rho=1000
g=9.81
nu=1/1000000
Dstar=dn*((rhos/rho-1)*g/(nu^2))^(1/3)
ws=(M*nu)/(N*dn)*(sqrt(1/4+((4*N)/(3*(M^2))*Dstar^3)^(1/n))-0.5)^n
return(ws)
}
#Eq.4.11
ws2<-function(d)
{
rhos=2650
rho=1000
g=9.81
nu=1/1000000
Dstar=d*((rhos/rho-1)*g/(nu^2))^(1/3)
ws2=(g*nu*(rhos-rho)/rho)^(1/3)*((18/(Dstar^2))^(0.898*(1+0.936*Dstar)/(1+Dstar))+(0.317/Dstar)^0.449)^(-1.114)
return(ws2)
}
#Define the inverse function of the ws2
f.inv <- inverse(ws2,lower=0.00001,upper=0.1)
dn=seq(from=0.00002,to=0.01,by=0.0005)
#sp=1
w1=0
d1=0
#sp=0.9
w2=0
d2=0
#sp=0.7
w3=0
d3=0
#sp=0.5
w4=0
d4=0
#sp=0.3
w5=0
d5=0
i=1
#for sp=0.9
for(i in 1:length(dn) )
{
w1[i]=ws(dn[i],1)
d1[i]=f.inv(w1[i])
w2[i]=ws(dn[i],0.9)
d2[i]=f.inv(w2[i])
w3[i]=ws(dn[i],0.7)
d3[i]=f.inv(w3[i])
w4[i]=ws(dn[i],0.5)
d4[i]=f.inv(w4[i])
w5[i]=ws(dn[i],0.3)
d5[i]=f.inv(w5[i])
}Plot
#ploting the results
plot((d1*1000),(dn*0.9*1000),type="l",pch=8,col="black",lwd=1,xlab="Fall Diameter [mm]",ylab="Sieve Diameter [mm]",log="xy")
lines((d2*1000),(dn*0.9*1000),type="l",lwd=1.2,pch=19,col="red",lty=2,cex=0.5,log="xy")
lines((d3*1000),(dn*0.9*1000),type="l",lwd=1.2,pch=13,col="blue",lty=3,cex=0.5,log="xy")
lines((d4*1000),(dn*0.9*1000),type="l",lwd=1.2,pch=17,col="darkgreen",lty=4,cex=0.5,log="xy")
lines((d5*1000),(dn*0.9*1000),type="l",lwd=1.2,pch=17,col="orange",lty=5,cex=0.5,log="xy")
#legend
legend("topleft",legend=c("Sp=1","Sp=0.9","Sp=0.7","Sp=0.5","Sp=0.3"),col=c("black","red","blue","darkgreen","orange"), lty=c(1,2,3,4,5),lwd=c(1.2,1.2,1.2,1.2,1.2),cex=1.2,border="black",bty="o",box.lwd = 1,box.col = "black",bg = "white")