First proper commit for Orifice Sizing .ipynb notebook with interactive plots

Author: ScoobyDew

HEMCalcs.py

@@ -0,0 +1,87 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import CoolProp.CoolProp as CP
+from matplotlib.cm import ScalarMappable
+from matplotlib.colors import Normalize
+from mpl_toolkits.mplot3d import Axes3D
+
+def A(d):
+    return np.pi * (d / 2) ** 2
+
+def mHEM(A, rho2, h1, h2, N, Cd=0.66):
+        return N * A * Cd * rho2 * np.sqrt(2*(h1 - h2))
+
+
+T1 = 10 # Upstream temperature (C)
+P2 = 20e5 # Downstream pressure (Pa)
+N = 12 # Number of orrifaces
+
+def HEM_CP(T1, P2, subst='NitrousOxide'):
+    h1 = CP.PropsSI('H', 'T', T1 + 273.15, 'Q', 0, subst)
+    s1 = CP.PropsSI('S', 'T', T1 + 273.15, 'Q', 0, subst)
+
+    # Find downstream enthalpy for liquid with upstream entropy and downstream pressure
+    h2 = CP.PropsSI('H', 'P', P2, 'S', s1, subst)
+
+    rho2 = CP.PropsSI('D', 'P', P2, 'S', s1, subst)
+    return h1, h2, rho2
+
+def plotting(d, T, N):
+    h1, h2, rho2 = HEM_CP(T, P2)
+
+    m_HEM_vectorized = np.vectorize(mHEM)
+    return m_HEM_vectorized(A(d), rho2, h1, h2, N)
+
+def HEMmassflowrate():
+    temps = np.linspace(-10, 32, 500)
+    d = np.linspace(0.1, 2.5, 100) / 1000
+
+
+    plt.figure()
+    # set colors to viridis colormap
+    colormap = plt.cm.viridis
+    color = iter(colormap(np.linspace(0, 1, len(temps))))
+
+    for T in temps:
+        m = plotting(d, T, N)
+        plt.plot(d * 1000, m, label=f'{T:.1f} C', color=next(color), linewidth=.3, alpha=1)
+
+    sm = ScalarMappable(cmap=colormap, norm=Normalize(vmin=min(temps), vmax=max(temps)))
+    sm.set_array([])
+    ax = plt.gca()
+    cbar = plt.colorbar(sm, ax=ax)
+    cbar.set_label('Upstream Temperature (C)')
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Mass Flow Rate (kg/s)')
+    plt.title('Mass Flow Rate vs Diameter vs Temperature')
+    plt.show()
+
+
+def NormalisedHEM():
+    # Plot Mass Flow Rate vs Upstream temperature for d =0.5mm and N = 1
+    temps = np.linspace(-10, 70, 500)
+    d = 1
+    N = 1
+
+    m = plotting(d, temps, N) / 0.66
+    plt.figure()
+    plt.plot(temps, m/1e4)
+    plt.xlabel('Upstream Temperature (C)')
+    plt.ylabel(r'Normalised mass Flow Rate $\left( \frac{kg}{m^2} \right)$ x10$^{-4}$')
+    plt.title('Mass Flow Rate vs Upstream Temperature')
+    plt.tight_layout()
+    plt.show()
+
+# Plot Mass Flow Rate vs diameter(0 to 2) for ambient temperature of 20degC and N = 12
+temps = 20
+d = np.linspace(0.1, 2.5, 100) / 1000
+N = 12
+
+m = plotting(d, temps, N)
+plt.figure()
+plt.plot(d * 1000, m)
+plt.xlabel('Diameter (mm)')
+plt.ylabel('Mass Flow Rate (kg/s)')
+plt.title('Mass Flow Rate vs Diameter')
+plt.tight_layout()
+plt.show()

NHNECalcs.py

