tuse momentum conservation eq from Kingslake and Ng 2013 - 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 736f50c0a19c6d0ab65821b49e60a0fc010e0f84
(DIR) parent 5b2de54976269ae0d0883f728b9a6337610505f2
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 8 Mar 2017 15:27:06 -0800
use momentum conservation eq from Kingslake and Ng 2013
Diffstat:
M 1d-channel.py | 121 +++++++++++++++++--------------
1 file changed, 66 insertions(+), 55 deletions(-)
---
(DIR) diff --git a/1d-channel.py b/1d-channel.py
t@@ -29,9 +29,10 @@ 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 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
-plot_interval = 20 # Time steps between plots
+print_output_convergence = False # Display convergence statistics during run
+safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
+# plot_interval = 20 # Time steps between plots
+plot_interval = 1 # Time steps between plots
# Physical parameters
rho_w = 1000. # Water density [kg/m^3]
t@@ -44,23 +45,23 @@ 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 = 1e-3 # Sand transported near margin
+m_dot = 2.5e-8
mu_w = 1.787e-3 # Water viscosity [Pa*s]
friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
-# Hewitt 2011 channel flux parameters
+# Channel hydraulic properties
manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
-F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
+#F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
# Channel growth-limit parameters
c_1 = -0.118 # [m/kPa]
c_2 = 4.60 # [m]
# Minimum channel size [m^2], must be bigger than 0
-# S_min = 1e-2
+S_min = 1e-2
# S_min = 1e-1
-S_min = 1.
+# S_min = 1.
# # Initialize model arrays
t@@ -112,9 +113,11 @@ def avg_midpoint(arr):
return (arr[:-1] + arr[1:])/2.
-def channel_water_flux(S, hydro_pot_grad):
- # Hewitt 2011
- return numpy.sqrt(1./F*S**(8./3.)*-hydro_pot_grad)
+def channel_hydraulic_roughness(manning, S, W, theta):
+ # Determine hydraulic roughness assuming that the channel has no channel
+ # floor wedge.
+ l = W + W/numpy.cos(numpy.deg2rad(theta)) # wetted perimeter
+ return manning**2.*(l**2./S)**(2./3.)
def channel_shear_stress(Q, S):
t@@ -263,11 +266,13 @@ def channel_growth_rate_sedflux(Q_t, porosity, s_c):
def update_channel_size_with_limit(S, S_old, dSdt, dt, N_c):
# Damsgaard et al, in prep
S_max = numpy.maximum(
- numpy.maximum((c_1*numpy.maximum(N_c, 0.)/1000. + c_2), 0.)**2./4. *
+ numpy.maximum(
+ 1./4.*(c_1*numpy.maximum(N_c, 0.)/1000. + c_2), 0.)**2. *
numpy.tan(numpy.deg2rad(theta)), S_min)
S = numpy.maximum(numpy.minimum(S_old + dSdt*dt, S_max), S_min)
W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
- return S, W, S_max
+ dSdt = S - S_old # adjust dSdt for clipping due to channel size limits
+ return S, W, S_max, dSdt
def flux_solver(m_dot, ds):
t@@ -286,7 +291,7 @@ def flux_solver(m_dot, ds):
Q[1:] = m_dot*ds[1:] + Q[:-1]
max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
- if output_convergence:
+ if print_output_convergence:
print('it = {}: max_res = {}'.format(it, max_res))
# import ipdb; ipdb.set_trace()
t@@ -298,30 +303,34 @@ def flux_solver(m_dot, ds):
return Q
-def pressure_solver(psi, F, Q, S):
+def pressure_solver(psi, f, Q, S):
# Iteratively find new water pressures
- # dP_c/ds = psi - FQ^2/S^{8/3}
+ # dP_c/ds = psi - f*rho_w*g*Q^2/S^{8/3} (Kingslake and Ng 2013)
it = 0
max_res = 1e9 # arbitrary large value
- while max_res > tol_P_c or it < Ns*40:
+ while max_res > tol_P_c or it < Ns:
P_c_old = P_c.copy()
- # Upwind finite differences
- P_c[:-1] = -psi[:-1]*ds[:-1] \
- + F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
- + P_c[1:] # Upstream
+ # P_downstream = P_upstream + dP
+ # P_c[1:] = P_c[:-1] \
+ # + psi[:-1]*ds[:-1] \
+ # - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
# Dirichlet BC (fixed pressure) at terminus
P_c[-1] = 0.
- # von Neumann BC (no gradient = no flux) at s=0
- P_c[0] = P_c[1]
+ # P_upstream = P_downstream - dP
+ P_c[:-1] = P_c[1:] \
+ - psi[:-1]*ds[:-1] \
+ + f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
+ # + psi[:-1]*ds[:-1] \
+ # - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
- if output_convergence:
+ if print_output_convergence:
print('it = {}: max_res = {}'.format(it, max_res))
if it >= max_iter:
t@@ -335,9 +344,9 @@ def pressure_solver(psi, F, Q, S):
def plot_state(step, time, S_, S_max_, title=True):
# Plot parameters along profile
fig = plt.gcf()
- fig.set_size_inches(3.3*1.1, 3.3*1.1)
+ fig.set_size_inches(3.3*1.1, 3.3*1.1*1.5)
- ax_Pa = plt.subplot(2, 1, 1) # axis with Pascals as y-axis unit
+ ax_Pa = plt.subplot(3, 1, 1) # axis with Pascals as y-axis unit
# ax_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
# ax_Pa.plot(s/1000., N/1000., '--r', label='$N$')
ax_Pa.plot(s_c/1000., N_c/1e6, '-k', label='$N$')
t@@ -355,25 +364,29 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_Pa.set_ylabel('[MPa]')
ax_m3s.set_ylabel('[m$^3$/s]')
- ax_m2 = plt.subplot(2, 1, 2, sharex=ax_Pa)
+ ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
+ ax_m3s_sed.plot(s_c/1000., Q_g, ':', label='$Q_g$')
+ ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_s$')
+ ax_m3s_sed.plot(s_c/1000., Q_t, '--', label='$Q_t$')
+ ax_m3s_sed.set_ylabel('[m$^3$/s]')
+ ax_m3s_sed.legend(loc=2)
+
+ ax_m2 = plt.subplot(3, 1, 3, sharex=ax_Pa)
ax_m2.plot(s_c/1000., S_, '-k', label='$S$')
ax_m2.plot(s_c/1000., S_max_, '--', color='#666666', label='$S_{max}$')
# ax_m.semilogy(s_c/1000., S, '-k', label='$S$')
# ax_m.semilogy(s_c/1000., S_max, '--k', label='$S_{max}$')
- ax_ms = ax_m2.twinx()
+ ax_m2s = 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_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_m2s.plot(s_c/1000., dSdt, ':', label='$\partial S/\partial t$')
ax_m2.legend(loc=2)
- ax_ms.legend(loc=3)
+ ax_m2s.legend(loc=3)
ax_m2.set_xlabel('$s$ [km]')
ax_m2.set_ylabel('[m$^2$]')
- ax_ms.set_ylabel('[m/s]')
+ ax_m2s.set_ylabel('[m$^2$/s]')
ax_Pa.set_xlim([s.min()/1000., s.max()/1000.])
t@@ -412,8 +425,8 @@ hydro_pot_grad = gradient(hydro_pot, s)
N_c = avg_midpoint(N)
H_c = avg_midpoint(H)
-# Find fluxes in channel segments [m^3/s]
-Q = channel_water_flux(S, hydro_pot_grad)
+# Find water flux
+Q = flux_solver(m_dot, ds)
# Water-pressure gradient from geometry [Pa/m]
psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
t@@ -428,18 +441,22 @@ time = 0.
step = 0
while time <= t_end:
- # print('@ @ @ step ' + str(step))
-
+ # Determine time step length from water flux
dt = find_new_timestep(ds, Q, S)
+ # Display current simulation status
print_status_to_stdout(time, dt)
it = 0
- max_res = 1e9 # arbitrary large value
+
+ # Initialize the maximum normalized residual for S to an arbitrary large
+ # value
+ max_res = 1e9
S_old = S.copy()
# Iteratively find solution, do not settle for less iterations than the
- # number of nodes
+ # number of nodes to make sure information has had a chance to pass through
+ # the system
while max_res > tol_Q or it < Ns:
S_prev_it = S.copy()
t@@ -447,43 +464,37 @@ while time <= t_end:
# Find average shear stress from water flux for each channel segment
tau = channel_shear_stress(Q, S)
- # Find sediment erosion and deposition rates for each channel segment
- # e_dot = channel_erosion_rate(tau)
- # d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
-
+ # Determine sediment fluxes for each size fraction
f_g = 1./f_s # gravel volume fraction is reciprocal to sand
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
- # TODO: Update f_s from fluxes
-
- # Determine change in channel size for each channel segment
- # dSdt = channel_growth_rate(e_dot, d_dot, W)
- # Use backward differences and assume dS/dt=0 in first segment
+ # Determine change in channel size for each channel segment.
+ # Use backward differences and assume dS/dt=0 in first segment.
dSdt[1:] = channel_growth_rate_sedflux(Q_t, porosity, s_c)
# Update channel cross-sectional area and width according to growth
# rate and size limit for each channel segment
- S, W, S_max = update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
+ S, W, S_max, dSdt = \
+ update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
# Find new water fluxes consistent with mass conservation and local
# meltwater production (m_dot)
Q = flux_solver(m_dot, ds)
- # Find the corresponding sediment flux
- # Q_b = bedload_sediment_flux(
- # Q_s = suspended_sediment_flux(c_bar, Q, S)
+ # Find hydraulic roughness
+ f = channel_hydraulic_roughness(manning, S, W, theta)
# Find new water pressures consistent with the flow law
- P_c = pressure_solver(psi, F, Q, S)
+ P_c = pressure_solver(psi, f, Q, S)
# Find new effective pressure in channel segments
N_c = rho_i*g*H_c - P_c
# Find new maximum normalized residual value
max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
- if output_convergence:
+ if print_output_convergence:
print('it = {}: max_res = {}'.format(it, max_res))
# import ipdb; ipdb.set_trace()