Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
import numpy as np
def get_chezy_coeff(h, ks):
# Dimensionless Chézy coefficient for rough turbulent flow (=1/sqrt(Cd))
# based on the empirical White-Colebrook eq. for large estuaries (high Re)
# In the hydro, the bulk drag Cd = 2.5-3 => min value of Chezy is 1/sqrt(2.5e-3)
x = np.array(2.5 * np.log(11 * h / ks))
x = np.append(x, 1.0 / np.sqrt(2.5e-3))
return np.min(x)
def get_skin_roughness_stress(U, h, rho_w, D90):
# Bottom shear stress due to skin roughness
ks = 3 * D90
Cs = get_chezy_coeff(h, ks)
return rho_w * (U ** 2) / (Cs ** 2)
def get_form_roughness_stress(U, h, rho_w, eta):
# Bottom shear stress due to bed form roughness
# eta = bed form roughness height
Cs = get_chezy_coeff(h, eta)
return rho_w * (U ** 2) / (Cs ** 2)
def get_etab(M, h, D50):
# Equilibrium bed form height eta_b according to Van Rijn (1984)
# eta_b = 0 if D50 < 0.05mm, if M<1 or if M>24
# if (D50<5e-5 or M<=1 or M>24):
# return 0.0
# else:
# return 0.11 * h * ((D50/h)**0.3) * (1-np.exp(-0.5*(M-1))) * (24-M)
res = np.zeros_like(M)
if D50 < 5e-5:
return res + 1e-2
else:
res = 0.11 * h * ((D50 / h) ** 0.3) * (1 - np.exp(-0.5 * (M - 1))) * (24 - M)
res[np.logical_or(M <= 1, M > 24)] = 0
return res + 1e-2
def get_ws(Dgr, D, nu):
# Particle fall velocity ws (m/s) from Soulsby (1997)
if (np.array(Dgr)).size == 1:
if Dgr >= 0.672:
tmp = np.sqrt(107.33 + 1.049 * (Dgr ** 3)) - 10.36
else:
tmp = 0.0077 * (Dgr ** 2)
else:
tmp = 0.0077 * (Dgr ** 2)
index = Dgr >= 0.672
tmp[index] = np.sqrt(107.33 + 1.049 * (Dgr[index] ** 3)) - 10.36
return nu * tmp / D
def get_zs(h, ws, ust):
# Height of the centroid of the suspended load
# ws = particle fall get_bed_shear_velocity
# ust = shear velocity from both shear and form drag
kappa = 0.4
p = 1.2 * np.log(ws / (kappa * ust)) - 0.4
return 0.0398 * h * 10 ** (-1.08 * np.tanh(p))
def get_us(ust, zs, ks):
# Advection speed of suspended load us (m/s)
# ust = shear velocity associated with form drag only
# ks = bed form roughness height
return 2.5 * ust * np.log(30 * zs / ks)
########
def get_theta(tau, rho_w, rho_s, D):
# Shields parameter
# tau = skin shear stress
g = 9.81
s = rho_s / rho_w
return tau / (rho_w * g * (s - 1) * D)
def get_Dgr(rho_w, rho_s, nu, D50):
# Sediment dimensionless grain size
# D50 = grain size at which 50% of the sediment is finer
g = 9.81
s = rho_s / rho_w
return D50 * ((s - 1) * g / (nu ** 2)) ** (1 / 3)
def get_thetac(Dgr):
# Critical Shields parameter from Soulsby and Whitehouse (1997)
return 0.3 / (1 + 1.2 * Dgr) + 0.055 * (1 - np.exp(-0.02 * Dgr))
def get_bed_shear_velocity(U, Cs):
return U / Cs
def get_ub(ust, M):
# Advection speed of the bedload from Engelund and Fredsoe (1976)
# ust = shear velocity associated with skin friction only
# M = dimensionless mobility parameter
ub = ust * (10 - 7 / np.sqrt(M))
return np.maximum(ust * (10 - 7 / np.sqrt(M)), np.zeros_like(M))
def get_sediment_param(sed_str="medium_sand"):
mu = 1e-6 # micrometer
if sed_str == "coarse_sand":
D90 = 1000 * mu
D50 = 750 * mu
D35 = 500 * mu
rho_s = 1900.0
elif sed_str == "medium_sand":
D90 = 500 * mu
D50 = 375 * mu
D35 = 250 * mu
rho_s = 1900.0
elif sed_str == "fine_sand":
D90 = 250 * mu
D50 = 190 * mu
D35 = 125 * mu
rho_s = 1900.0
elif sed_str == "coarse_silt":
D90 = 31 * mu
D50 = 23 * mu
D35 = 16 * mu
rho_s = 1700.0
elif sed_str == "medium_silt":
D90 = 16 * mu
D50 = 12 * mu
D35 = 8 * mu
rho_s = 1700.0
elif sed_str == "fine_silt":
D90 = 8 * mu
D50 = 6 * mu
D35 = 4 * mu
rho_s = 1700.0
elif sed_str == "dissolved":
D90 = 0
D50 = 0
D35 = 0
rho_s = 1000.0
else:
print("Sediment type unknown")
exit(0)
return D90, D50, D35, rho_s
def get_sed_size_density(min_size, max_size, n):
"""generate random size and related density in given range for each particle"""
mu = 1e-6
D = np.random.uniform(min_size * mu, max_size * mu, n) # in micrometers
phi = -np.log2(D * 1e3) # This equation expects D in mm
rho_s = (
2.374 - 0.175 * phi + 0.008 * phi * phi
) * 1e3 # multiplied by 1e3 to get kg/m3 instead of g/cm3
return D, rho_s
def resuspensionProbability(theta, theta_crit):
"""Return a probability of resupsension for a settled particle for which M>1."""
proba = np.random.uniform(0, (theta / theta_crit) ** 2)
if proba > 0.8:
return True
else:
return False