+1 (315) 557-6473 

Python Program to Implement Steady State Shock Equations Assignment Solution.


Instructions

Objective
Write a program to implement steady state shock equations in python.

Requirements and Specifications

program to implement steady state shock equations in python

Source Code

import numpy as np

import matplotlib.pyplot as plt

M = 6.0

gam = 5./3.

# ==========================================================================

# ==========================================================================

def ddx(x, F):

    Ns = F.size - 1

    dFdx = np.zeros(Ns + 1, dtype='double')

    dFdx[0] = (F[1] - F[0]) / (x[1] - x[0])

    for i in range(1, Ns):

        dFdx[i] = (F[i + 1] - F[i - 1]) / (x[i + 1] - x[i - 1])

    dFdx[Ns] = (F[Ns] - F[Ns - 1]) / (x[Ns] - x[Ns - 1])

    return dFdx

# ==========================================================================

# ==========================================================================

if __name__ == '__main__':

    electron = np.loadtxt('./electron_quantites.txt')

    ion = np.loadtxt('./ion_quantites.txt')

    x = electron[:,0] # positions

    rho = ion[:,1] # this is rho/rho0.

    ke = electron[:,2]

    Te = electron[:,1]

    Ti = ion[:,2]

    ki = ion[:,3]

    eta = ion[:,4]

    u = ion[:,5]

    p = ion[:,6]

    E = ion[:,7]

    # First, calculate du/dx

    dudx = ddx(x, u)

    # Calculate dTi/dx

    dTidx = ddx(x, Ti)

    # Calculate dTe/dx

    dTedx = ddx(x, Te)

    # From equation (1b), calculate p0u0^2 + p0

    eq1b = p + gam*M**2.*rho*u**2. - eta*dudx

    eq1b_right = gam*M**2 + p[0]

    # Equation (1c)

    eq1c = E*u + p*u - ki*dTidx - ke*dTedx - eta*u*dudx

    eq1c_right = E[0]/rho[0] + p[0]

    # Eq. 1a

    eq1a = rho*u

    #==========================================================================

    #==========================================================================

    # create plot

    fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(7,7))

    # Plot dTidx

    axes[0,0].plot(x, rho)

    axes[0,0].set_title(r'$\rho$')

    axes[0,0].set_xlabel('Position')

    axes[0,0].set_ylabel('Plasma mass density')

    axes[0,0].grid(True)

    axes[0,1].plot(x, ke, label = 'ke')

    axes[0,1].set_title(r'$k_{e}$')

    axes[0,1].set_xlabel('Position')

    axes[0,1].set_ylabel('Electron Conductivity')

    axes[0, 1].grid(True)

    axes[0,2].plot(x, Te, label = 'Te')

    axes[0,2].set_title(r'$T_{e}$')

    axes[0,2].set_xlabel('Position')

    axes[0,2].set_ylabel('Electron Temperature')

    axes[0, 2].grid(True)

    axes[1,0].plot(x, Ti, label = 'Ti')

    axes[1,0].set_title(r'$T_{i}$')

    axes[1,0].set_xlabel('Position')

    axes[1,0].set_ylabel('Ion Conductivity')

    axes[1, 0].grid(True)

    axes[1,1].plot(x, ki, label = 'ki')

    axes[1,1].set_title(r'$k_{i}$')

    axes[1, 1].set_xlabel('Position')

    axes[1, 1].set_ylabel('Ion conductivity')

    axes[1, 1].grid(True)

    axes[1,2].plot(x, eta, label = 'eta')

    axes[1,2].set_title(r'$\eta$')

    axes[1, 2].set_xlabel('Position')

    axes[1, 2].set_ylabel('Ion viscosity')

    axes[1, 2].grid(True)

    plt.show()

    fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(7,7))

    axes[0,0].plot(x, u, label = 'u')

    axes[0,0].set_title(r'$u$')

    axes[0, 0].set_xlabel('Position')

    axes[0, 0].set_ylabel('Bulk fluid velocity')

    axes[0, 0].grid(True)

    axes[0,1].plot(x, p, label = 'p')

    axes[0,1].set_title(r'$p$')

    axes[0, 1].set_xlabel('Position')

    axes[0, 1].set_ylabel('Plasma pressure')

    axes[0, 1].grid(True)

    axes[0,2].plot(x, E, label = 'E')

    axes[0,2].set_title(r'$E$')

    axes[0, 2].set_xlabel('Position')

    axes[0, 2].set_ylabel('Total Plasma Energy')

    axes[0, 2].grid(True)

    axes[1,0].plot(x, dudx, label = 'dudx')

    axes[1,0].set_title(r'$\frac{\partial u}{\partial x}$')

    axes[1, 0].set_xlabel('Position')

    axes[1, 0].set_ylabel('Rate of change of bulk velocity')

    axes[1, 0].grid(True)

    axes[1,1].plot(x, dTidx, label = 'dTidx')

    axes[1,1].set_title(r'$\frac{\partial T_{i}}{\partial x}$')

    axes[1, 1].set_xlabel('Position')

    axes[1, 1].set_ylabel('Rate of change of Ion Temp.')

    axes[1, 1].grid(True)

    axes[1,2].plot(x, dTedx, label = 'dTedx')

    axes[1,2].set_title(r'$\frac{\partial T_{e}}{\partial x}$')

    axes[1, 2].set_xlabel('Position')

    axes[1, 2].set_ylabel('Rate of change of Electron Temp.')

    axes[1, 2].grid(True)

    plt.show()

    fig = plt.figure(figsize=(7, 5), dpi=120)

    plot = plt.plot(x, (eq1b-eq1b_right)/eq1b_right*100., 'royalblue', linestyle='--', label='Eq. (1b) percent error')

    plot = plt.plot(x, (eq1c-eq1c_right)/eq1c_right*100., 'tomato', linestyle='-', label='Eq. (1c) percent error')

    plot = plt.plot(x, (eq1a - 1.0)/1.0 *100, 'purple', linestyle='-', label = 'Eq. (1a) percent error')

    plt.ylabel(r'Percent error [%]')

    plt.xlabel(r'Position $\hat{x}$')

    #plt.xlim([25.0, 32.5])

    plt.xlim([min(x), max(x)])

    legend = plt.legend(loc='best', shadow=False, fontsize='small')

    #plt.savefig('equation_error.eps', format='eps', dpi=1000)

    plt.grid(True)

    plt.show()

    plt.close()

    #==========================================================================

    #==========================================================================