commit 26bc00c665916c00d8fc195170f6df2062d77332
parent 8119c0ae2e91d175f51c319d32a10abc108ebcd3
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Thu, 2 Feb 2017 12:31:15 -0800
try different deposition laws, d_dot > e_dot
Diffstat:
| M | 1d-channel.py | | | 67 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------- |
1 file changed, 56 insertions(+), 11 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -28,6 +28,7 @@ tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
max_iter = 1e2*Ns # Maximum number of solver iterations before failure
output_convergence = False # Display convergence statistics during run
+safety = 0.1 # Safety factor ]0;1] for adaptive timestepping
# Physical parameters
rho_w = 1000. # Water density [kg/m^3]
@@ -38,20 +39,22 @@ theta = 30. # Angle of internal friction in sediment [deg]
# Water source term [m/s]
#m_dot = 7.93e-11
-#m_dot = 4.5e-8
-m_dot = 5.79e-5
+m_dot = 4.5e-7
+#m_dot = 5.79e-5
# Walder and Fowler 1994 sediment transport parameters
K_e = 0.1 # Erosion constant [-], disabled when 0.0
-K_d = 6.0 # Deposition constant [-], disabled when 0.0
+#K_d = 6.0 # Deposition constant [-], disabled when 0.0
+K_d = 1e-1 # Deposition constant [-], disabled when 0.0
+alpha = 2e5 # Geometric correction factor (Carter et al 2017)
#D50 = 1e-3 # Median grain size [m]
#tau_c = 0.5*g*(rho_s - rho_i)*D50 # Critical shear stress for transport
d15 = 1e-3 # Characteristic grain size [m]
-tau_c = 0.025*d15*g*(rho_s - rho_i) # Critical shear stress (Carter 2016)
+tau_c = 0.025*d15*g*(rho_s - rho_i) # Critical shear stress (Carter 2017)
#tau_c = 0.
mu_w = 1.787e-3 # Water viscosity [Pa*s]
froude = 0.1 # Friction factor [-]
-v_s = d15**2.*g*2.*(rho_s - rho_i)/(9.*mu_w) # Settling velocity (Carter 2016)
+v_s = d15**2.*g*2.*(rho_s - rho_i)/(9.*mu_w) # Settling velocity (Carter 2017)
# Hewitt 2011 channel flux parameters
manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
@@ -91,7 +94,7 @@ N_c = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
e_dot = numpy.zeros_like(S) # Sediment erosion rate in channel segments [m/s]
d_dot = numpy.zeros_like(S) # Sediment deposition rate in chan. segments [m/s]
-c_bar = numpy.zeros_like(S) # Vertically integrated sediment content [m]
+c_bar = numpy.zeros_like(S) # Vertically integrated sediment concentration [-]
tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
res = numpy.zeros_like(S) # Solution residual during solver iterations
@@ -118,25 +121,64 @@ def channel_shear_stress(Q, S):
def channel_erosion_rate(tau):
# Parker 1979, Walder and Fowler 1994
- return K_e*v_s*(tau - tau_c).clip(0.)/(g*(rho_s - rho_w)*d15)
+ #return K_e*v_s*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)
+ # Carter et al 2017
+ return K_e*v_s/alpha*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)
def channel_deposition_rate_kernel(tau, c_bar, ix):
# Parker 1979, Walder and Fowler 1994
- return K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+ #result = K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+
+ # Carter et al. 2017
+ result = K_d*v_s/alpha*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+
+ print('tau[{}] = {}'.format(ix, tau[ix]))
+ print('c_bar[{}] = {}'.format(ix, c_bar[ix]))
+ print('e_dot[{}] = {}'.format(ix, e_dot[ix]))
+ print('d_dot[{}] = {}'.format(ix, result))
+ print('')
+
+ return result
+
+def channel_deposition_rate_kernel_ng(c_bar, ix):
+ # Ng 2000
+ h = W[ix]/2.*numpy.tan(numpy.deg2rad(theta))
+ epsilon = numpy.sqrt((psi[ix] - (P_c[ix] - P_c[ix - 1])/ds[ix])\
+ /(rho_w*froude))*h**(3./2.)
+ return v_s/epsilon*c_bar[ix]
def channel_deposition_rate(tau, c_bar, d_dot, Ns):
# Parker 1979, Walder and Fowler 1994
# Find deposition rate from upstream to downstream, margin at is=0
+
+ print("\n## Before loop:")
+ print(c_bar)
+ print(d_dot)
+ print('')
+
+
# No sediment deposition at upstream end
c_bar[0] = 0.
d_dot[0] = 0.
for ix in numpy.arange(1, Ns - 1):
# Net erosion in upstream cell
- c_bar[ix] += numpy.maximum((e_dot[ix - 1] - d_dot[ix - 1])*dt, 0.)
+ #c_bar[ix] = numpy.maximum((e_dot[ix-1] - d_dot[ix-1])*dt*ds[ix-1], 0.)
+ c_bar[ix] = numpy.maximum(
+ W[ix - 1]*ds[ix - 1]*rho_s/rho_w*
+ (e_dot[ix - 1] - d_dot[ix - 1])/Q[ix - 1]
+ , 0.)
+
+ #d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
+ d_dot[ix] = channel_deposition_rate_kernel_ng(c_bar, ix)
+
+
+ print("\n## After loop:")
+ print(c_bar)
+ print(d_dot)
+ print('')
- d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
return d_dot, c_bar
@@ -260,7 +302,6 @@ def plot_state(step, time):
def find_new_timestep(ds, Q, S):
# Determine the timestep using the Courant-Friedrichs-Lewy condition
- safety = 0.2
dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
if dt < 1.0:
@@ -299,6 +340,8 @@ plot_state(-1, 0.0)
time = 0.; step = 0
while time <= t_end:
+ #print('@ @ @ step ' + str(step))
+
dt = find_new_timestep(ds, Q, S)
print_status_to_stdout(time, dt)
@@ -333,8 +376,10 @@ while time <= t_end:
plot_state(step, time)
+ #import ipdb; ipdb.set_trace()
if step > 0:
break
+
# Update time
time += dt
step += 1
