tworking example with Wilcock two-phase sediment transport - granular-channel-hydro - subglacial hydrology model for sedimentary channels
(HTM) git clone git://src.adamsgaard.dk/granular-channel-hydro
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit 5b2de54976269ae0d0883f728b9a6337610505f2
(DIR) parent 5de4d190454b1c94b1982063f442efcfc30ee649
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 6 Mar 2017 15:07:13 -0800
working example with Wilcock two-phase sediment transport
Diffstat:
M 1d-channel.py | 79 ++++++++++++-------------------
1 file changed, 31 insertions(+), 48 deletions(-)
---
(DIR) diff --git a/1d-channel.py b/1d-channel.py
t@@ -21,16 +21,16 @@ import sys
# # Model parameters
-Ns = 25 # Number of nodes [-]
-# Ls = 100e3 # Model length [m]
-Ls = 1e3 # Model length [m]
-total_days = 60. # Total simulation time [d]
+Ns = 25 # Number of nodes [-]
+# Ls = 100e3 # Model length [m]
+Ls = 1e3 # Model length [m]
+total_days = 60. # Total simulation time [d]
t_end = 24.*60.*60.*total_days # Total simulation time [s]
-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
+tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
+tol_P_c = 1e-3 # Tolerance criteria for the norm. 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
+safety = 0.1 # Safety factor ]0;1] for adaptive timestepping
plot_interval = 20 # Time steps between plots
# Physical parameters
t@@ -44,27 +44,10 @@ D_g = 1. # Mean grain size in gravel fraction (> 2 mm)
D_s = 0.01 # Mean grain size in sand fraction (<= 2 mm)
# Water source term [m/s]
-m_dot = 7.93e-11
-# m_dot = 1.0e-7
-# m_dot = 2.0e-6
-# m_dot = 4.5e-7
-# m_dot = 5.79e-5
-# m_dot = 5.0e-6
-# m_dot = 1.8/(1000.*365.*24.*60.*60.) # Whillan's melt rate from Joughin 2004
-
-# Walder and Fowler 1994 sediment transport parameters
-K_e = 6.0 # Erosion constant [-], disabled when 0.0
-# K_d = 6.0 # Deposition constant [-], disabled when 0.0
-K_d = 0.01*K_e # Deposition constant [-], disabled when 0.0
-alpha = 1e5 # 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 2017)
-# tau_c = 0.
+m_dot = 1e-3 # Sand transported near margin
+
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 2017)
+friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
# Hewitt 2011 channel flux parameters
manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
t@@ -75,8 +58,8 @@ c_1 = -0.118 # [m/kPa]
c_2 = 4.60 # [m]
# Minimum channel size [m^2], must be bigger than 0
-# S_min = 1e-1
# S_min = 1e-2
+# S_min = 1e-1
S_min = 1.
t@@ -111,7 +94,7 @@ 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
-Q_t = numpy.zeros_like(S) # Sediment flux where D <= 2 mm [m3/s]
+Q_t = numpy.zeros_like(S) # Total sediment flux [m3/s]
Q_s = numpy.zeros_like(S) # Sediment flux where D <= 2 mm [m3/s]
Q_g = numpy.zeros_like(S) # Sediment flux where D > 2 mm [m3/s]
f_s = numpy.ones_like(S)*sand_fraction # Initial sediment fraction of sand [-]
t@@ -135,9 +118,9 @@ def channel_water_flux(S, hydro_pot_grad):
def channel_shear_stress(Q, S):
- # Weertman 1972, Walder and Fowler 1994
+ # Determine mean wall shear stress from Darcy-Weisbach friction loss
u_bar = Q/S
- return 1./8.*froude*rho_w*u_bar**2.
+ return 1./8.*friction_factor*rho_w*u_bar**2.
def channel_erosion_rate(tau):
t@@ -164,7 +147,7 @@ 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.)
+ / (rho_w*friction_factor))*h**(3./2.)
return v_s/epsilon*c_bar[ix]
t@@ -223,9 +206,15 @@ def channel_sediment_flux_sand(tau, W, f_s, D_s):
shields_stress = tau/((rho_s - rho_w)*g*D_s)
# import ipdb; ipdb.set_trace()
- return 11.2*f_s*W/((rho_s - rho_w)/rho_w*g) \
+ Q_c = 11.2*f_s*W/((rho_s - rho_w)/rho_w*g) \
* (tau/rho_w)**1.5 \
- * (1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))**4.5
+ * numpy.maximum(0.0,
+ (1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))
+ )**4.5
+
+ # The above relation gives 'nan' values for low values of tau
+
+ return Q_c
def channel_sediment_flux_gravel(tau, W, f_g, D_g):
t@@ -249,7 +238,8 @@ def channel_sediment_flux_gravel(tau, W, f_g, D_g):
# From Wilcock 2001, eq. 3
Q_g = 11.2*f_g*W/((rho_s - rho_w)/rho_w*g) \
* (tau/rho_w)**1.5 \
- * (1.0 - 0.846*ref_shear_stress/shields_stress)**4.5
+ * numpy.maximum(0.0,
+ (1.0 - 0.846*ref_shear_stress/shields_stress))**4.5
# From Wilcock 2001, eq. 4
I = numpy.nonzero(ref_shear_stress/shields_stress < 1.)
t@@ -301,19 +291,13 @@ def flux_solver(m_dot, ds):
# import ipdb; ipdb.set_trace()
if it >= max_iter:
- raise Exception('t = {}, step = {}:'.format(time, step) +
+ raise Exception('t = {}, step = {}: '.format(time, step) +
'Iterative solution not found for Q')
it += 1
return Q
-def suspended_sediment_flux(c_bar, Q, S):
- # Find the fluvial sediment flux through the system
- # Q_s = c_bar * u * S, where u = Q/S
- return c_bar*Q
-
-
def pressure_solver(psi, F, Q, S):
# Iteratively find new water pressures
# dP_c/ds = psi - FQ^2/S^{8/3}
t@@ -380,9 +364,9 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_ms = ax_m2.twinx()
# ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
# ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
- ax_ms.plot(s_c/1000., Q_t, label='$Q_t$')
ax_ms.plot(s_c/1000., Q_g, label='$Q_g$')
ax_ms.plot(s_c/1000., Q_s, label='$Q_s$')
+ ax_ms.plot(s_c/1000., Q_t, '--', label='$Q_t$')
# TODO: check units on sediment fluxes: m2/s or m3/s ?
ax_m2.legend(loc=2)
t@@ -471,7 +455,6 @@ while time <= t_end:
Q_s = channel_sediment_flux_sand(tau, W, f_s, D_s)
Q_g = channel_sediment_flux_gravel(tau, W, f_g, D_g)
Q_t = Q_s + Q_g
- break
# TODO: Update f_s from fluxes
t@@ -490,7 +473,7 @@ while time <= t_end:
# Find the corresponding sediment flux
# Q_b = bedload_sediment_flux(
- Q_s = suspended_sediment_flux(c_bar, Q, S)
+ # Q_s = suspended_sediment_flux(c_bar, Q, S)
# Find new water pressures consistent with the flow law
P_c = pressure_solver(psi, F, Q, S)
t@@ -505,8 +488,8 @@ while time <= t_end:
# import ipdb; ipdb.set_trace()
if it >= max_iter:
- raise Exception('t = {}, step = {}:'.format(time, step) +
- 'Iterative solution not found for Q')
+ raise Exception('t = {}, step = {}: '.format(time, step) +
+ 'Iterative solution not found')
it += 1
# Generate an output figure for every n time steps