noneq1dvibration.py

# Quasi-1D nozzle with HLL and nonequilibrium vibration of N2

# R.Deiterding Dec. 2016
# Modified by Matthew Beckett 2017

import numpy as np
import time
import math as m
import matplotlib
import matplotlib.pyplot as plt

print("start CPU time = {} s".format(time.clock()))

# Input data (or via function arguments)

L = 0.0294
P0 = 10*101325
theta_vib = 3340.0
C = 721.434e-6
K = 1.91e6 
R = 297
T0 = 4000
rho0 = P0/(R * T0)
gamma = 1.4
a0 = m.sqrt(gamma*R*T0)
nstep = 2000
dt = 0.0025


Nx = 121
U = np.zeros((Nx,4))
F, Ftil = np.zeros((Nx,4)), np.zeros((Nx,4))
J, J2, pA = np.zeros(Nx), np.zeros(Nx), np.zeros(Nx)



##############   NEW GEOMETRY  ##################
Nx = 121

dY = 1.0 / (Nx - 1.0)
XoverL = np.arange(-0.5, 0.5+dY, dY)

evin = 0.61
line = ((0.56-evin)/0.4)*(XoverL + 0.5) + evin
line2 = ((0.425-0.56)/(.1))*XoverL + 0.425
line3 = ((0.36-0.41)/(.5))*XoverL + 0.425



hit_point1 = 0
hit_point2 = 0
point1hit = False
for i in range(0, len(XoverL)):
    if float(XoverL[i]) > -0.1 and not point1hit:
        hit_point1 = i
        point1hit = True
    if float(XoverL[i]) > 1.0E-14:
        hit_point2 = i-1
        break

line[hit_point1:] = line2[hit_point1:]
line[hit_point2:] = line3[hit_point2:]

line = line / gamma

flip_point = 0
Y = XoverL + 0.15

A_ratio = 20*(0.15-0.33768675*(Y**(1/3))+0.10129268*(Y**(1/2))+0.26794919*Y)
A_ratio2 = 20*(0.09791353-Y)


for i in range(0, len(Y)):
    if Y[i] >= 0.0178699:
        flip_point = i
        break

A_ratio[0:flip_point] = A_ratio2[0:flip_point]

Afactor1 = (5.35898 - (2.25125/(Y**(2/3))) + (1.01293*(Y**(-1/2))))/A_ratio
Afactor2 = -20.*np.ones(Nx)/A_ratio

AFACTOR = Afactor2
AFACTOR[flip_point:] = Afactor1[flip_point:]

#################################################

# Initial condition

rho = np.zeros(Nx)
T = np.zeros(Nx)
V = np.zeros(Nx)


rho = 1-0.7*(XoverL + 0.5)
T = 1-0.5*(XoverL + 0.5)
V = (0.1+1.09*(XoverL+0.5))*np.sqrt(T)


z = 0

test = np.zeros(Nx)

for i in range(0, Nx):
    U[i][0] = rho[i]*A_ratio[i]
    U[i][1] = U[i][0]*V[i]
    U[i][2] = U[i][0]*((T[i]/((2./5.)*gamma)) + 0.5*V[i]**2 + line[i])
    U[i][3] = U[i][0] * line[i]
    



# Main loop

