sediment_functions.py

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