In [1]:

    

from math import sqrt, exp
from scipy.constants import g
import numpy as np
import matplotlib.pyplot as plt


    

In [2]:

    

def kinematic_viscosity(t=10):
    return 1.792e-6 / (1 + pow(t / 25, 1.165))


    

In [3]:

    

def shields_curve_equation(d, r_p=2.65, r_w=1.0, v=1.40e-6):
    """
    Calculates shields curve equation mathematically.
    :param d: Particle diameter [m]
    :param r_p: Density of particle [t/m3]
    :param r_w: Density of water [t/m3]
    :param v: Kinematic viscosity [m2/s]
    :return: array of (critical shields number, critical shear stress, shear velocity, reynolds number)
    """
    dimensionless_dia = pow((r_p / r_w - 1.0) * g / pow(v, 2.0), 1.0 / 3.0) * d
    critical_shields_number = 0.23 / dimensionless_dia + 0.054 * (1.0 - exp(-pow(dimensionless_dia, 0.85) / 23.0))
    critical_shear_stress = critical_shields_number * (r_p / r_w - 1.0) * g * r_w * d
    shear_velocity = sqrt(critical_shear_stress / r_w)
    reynolds_number = shear_velocity * d / v
    return critical_shields_number, critical_shear_stress, shear_velocity, reynolds_number


    

In [4]:

    

print(shields_curve_equation(0.001, 2.65, 1.00,1.33345e-6))


    

    




(0.034643948029674375, 0.0005605727703595902, 0.02367641802215002, 17.755759887622347)

In [5]:

    

def shields_diagram(r_p=2.65, r_w=1.0, v=1.40e-6):
    diameters = np.linspace(0.01, 100.0, 2000, dtype=np.float)
    print(diameters)
    critical_shields_number = np.fromiter((shields_curve_equation(d / 1000.0, r_p, r_w, v)[0] for d in diameters),
                                          np.float)
    reynolds_number = np.fromiter((shields_curve_equation(d / 1000.0, r_p, r_w, v)[3] for d in diameters), np.float)
    fig, ax = plt.subplots()
    ax.set_xlabel(r'$Reynolds\;Sayısı\;[Re_{*c}]$')
    ax.set_ylabel(r'$Shields\;Sayısı\;[\tau_{*c}]$')
    ax.loglog(reynolds_number, critical_shields_number, basex=10)
    ax.grid(True, which='both')
    ax.set_ylim(bottom=0.01)
    ax.set_ylim(top=1.0)
    ax.set_xlim(left=0.1)
    ax.set_xlim(right=1000.0)
    plt.show()


    

In [6]:

    

shields_diagram(2.65, 1.00,1.33345e-6)


    

    




[1.0000000e-02 6.0020010e-02 1.1004002e-01 ... 9.9899960e+01 9.9949980e+01
 1.0000000e+02]






    

In [7]:

    

def shields_diagram_critical_shear_stress_velocity(r_p=2.65, r_w=1.0, v=1.40e-6):
    diameters = np.linspace(0.01, 100.0, 2000, dtype=np.float)
    shear_stress = np.fromiter((shields_curve_equation(d / 1000.0, r_p, r_w, v)[1] * g * 1000 for d in diameters),
                               np.float)
    shear_velocity = np.fromiter((shields_curve_equation(d / 1000.0, r_p, r_w, v)[2] for d in diameters), np.float)
    fig, ax1 = plt.subplots()
    ax1.set_xlabel('$D\;[mm]$')
    ax1.set_ylabel('$u_{*c}\;[m/s]$')
    ax1.loglog(diameters, shear_velocity, basex=10, color='r', linestyle='--', label='$u_{*c}$')
    ax1.tick_params(axis='y')
    ax1.grid(True, which='major')
    ax1.grid(True, which='minor')
    ax1.set_ylim(bottom=0.001)
    ax1.set_ylim(top=1.0)
    ax1.legend(loc=2)

    ax2 = ax1.twinx()
    ax2.set_ylabel(r'$\tau_{*c}\,[\dfrac{N}{m^2}]$')
    ax2.loglog(diameters, shear_stress, basex=10, label=r'$\tau_{*c}$')
    ax2.tick_params(axis='y')
    ax2.set_ylim(bottom=1.0)
    ax2.set_ylim(top=1000.0)
    ax2.legend(loc=4)
    fig.tight_layout()
    plt.show()


    

In [8]:

    

shields_diagram_critical_shear_stress_velocity(2.65, 1.00,1.33345e-6)