for n in range(1,nstep):

    # Predictor

    for i in range(0, Nx):
        T_tr = (2/5)*(gamma)*((U[i][2] / U[i][0]) - 0.5*((U[i][1] / U[i][0])**2) - (U[i][3] / U[i][0]))
        pA[i] = U[i][0]*(1/gamma)*T_tr
        evib_eq = ((R*theta_vib)/(m.exp(theta_vib / (T_tr*T0))-1))/(gamma*R*T0)
        tau_vib = ((C/((rho0*(a0**2))*pA[i]/A_ratio[i]))*m.exp((K/(T_tr*T0))**(1/3)))*(a0/L)
        #tau_vib = (a*m.exp(b*((T_tr*T0)**(1/3)))/(rho0*U[i][0]/A_ratio[i]))*(a0/L)
        F[i][0] = U[i][1]
        F[i][1] = U[i][1]**2/U[i][0] + pA[i]
        F[i][2] = (U[i][2]+pA[i])*U[i][1]/U[i][0]
        F[i][3] = U[i][3] * U[i][1] / U[i][0]
        J[i] = pA[i]*AFACTOR[i]
        J2[i] = U[i][0] * ((evib_eq - (U[i][3]/U[i][0])) / tau_vib)


    for i in range(0,Nx-1):
        Vl = U[i][1]/U[i][0]
        al = m.sqrt(gamma*pA[i]/U[i][0])
        Vr = U[i+1][1]/U[i+1][0]
        ar = m.sqrt(gamma*pA[i+1]/U[i+1][0])
        sl = min([Vl-al,Vr-ar])
        sr = max([Vl+al,Vr+ar])
        if sr < 0.:
            Ftil[i] = F[i+1]
        elif sl > 0.:
            Ftil[i] = F[i]
        else:
            Ftil[i] = (sr*F[i]-sl*F[i+1]+sl*sr*(U[i+1]-U[i]))/(sr-sl)

    for i in range(1, Nx-1):
        U[i] = U[i] + (-(Ftil[i] - Ftil[i-1])/dY)*dt
        U[i][1] = U[i][1] + J[i]*dt
        U[i][3] = U[i][3] + J2[i]*dt

    # Inflow condition (based on isentropic flow from stagnation conditions)

    Vin = U[1][1]/U[1][0]   # extrapolate velocity in inflow
    Tin = 1.0 - 0.5*(gamma - 1.0)*Vin**2
    rhoin = (1.0 + 0.5*(gamma - 1.0)*Vin**2/Tin)**(-1.0/(gamma - 1.0))
    #print("rhoin:")
    #print(rhoin)
    evib_in = ((R*theta_vib)/(m.exp(theta_vib / (T0*Tin))-1))/(gamma*R*T0)
    U[0][0] = rhoin*A_ratio[0]
    U[0][1] = U[0][0]*Vin
    U[0][2] = U[0][0]*((Tin/((2./5.)*gamma)) + 0.5*Vin**2 + evib_in)
    U[0][3] = U[0][0] * evib_in
    ##print(U[0])

    # Outflow condition (fixed pressure, extrapolate other quantitites)

    U[Nx-1] = U[Nx-2]


# Table output

derived_total, total_energy, internal_energy, kinetic_energy, rhop, Vp, Tp, pp, Mp, evibp, pressure = np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx), np.zeros(Nx)
for i in range(0, Nx):
    rhop[i] = U[i][0] / A_ratio[i]
    Vp[i] = U[i][1]/U[i][0]
    Tp[i] = (2./5.)*gamma * (U[i][2] / U[i][0] - (0.5 * Vp[i]*Vp[i]) - (U[i][3] / U[i][0]));
    Mp[i] = Vp[i]/m.sqrt(Tp[i])
    pp[i] = Tp[i]*rhop[i]
    evibp[i] = (U[i][3] / U[i][0])*(a0**2)/(R*T0)
    pressure[i] = ((pA[i]/A_ratio[i])*rho0*(a0**2))/P0
    total_energy[i] = U[i][2] / U[i][0]
    kinetic_energy[i] = 0.5*(Vp[i]*Vp[i])
    internal_energy[i] = (Tp[i]/((2./5.)*gamma))
    print ("%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f" % (XoverL[i], A_ratio[i], rhop[i], Vp[i], Tp[i], evibp[i], Mp[i]))

print("end CPU time = {} s".format(time.clock()))

derived_total = internal_energy + kinetic_energy + evibp

# Graphic output

matplotlib.rc('text', usetex = True)

params = {'text.latex.preamble' : [r'\usepackage{siunitx}', r'\usepackage{amsmath}', r'\usepackage{gensymb}']}
plt.rcParams.update(params)
fontProperties = {'family':'serif', 'weight': 'bold', 'size': 20}

fig = plt.figure(num=None, figsize=(8, 6), dpi=150, facecolor='w', edgecolor='k')
ax1 = fig.add_subplot(1, 1, 1)
ax1.set_xlabel(r"Distance along nozzle, $x/L$", fontProperties)
ax1.set_ylabel(r"Vibrational energy, $e_{vib}/RT_0$", fontProperties)
ax1.set_ylim(0.35, 0.70)
plt.plot(XoverL, evibp, color='black', lw=2)
plt.plot(XoverL, line*gamma, '--', color='black', lw=2)
ax1.set_xticklabels(ax1.get_xticks(), fontProperties)
ax1.set_yticklabels(ax1.get_yticks(), fontProperties)

plt.show()