Skip to content
Extraits de code Groupes Projets
sediment_functions.py 4,81 ko
Newer Older
  • Learn to ignore specific revisions
  • 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