@@ -0,0 +1,72 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import CoolProp.CoolProp as CP
+import ipywidgets as widgets
+
+from HEMCalcs import plotting, A
+from SPICalcs import CalcCPI, m_CPI
+
+def NHNEPlot(T, N, kap):
+    d = np.linspace(0.1, 6, 1000) / 1000
+    Pc = 20e5
+
+    mNom = 1.36
+    mHEM = plotting(d, T, N)
+    mSPI = CalcCPI(T, 'NitrousOxide', Pc, N, d)[0]
+    mkappa = (1 / (1 + kap)) * mSPI + (1 - (kap / (1 + kap))) * mHEM
+
+    plt.plot(d * 1000, mHEM, label=r'$\dot{m}_{HEM}$')
+    plt.plot(d * 1000, mSPI, label=r'$\dot{m}_{SPI}$')
+    plt.plot(d * 1000, mkappa, label=f'$\\dot{{m}}_{{\\kappa}}$, $\\kappa={kap}$')
+    plt.axhline(mNom, color='r', linestyle='--', label=f'$\\dot{{m}}_{{Design}}={mNom}$ kg/s')
+
+    plt.fill_between(d * 1000, mHEM, mSPI, color='lightgrey', alpha=0.5)
+
+    idx_HEM = np.argmin(np.abs(mHEM - mNom))
+    idx_SPI = np.argmin(np.abs(mSPI - mNom))
+    idx_kappa = np.argmin(np.abs(mkappa - mNom))
+
+    d_HEM = d[idx_HEM] * 1000
+    d_SPI = d[idx_SPI] * 1000
+    d_kappa = d[idx_kappa] * 1000
+
+    plt.plot(d_HEM, mNom, color='tab:blue', marker='o', label=f'$d_{{HEM}}={d_HEM:.2f}$ mm')
+    plt.plot(d_SPI, mNom, color='tab:orange', marker='o', label=f'$d_{{SPI}}={d_SPI:.2f}$ mm')
+    plt.plot(d_kappa, mNom, color='tab:green', marker='o', label=f'$d_{{\\kappa}}={d_kappa:.2f}$ mm')
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Mass Flow Rate (kg/s)')
+    plt.title(f'Mass Flow Rate vs Diameter\n {N} Orifices at $T_{{tank}} =${T}°C')
+    plt.xlim(0, 2.5)
+    plt.ylim(0, 2)
+    plt.tight_layout()
+    plt.legend()
+
+
+# Create sliders for the temperature, the number of orifices and the kappa value
+def NHNESliders():
+    T_slider = widgets.FloatSlider(
+        value=30,
+        min=-10,
+        max=60,
+        step=0.5,
+        description=r'$T_{\text{tank}}$ (°C)',
+        continuous_update=False
+    )
+    N_slider = widgets.IntSlider(
+        value=12,
+        min=1,
+        max=50,
+        step=1,
+        description='N',
+        continuous_update=False
+    )
+    kap_slider = widgets.FloatSlider(
+        value=1.4,
+        min=0,
+        max=2,
+        step=0.01,
+        description=r'$\kappa$',
+        continuous_update=False
+    )
+    return T_slider, N_slider, kap_slider
+

NitrousDensityPlots.py

