ihmegroup / stanshock Goto Github PK

from StanShock.stanShock import stanShock
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
import time
import cantera as ct
from StanShock.utils import getPressureData


def main(mech_filename: str = "N2O2HeAr.xml",
         show_results: bool = True) -> None:

    # Tailored conditions for IIST shock tube
    
    # subscript 1 is for driven section initial conditions and subscript 4 is for driver section
    T1 = 300 
    T4 = T1
    p1 = 0.3566e5
    p4 = 27.3566e5

    tFinal = 12e-3 #time of simulations

    # plotting parameters
    fontsize = 12

    # provided geometry
    DDriven = 59.5e-3
    DDriver = DDriven
    LDriver = 3.2
    LDriven = 5.235

    # Set up gasses and determine the initial pressures
    u1 = 0.0
    u4 = 0.0  # initially 0 velocity
    gas1 = ct.Solution(mech_filename)
    gas4 = ct.Solution(mech_filename)
    gas1.TPX = T1,p1,"O2:0.20,N2:0.79,AR:0.01"
    gas4.TPX = T4,p4,"He:1"

    # set up geometry
    nX = 1500  # mesh resolution
    xLower = -LDriver
    xUpper = LDriven
    xShock = 0.0
    geometry = (nX, xLower, xUpper, xShock)
    DeltaD = DDriven - DDriver
    DeltaX = (xUpper - xLower) / float(nX) * 10  # diffuse area change for numerical stability

    def D(x):
        diameter = DDriven + (DeltaD / DeltaX) * (x - xShock)
        diameter[x < (xShock - DeltaX)] = DDriver
        diameter[x > xShock] = DDriven
        return diameter

    def dDdx(x):
        dDiameterdx = np.ones(len(x)) * (DeltaD / DeltaX)
        dDiameterdx[x < (xShock - DeltaX)] = 0.0
        dDiameterdx[x > xShock] = 0.0
        return dDiameterdx

    A = lambda x: np.pi / 4.0 * D(x) ** 2.0
    dAdx = lambda x: np.pi / 2.0 * D(x) * dDdx(x)
    dlnAdx = lambda x, t: dAdx(x) / A(x)

    # set up solver parameters
    print("Solving with boundary layer terms")
    boundaryConditions = ['reflecting', 'reflecting']
    state1 = (gas1, u1)
    state4 = (gas4, u4)

    ssbl = stanShock(gas1, initializeRiemannProblem=(state4, state1, geometry),
                     boundaryConditions=boundaryConditions,
                     cfl=0.5,
                     outputEvery=25,
                     includeBoundaryLayerTerms=True,
                     reacting = False,
                     diffusion = True,
                     DOuter=D,
                     Tw=T1,  # assume wall temperature is in thermal eq. with gas
                     dlnAdx=dlnAdx)
 
    ssbl.addProbe(5.117)  # pressure sensor 1
    ssbl.addProbe(5.012)  # pressure sensor 2
    ssbl.addProbe(4.907)  # pressure sensor 3
    ssbl.addProbe(4.532)  # pressure sensor 4
    
    ssbl.addXTDiagram("density")

    # Solve
    t0 = time.perf_counter()
    ssbl.advanceSimulation(tFinal)
    t1 = time.perf_counter()
    print("The process took ", t1 - t0)

    mpl.rcParams['font.size'] = fontsize
    
    fig1 = plt.figure("Figure 1")

    plt.plot(np.array(ssbl.probes[0].t) * 1000.0, np.array(ssbl.probes[0].p) / 1.0e5,label = 'pcb1')
    
    plt.plot(np.array(ssbl.probes[1].t) * 1000.0, np.array(ssbl.probes[1].p) / 1.0e5,label = 'pcb2')

    plt.plot(np.array(ssbl.probes[2].t) * 1000.0, np.array(ssbl.probes[2].p) / 1.0e5,label = 'pcb3')

    plt.plot(np.array(ssbl.probes[3].t) * 1000.0, np.array(ssbl.probes[3].p) / 1.0e5,label = 'pcb4')

    plt.title("PCB-pressure")
    plt.xlabel("Time (ms)")
    plt.legend(loc='upper left')
    plt.ylabel("Pressure (bar)")
    plt.grid()

    ssbl.plotXTDiagram(ssbl.XTDiagrams["density"],5)
    plt.title("density (kg/m3)")
    plt.tick_params(axis='y', which='both', labelleft='off', labelright='on')
    plt.yticks(range(int(tFinal*1000)+1))
    plt.grid(axis='y',linewidth=0.5, alpha=0.5)

    plt.show()

if __name__ == "__main__":
    main()