Python Program to Implement Steady State Shock Equations Assignment Solution

Instructions

Objective

Write a python homework to implement steady state shock equations.

Requirements and Specifications

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() #========================================================================== #==========================================================================