Not able to Simulate tailored interface conditions using viscous case (in inviscid simulation I am getting tailored conditions)

[darshan1664]

Hi, Professor @kgrogan87; and other prominent researchers @danmohad, @waitong94 , @eboigne. I am Darshan, I am doing research in shock tube systems at Indian Institute of Space Science and Technology (IIST), India. My project includes experiments in shock tube and tunnel. Before doing each experiment, I am simulating pre-decided conditions using StanShock (viscous mode), to get a rough idea of expected results, and it worked until now.

But now I have to run one tailored interface case in our shock tube. So the conditions (initial driver and driven pressures) for tailored interface I have calculated using standard equations (which are derived for ideal gas) available in literature (see eq. 2 and 3 in https://ntrs.nasa.gov/api/citations/19660020190/downloads/19660020190.pdf). So, as usual, I am doing simulation in Stanshock viscous model before experiment, but I am not getting tailored interface ( instead, I am getting interaction of reflected shock with contact surface results in expansion fan, also see the images attached for density XTdiagram); But when I tried inviscid simulation, I am getting tailored interface (see images attached for inviscid case). In viscous case, I also tried so many iterations of conditions which are near to the tailored conditions derived from theory, but none of them gave tailored interface, I am not able understand what is happening there. I have attached my case files for both viscous and invisicid case, please check them, suggest some action or shed some light on possible causes of this problem.

Here is the python input script I am using for viscous simulation of tailored condition for our shock tube. (I am using same script for inviscid case, but only difference is boolean value in includeBoundaryLayerTerms is False).

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

Here, first image is for viscous case (you can see expansion fan reflecting from contact surface interaction which is clearly not a tailored interface), and the second image for inviscid case (you can see that after interaction with reflected shock, contact surface is stationary and no waves emerges from this interaction).

[darshan1664]

Here is the input script and other related case files for tailored viscous and inviscid cases.
Inviscid_tailored.zip
viscous_tailored.zip

[corykinney]

@darshan1664 What you are seeing is not an issue with StanShock but a result of the fact that the tailoring equations do not take into effect viscous losses. As you noted, the inviscid simulation is properly tailored, which is expected since it matches the assumptions used to derive the tailoring equations. However, for the viscid simulation (and experimentally in shock tunnels), the shock attenuates as it travels down the tunnel; therefore, at the point where the reflected shock interacts with the contact surface, the conditions for the expanded driver gas (state 3) and the shocked driven gas (state 2) are no longer tailored.