@@ -0,0 +1,174 @@
+import CoolProp.CoolProp as CP
+import matplotlib.pyplot as plt
+import numpy as np
+import ipywidgets as widgets
+from IPython.display import display
+from matplotlib.ticker import FuncFormatter
+
+# Define the subst
+substance = 'NitrousOxide'
+
+# Define temperature range (in Kelvin)
+T_min = CP.PropsSI(substance, 'Tmin') + 0.01  # slightly above the minimum temperature
+T_max = CP.PropsSI(substance, 'Tcrit')        # critical temperature
+temperatures = np.linspace(T_min, T_max, int(1e3))
+
+def get_saturations(T):
+    liquid_density = []
+    vapor_density = []
+    saturation_pressures = []
+
+    # Calculate densities and pressures at each temperature
+    for T in temperatures:
+        ld = CP.PropsSI('D', 'T', T, 'Q', 0, substance)         # liquid density
+        vd = CP.PropsSI('D', 'T', T, 'Q', 1, substance)         # vapor density
+        psat = CP.PropsSI('P', 'T', T, 'Q', 0, substance) / 1e5 # saturation pressure
+
+        liquid_density.append(ld)
+        vapor_density.append(vd)
+        saturation_pressures.append(psat)
+
+    return liquid_density, vapor_density, saturation_pressures
+
+def plot_curves(isoT=15, isoP=60):
+    isoT += 273.15
+
+    # Create figure and axes
+    fig, axs = plt.subplots(1, 2, figsize=(15, 6), sharey=True)
+
+    # Get saturation data
+    liquid_density, vapor_density, saturation_pressures = get_saturations(temperatures)
+
+    """Density vs Temperature"""
+    # Plot saturation curves with Temperature on x-axis
+    axs[0].plot(temperatures, liquid_density, label='Saturated Liquid Density')
+    axs[0].plot(temperatures, vapor_density, label='Saturated Vapor Density')
+    axs[0].fill_between(temperatures, liquid_density, vapor_density, color='gray', alpha=0.2)
+
+    # Add critical point
+    T_critical = CP.PropsSI('Tcrit', substance)
+    D_critical = CP.PropsSI('rhocrit', substance)
+    # print(f'T_critical = {T_critical:.2f} K')
+    # print(f'D_critical = {D_critical:.2f} kg/m^3')
+    axs[0].plot(T_critical, D_critical, 'kx', label='Critical Point')
+
+    # Plot isobar
+    isoP *= 1e5
+    isoP_Ts = np.linspace(260,370, 500)
+    try:
+        rho = CP.PropsSI('D', 'P', isoP, 'T', isoP_Ts, substance)
+        axs[0].plot(isoP_Ts, rho, label=f'Isobar at {(isoP / 1e5):.1f} bar')
+        T_sat = CP.PropsSI('T', 'P', isoP, 'Q', 0, substance)
+        axs[0].axvline(T_sat, color='gray', linestyle='--', label=f'$T_{{sat}}$ = {T_sat-273.15:.1f} °C')
+        # Plot saturation points
+        D_V = CP.PropsSI('D', 'P', isoP, 'Q', 1, substance)
+        D_L = CP.PropsSI('D', 'P', isoP, 'Q', 0, substance)
+        T_V = CP.PropsSI('T', 'P', isoP, 'Q', 1, substance)
+        T_L = CP.PropsSI('T', 'P', isoP, 'Q', 0, substance)
+        axs[0].plot(T_V, D_V, 'ro', label=f'$\\rho_{{V, sat}}$ = {D_V:.2f} kg/m$^3$')
+        axs[0].plot(T_L, D_L, 'bo', label=f'$\\rho_{{L, sat}}$ = {D_L:.2f} kg/m$^3$')
+
+        # # Annotate the density of the saturated liquid and vapor use LaTeX formatting and put into boxes
+        # axs[0].annotate(f'{D_L:.2f} kg/m$^3$', xy=(T_L, 900))
+        # axs[0].annotate(f'{D_V:.2f} kg/m$^3$', xy=(T_V, 0))
+    except:
+        pass
+
+    # Label plot
+
+    axs[0].set_title('Density vs Temperature')
+
+    # Example values for the limits and increments in Celsius
+    x_min_celsius = -10  # Minimum value in Celsius
+    x_max_celsius = 60  # Maximum value in Celsius
+    increments_celsius = 10  # Increment step in Celsius
+
+    # Convert these values to Kelvin for setting limits and ticks
+    x_min_kelvin = x_min_celsius + 273.15
+    x_max_kelvin = x_max_celsius + 273.15
+    increments_kelvin = increments_celsius
+
+    # Set the limits for the x-axis
+    axs[0].set_xlim(x_min_kelvin, x_max_kelvin)
+
+    # Set the ticks for the x-axis
+    # np.arange creates an array from x_min_kelvin to x_max_kelvin with a step of increments_kelvin
+    axs[0].set_xticks(np.arange(x_min_kelvin, x_max_kelvin + 1, increments_kelvin))
+
+    # Modify the x-axis to show temperature in degrees Celsius
+    axs[0].xaxis.set_major_formatter(FuncFormatter(lambda val, pos: f'{(val - 273.15):.0f}'))
+
+    # Set the label for the x-axis
+    axs[0].set_xlabel('Temperature (°C)')
+
+
+
+    axs[0].set_ylabel(f'Density kg/m$^3$')
+    axs[0].legend()
+
+
+    """Density vs Pressure"""
+    # Plot saturation curves with Pressure on x-axis
+    axs[1].plot(saturation_pressures, liquid_density, label='Saturated Liquid Density')
+    axs[1].plot(saturation_pressures, vapor_density, label='Saturated Vapor Density')
+    axs[1].fill_between(saturation_pressures, liquid_density, vapor_density, color='gray', alpha=0.2)
+
+    P_critical = CP.PropsSI('Pcrit', substance) / 1e5
+    # print(f'P_critical = {P_critical:.2f} bar')
+    axs[1].plot(P_critical, D_critical, 'kx', label='Critical Point')
+
+    # Plot isotherm
+    isoT_Ps = np.linspace(1e4, 90e5, 500)
+    try:
+        rho = CP.PropsSI('D', 'T', isoT, 'P', isoT_Ps, substance)
+        isoT_Ps_bar = [pressure / 1e5 for pressure in isoT_Ps]  # Convert to bar
+        axs[1].plot(isoT_Ps_bar, rho, label=f'Isotherm at {isoT-273.15:.1f} °C')
+        P_sat = CP.PropsSI('P', 'T', isoT, 'Q', 0, substance)
+        axs[1].axvline(P_sat / 1e5, color='gray', linestyle='--', label=f'$P_{{sat}}$ = {P_sat / 1e5:.2f} bar')
+
+        # Plot saturation points
+        D_V = CP.PropsSI('D', 'T', isoT, 'Q', 1, substance)
+        D_L = CP.PropsSI('D', 'T', isoT, 'Q', 0, substance)
+        P_V = CP.PropsSI('P', 'T', isoT, 'Q', 1, substance) / 1e5
+        P_L = CP.PropsSI('P', 'T', isoT, 'Q', 0, substance) / 1e5
+        axs[1].plot(P_V, D_V, 'ro', label=f'$\\rho_{{V, sat}}$ = {D_V:.2f} kg/m$^3$')
+        axs[1].plot(P_L, D_L, 'bo', label=f'$\\rho_{{L, sat}}$ = {D_L:.2f} kg/m$^3$')
+
+        # # Annotate the density of the saturated liquid and vapor use LaTeX formatting and put into boxes
+        # axs[1].annotate(f'{D_L:.2f} kg/m$^3$', xy=(P_L+2, D_L+50))
+        # axs[1].annotate(f'{D_V:.2f} kg/m$^3$', xy=(P_V + 3, 0))
+    except:
+        pass
+
+    axs[1].set_title('Density vs Pressure')
+    axs[1].set_xlabel('Pressure (bar)')
+    axs[1].set_xlim(0, 90)
+    axs[1].legend(loc = 'upper right')
+
+    # Display plot
+    plt.suptitle(f'$N_{2}O$ Density Curves', fontsize=20)
+    plt.tight_layout()
+    plt.show()
+
+def density_sliders():
+    isoT_slider = widgets.FloatSlider(
+        value=15,
+        min=-10,
+        max=60,
+        step=0.5,
+        description= r'$T_{\text{const}}$ (°C)',
+        continuous_update=False
+    )
+
+    isoP_slider = widgets.FloatSlider(
+        value=40,
+        min=1,
+        max=90,
+        step=.5,
+        description= r'$P_{\text{const}}$ (bar)',
+        continuous_update=False
+    )
+    return isoT_slider, isoP_slider
+
+if __name__ == '__main__':
+    plot_curves()