Shock tube experiments are used to investigate nonequilibrium thermochemistry and radiative processes in reacting gas hypersonic flows. Boundary layer and shock structure are known to influence the spatial variation of the test slug state properties. This work derives and validates a novel method for viscous, quasionedimensional, nonequilibrium flow in a shock tube assuming constant shock speed. The proposed method fully resolves the shock structure. Mirels' estimate for boundary layer growth around the test slug determines the dilation rate at the centerline. This, along with relevant boundary conditions, appropriately models core flow in a shock tube. The flow equations are discretized by finite differences on a staggered grid. The resulting highly nonlinear set of algebraic equations is solved by Newton iterations. The Jacobian matrix is block tridiagonal with a Schur complement, allowing efficient inversion. This culminates in a unique and computationally efficient quasionedimensional method offering improved modeling of the physical characteristics of shock tube experiments. Results of a 3 km/s, 66.6 Pa argon test case solved by a viscous, axisymmetric Navier–Stokes solution had agreement with the proposed method in temperature and pressure profiles to within 2% and postshock velocity to within 15%. Reacting gas shock tube experiments in synthetic air and synthetic Titan atmospheres were analyzed. Radiance values in the nonequilibrium and equilibrium regions were compared under various assumptions for the shock structure and radial velocity distribution. These results highlight the necessity of a dedicated shock tube solver when analyzing shock tube thermochemistry, particularly when determining reaction rates and relaxation parameters.
I. INTRODUCTION
Nonequilibrium phenomena can significantly affect drag, heat flux, and communication processes^{1} of hypersonic vehicles. Vehicles designed to enter planetary atmospheres rely heavily on numerical calculations, which require accurate estimates of physical quantities such as reaction rate constants and intermolecular interaction parameters. The reaction rate constants and intermolecular parameters for transport quantities used in numerical predictions of hypersonic flows draw information from a broad range of sources, including experimental ground testing.^{2} However, the flow conditions typical of planetary entry are unattainable in continuous flow facilities due to extreme power requirements and are most commonly studied using shock tube facilities.
Shock tubes use a high sound speed driver gas to propagate a shock wave through the downstream tube filled with test gas. Experimenters typically match the freestream flight composition and density with the test gas fill condition, and match the shock speed to the flight vehicles speed.^{3} Various flow nonuniformities are present due to the complex processes occurring during these experiments, including effects from the driver operation,^{4} diaphragm rupture,^{5,6} and boundary layer effects.^{7} These processes determine the spatial and temporal variation of the test gas properties. As a result, an appropriate characterization of the test gas in a shock tube is not a trivial exercise.
A priori fullfacility simulations are unable to predict these nonuniformities due to experimental complexity and testtotest variation in both driver operation and diaphragm rupture processes. To date, no fullfacility simulations have been able to accurately match the experimentally measured shock speed profile nor the pressure histories down the tube.^{8–11} The use of these codes to trial or optimize thermochemical models is further compounded by the extremely large computational expense of sufficiently resolving the shock and boundary layer growth transiently as the shock propagates along the shock tube to the measurement location.^{8} Thus, the inherent computational expense required to replicate experimental conditions renders a priori fullfacility simulations unsuitable for quantitative analysis of shock tube experiments.
To allow experimental data to be matched a posteriori, simplifications to the flow are commonly introduced. As an example, rate coefficients for hightemperature reactions determined from spectroscopic data generated by shock tubes^{12–15} require assumptions about the temperature and pressure profiles of the test gas behind the shock. There are three a posteriori approaches to numerically model shock tube experiments in thermochemical nonequilibrium:

Postshock thermochemical relaxation from conditions found using the Rankine–Hugoniot equations for a normal shock

Stagnation line solutions of the flow over blunt bodies such spheres or cylinders

Lagrangian solvers that implement experimentally observed boundary conditions of shock trajectory and pressure measurements^{16,17}
The first approach uses the Rankine–Hugoniot equations to evaluate the chemically and vibrationally frozen postshock state of a normal shock, and then to allow thermochemical relaxation behind the shock.^{18–20} As the kinetic energy term dominates the energy equation, the postshock speed is effectively a constant value and therefore does not account for stagnation at the contact discontinuity. Another important limitation of this model is that the Rankine–Hugoniot jump condition is not realized due to diffusion, resulting in a misprediction of the reaction rates.
In the shock frame of reference, the postshock flow must stagnate against the contact discontinuity when the contact discontinuity and primary shock speeds are matched.^{21} Thus, the second approach is to use an analogy that the stagnation line of a shock tube is equivalent to the stagnation line flow into a blunt body such as a sphere.^{2} This approach has been used extensively^{22–25} as it allows for diffusive terms to be included in the numerical simulation. This allows the shock structure to be fully resolved (assuming mesh independence) and accounts for compression effects along the stagnation line.
Clarke et al.^{26} demonstrated that models based on the Rankine–Hugoniot relations and the blunt body stagnation line analogy are not equivalent to the stagnation of flow present in a shock tube. Therefore, they must undergo a spatial transformation to shock tube relevant coordinates when modeling reacting gas shock tube experiments. The two outlined approaches do not account for the rapid growth of the boundary layer behind the primary shock,^{21,27,28} which removes mass from the core flow (see Fig. 1). The associated postshock speed scales with the inverse of the square root of the postshock distance,^{28} compared with the linear decrease to stagnation associated with a blunt body and with the constant postshock velocity associated with a normal shock. This profile dictates the time of flight the fluid takes to move a certain distance behind the shock. For reacting flows, this determines how far thermochemical processes have proceeded and therefore the spatial profile of observed radiance for spectroscopy experiments. These effects have been found significant for modeling chemically reacting flows present in shock tube experiments, particularly for slower reacting flows.^{26,30,31}
Satchell et al. and Collen et al. demonstrated shock trajectory had a significant influence on the temperature profile of the test gas, and therefore on the prediction of radiance.^{16,20,32,33} Satchell et al. developed a solver (LASTA), which considers the influence of shock trajectory by a Lagrangian method coupled to an algebraic approximation for the pressure relaxation within the test slug.^{16} Each Lagrangian slice evolves at constant enthalpy and entropy under the action of sound waves propagating along the test slug. Mirels' boundary layer effects are also considered, further accounting for flow nonuniformities. This has been further developed by Steer et al. to include the spatial pressure profile at the end of a test as a boundary condition.^{17} This third approach utilizes the characteristic form of the Euler equations and consequently is unable to resolve viscous effects such as shock structure and thermal diffusion.
There is therefore a requirement for an analysis tool, which can appropriately account for the physical phenomena present in a shock tube in a computationally efficient manner. This does not currently exist in the literature, as a result the analysis of reaction rates or radiance using other models is flawed due to missing the physical behavior inherent in shock tube flows.
In this paper, we develop a novel quasionedimensional numerical method for simulating the core flow in shock tubes with a constant shock speed. This solution is the first to combine a viscous Navier–Stokes solution with appropriate shock tube boundary layer mass loss. Physically, this allows the state properties to be resolved through the shock layer, then convect behind the shock at the appropriate velocity with the relevant compression due to flow stagnation against the driver gas. These effects will be investigated through their influence on thermochemical kinetics by comparison of state properties. Comparison to experiments will be performed through emission spectroscopy measurements. The mass loss to the boundary layer is specified by an analytical Mirels' approach.^{21} Initially, a Park twotemperature model with reacting chemistry is incorporated into the model. This results in a computationally efficient reacting gas shock tube solver, which can be utilized as a test bed for various thermochemistry models, having direct application to reaction rate optimization and providing uncertainty analyses for shock tube experiments.
II. NUMERICAL FORMULATION
Therefore, for a flow containing n_{s} species, with n_{T} = 2 energy modes and a viscous state vector $q=[x,T],u,P$, we have $ns+nT+3$ equations thus ensuring closure of the problem.
A. Shock tube flow
The flow in a shock tube can be considered as quasionedimensional by assuming the core flow only varies axially, with the radial source term determined by a Mirels' analysis.^{21} This allows a system of quasionedimensional equations to be derived in cylindrical coordinates, with z the axial coordinate down the tube and r the radial distance from the centerline. z = 0 is located at the stagnation point, which is located at the contact discontinuity in fully developed shock tube flow (see Fig. 2).
Additionally, the rmomentum equation can be replaced with Eq. (15). As there is an additional unknown α, the maximal Mirels' length δ_{m} must be specified. A definition of shock location is also prescribed to ensure the distance from the shock to the stagnation point matches the Mirel's length [see Eqs. (24) and (25)].
By utilizing a thermochemical database containing derivatives of all dependent variables with respect to the state variables, all components of F, g, and h can be analytically partially differentiated with respect to q, α, and $zshock$. This allows the exact Jacobian J to be found.
B. Discretization and solution method
The system of equations is discretized using upwind schemes for firstorder derivatives, and central difference for secondorder derivatives. Temperatures, pressures, densities, and radial velocity components are evaluated at cell centers, while axial velocities are evaluated at cell edges (i.e., at $k+12$) (see Fig. 3).
This culminates in a discretized system of equations for a given cell k, found in entirety in Appendix A. Note that by substituting Eqs. (23) and (25), blunt body stagnation line problems and normal shock problems can be solved. This is described in Appendix B.
The grid is recursively refined by determining the difference instate variables between cells. If the difference between cells is above a specified convergence criterion, cells are split and the steady solution to the new grid found. This process iterates until the convergence criterion is met across the domain.
The solution approach utilized is a novel approach to solving a problem of this kind, where numerical derivatives are typically used resulting in slower and less stable solution convergence. Furthermore, the exact Jacobian found in NESS offers direct evaluation of the sensitivity of state variables to perturbations in any of the model inputs. Therefore, NESS offers improved convergence performance for optimization problems, with a direct evaluation of the solution's sensitivity to model inputs. The block matrix solution approach utilized by NESS allows additional constraints to be simply added, for example, allowing coupling with detailed models at the stagnation point,^{35} or adding more complex coupling to a boundary layer solution.
III. VALIDATION
The ability of the proposed method to evaluate thermochemical model performance is predicated on three key aspects:

Resolve state properties through a viscous shock layer

Appropriately model shock tube fluid mechanics, primarily

Modeling the postshock velocity profile

Stagnation of the flow against the driver gas


Appropriate implementation of the nonequilibrium thermochemistry model
Therefore, the validation activities undertaken aim to validate each of these individual aspects, giving confidence in the results produced by NESS. All problems were solved within two minutes using a single core of an Intel i71165G7 processor, with inputs of the fill conditions (pressure, temperature, and mole fractions of chemical species considered), shock speed, and shock tube radius. These inputs allow the maximal Mirels' length to be found and subsequently used by the solver. To facilitate the validation activities, NESS was used to evaluate relaxing Rankine–Hugoniot shock problems (denoted as RH), blunt body stagnation line problems, and the shock tube problem (denoted ST) (see Appendix B). There are five validation activities considered as follows:

Evaluation of shock structure in a constant Prandtl number flow

Two test cases related to the correct implementation of a reacting thermochemistry library.

Postshock relaxation of a singletemperature model with reacting chemistry.

Stagnation line evaluation of a cylinder in a nitrogen flow using a twotemperature, nonreacting gas model


A nonreacting argon shock tube test

A reacting nitrogen shock tube test.
A. Shock structure
The structure of a shock is dependent on the viscous effects present in the flow, hence why NESS must contain the viscous component of the Navier–Stokes equation to appropriately resolve this region. The greatest change instate variable properties occurs across the shock and therefore is also the limiting region for grid convergence. Morduchow and Libby evaluated the analytical solution to the shock structure resulting from a Mach 2 shock in air with a Prandtl number of 0.75.^{36} The velocity and pressure profiles produced by NESS for the same flow were compared against the analytical ideal gas solution of Morduchow and Libby to validate the correct implementation of the viscous Navier–Stokes equations in NESS (see Fig. 4).
The same problem was evaluated using an increasingly demanding convergence grid refinement criterion, with the discrepancy between the analytical and numerical results displayed in Fig. 5. The discrepancy in postshock velocity can be attributed to Morduchow and Libby using an ideal gas model, whereas the thermodynamic properties of the gas were evaluated using the JANAF polynomial fits.^{37} The jagged edges in the error profile correspond to changes in spacing between the grid points, and this phenomenon is not visible once the grid is sufficiently resolved.
The results in Fig. 5(b) can be interrogated further by comparing the value for the maximum discrepancy, shown as a logarithmic plot in Fig. 6. Firstorder numerical convergence is achieved, with a constant 0.15% error to the analytical result once convergence occurs. This corresponds to using a convergence criterion resulting in n > 5563.
B. Thermochemistry
NESS has been developed to test a variety of thermochemistry models and their performance compared against shock tube experiments. In this work, the solver uses an inhouse thermochemistry library developed at the University of Oxford (OCEAN). The library provides functions to evaluate thermodynamic properties, transport properties, reaction, and relaxation rates in nonequilibrium twotemperature plasmas. Each function evaluates both the quantities of interest and their exact Jacobians with respect to the flow variables. Two initial problems were considered in the validation, one to evaluate the reacting chemistry model, and the second to validate the translational–vibrational relaxation time model.
1. Chemistry
The reacting chemistry gas model was validated using a comparison case developed by Gollan,^{38} comparing the kinetic schemes of Marrone et al.^{39} and Gupta et al.^{40} This was implemented in the University of Queensland's POSHAX3 algorithm,^{19} which is a onedimensional Rankine–Hugoniot normal shock solver, which applies a chemical relaxation scheme behind a frozen postshock state. The flow considered is a Mach 12.28 shock in synthetic air at 300 K and 133.32 Pa, using a singletemperature model. The singletemperature model was implemented in NESS by specifying that T = T_{V} instead of the vibrational energy equation. In order to make results comparable to the results of Gollan,^{38} the radial velocity is set to zero and a fixed domain length is enforced. The comparisons between NESS and POSHAX3 are shown in Fig. 7, with NESS comparable with the POSHAX3 data to within 4% for the major species by 0.1 mm postshock. This is expected, as the same kinetic reaction scheme from Gupta et al. is used,^{40} and however, differences are found closer to the shock due to the viscous effects modeled by NESS. The results in this section validate the reacting chemistry code developed alongside NESS.
2. Translational–vibrational relaxation
The translational–vibrational relaxation of the thermochemistry used by NESS was validated using a nonreacting flow of nitrogen past a one meter infinite cylinder originally considered by Giordano et al.^{41} A steady inflow of chemically inert, Mach 6.5 nitrogen at 500 Pa was evaluated with NESS using a twotemperature model, with the radial outflow equation in B1 and stagnation conditions at the body surface. The temperature profiles along the stagnation line found by using translational–vibrational relaxation time parameters from Blackman^{42} and Millikan and White^{43} are shown in Fig. 8 alongside the results from Giordano using the same parameters.^{41} The results using both sets of relaxation parameters are comparable; although they do indicate Giordano et al. used an insufficiently converged mesh at the shock, this is also visible in the work by Gollan.^{38}
C. Nonreacting shock tube test
An argon shock tube test was used to isolate the unique hydrodynamics present in a shock tube from any thermochemical effects, taking advantage of argon's high ionization temperature to ensure the flow is not reacting. This offers direct insight into the ability of NESS to account for postshock convection appropriately, in addition to the pressure and temperature rise from stagnation of the flow against the driver gas. For this purpose, the Duff canonical case^{27} of a 3 km/s shock in 66.6 Pa argon was used to evaluate the performance of NESS in modeling the fluid dynamics in a shock tube via comparison with a fullfacility simulation using the FROSST code.^{8} A 50mmradius shock tube was specified, with the centerline of the flow taken when the shock reached 8 m down the tube. The maximal Mirels' length of the problem was found to be 605 mm and was set as the stagnation length for the NESS simulation. The shock trajectory of the simulated problem is shown in Fig. 9, which shows a decreasing shock velocity from 3.1 km/s at 3 m after the diaphragm to 3 km/s at the finishing location at 8 m. The temperature profile comparison between NESS and FROSST is shown in Fig. 10, with a maximum temperature discrepancy in the postshock test slug of less than 2%. The superior grid resolution of NESS allows the shock structure to be properly resolved, with the temperature profile matching satisfactorily for the first 120 mm. The increase in discrepancy at greater postshock distances is explained by the influence of shock speed variation. Satchell et al.^{32} and Collen et al.^{20} found that if gas further behind the shock has been processed by a faster shock, it also must have increased entropy and hence temperature. The updated version of LASTA by Steer et al.^{17} was run using the shock speed profile in Fig. 9, producing the shock speed corrected temperature profile shown in Fig. 10. This highlights a significant limitation of NESS as a shock tube solver, due to the implicit assumption of constant shock speed in the steady solution of the reacting Navier–Stokes equations.
The pressure and velocity profiles of NESS also show better resolution of the shock compared to FROSST (see Figs. 11 and 12). Away from the shock, the two codes produce comparable profiles, again with the discrepancies attributable to the shock speed variance. The agreement in the postshock velocity profile is especially important for reacting nonequilibrium flows.^{26} Additionally, the compression of the test slug along the stagnation line increases the pressure from 7.1 to 8.5 kPa, this combined with the temperature rise results in a test slug prediction differing to that predicted by only using normal shock relationships. The NESS pressure profile is within 1% of the FROSST value, and the axial velocity in the shock frame of reference being within 15%.
As Mirels' analysis produces a singularity at the shock, a smoothing function was introduced to ensure the outflow increased over the shock thickness proportionally with the axial velocity profile. This increased to a cutoff value which is related to the curvature of the shock [see Fig. 13(a)]. This parameter should not be interpreted as a value of radial velocity attained within the flow, rather it is a representation of the total convection away from the stagnation line due to shock curvature. Using an analytical shock curvature profile from De Boer^{44} and normal shock relations, it was found the maximal radial velocity due to shock curvature was 760 m/s. This value was used as the value for maximal radial velocity. This enables stable convergence of the problem, without introducing instabilities from the mass loss to the boundary layer. Additionally, the validity of Mirels' assumption of boundary layer growth can be seen in Fig. 13(b). Mirels found the density–velocity product should be directly proportional to a nondimensional length given by $1\u2212(zshock\u2212z)/\delta m$.^{21} By extrapolating to zero, the maximal Mirels' length of FROSST is longer than that calculated for a 3 km/s shock due to the shock speed variance. However, due to mixing at the contact discontinuity [from 0 to 0.02 in Fig. 13(b)], the last portion of the test slug is disturbed resulting in divergence from Mirels' approximation.
The agreement in temperature, pressure, and velocity profiles between NESS and FROSST validate the fluid dynamics of NESS and its applicability to shock tube modeling, especially when shock speed variation is limited.
D. Reacting shock tube case
To this point, validation activities have focused on isolated features in NESS (Sec. III A for shock structure, Sec. III B for thermochemistry, and Sec. III C for the unique fluid mechanics present in shock tubes). Further confidence in the results from NESS can be gained by combining the effect of these phenomena. Given the uniqueness of NESS in the way it models shock tube flows, a direct comparison was not possible to reacting gas numerical results found in the literature. Instead, a simplified hydrodynamic problem in thermochemical nonequilibrium is considered, a 6.2 km/s shock in 133.3 Pa pure nitrogen, in a 100mmdiameter shock tube.^{45} This case was examined by Kim and Boyd^{46} using a variety of welldetailed temperature models and assuming relaxation from the conditions immediately post a chemically and vibrationally frozen normal shock using the Rankine–Hugoniot equations. The twotemperature results of Kim and Boyd are compared to NESSRH and made as both models assume relaxation from conditions immediately after a chemically and vibrationally frozen normal shock. The primary difference between the models is that NESS uses the preferential dissociation model implemented in Gnoffo et al.,^{34} thus differing in the way dissociation energy source terms and electronic–rotational energy transfer are evaluated. Additionally, NESS is a viscous solver, compared to the nonviscous solver used by Kim and Boyd. Therefore, there is sufficient agreement in the results of Fig. 14 to indicate the combined reacting gas and hydrodynamic implementation of NESS is valid, while also indicating the difference in implemented thermochemistry models.
IV. IMPACT OF SHOCK STRUCTURE AND BOUNDARY LAYER EFFECTS
The same test case used in Sec. III D is used in conjunction with a lowpressure condition to evaluate the importance of using a realistic representation for shock structure and wall boundary layer. By reducing the fill pressure of the test case, Mirels' length will be reduced^{21} and the distance for reactions to occur will be increased.^{26} The chosen shock speed is sufficiently high such that there exists dissociation and ionization processes within the flow. This produces a generic test case condition, which can be used to evaluate the applicability of the Rankine–Hugoniot model to shock tube modeling. Both problems use a 100mmdiameter shock tube filled with pure nitrogen and identical thermochemistry models. NESSST resolves the flow through the shock and accounts for mass loss to the boundary layer, while NESSRH accounts for neither. Therefore, the following analyses illustrate the relative importance of including hydrodynamic effects in simulations, and their impact on understanding the thermochemistry in a shock tube experiment.
A. 6.2 km/s shock through 133.3 Pa nitrogen
The first condition is the identical condition used by Kim and Boyd^{46} in Sec. III D, i.e., a 6.2 km/s shock through 133.3 Pa nitrogen gas in a 100mm shock tube. This condition has a maximal Mirels' length of 744 mm and an estimated maximum radial velocity of 774 m/s. Figure 15 demonstrates the effect of considering shock structure, as time of flight and compression effects are small for this region when postshock pressure is high, especially when the maximal Mirels' length is large.^{26} Due to the resolution of shock thickness found in NESSST, there is an appreciable rise distance for the translational–rotational temperature, resulting in an offset between the models during the decay toward equilibrium values. This also results in a decrease in peak translational temperature by approximately 4%. However, the same offset in vibroelectronic temperature is not observed, resulting in a spatially different combination of the two temperatures.
Furthermore, detail regarding the discrepancies between the solutions is displayed in Figs. 15(c) and 15(d). Immediately behind the shock, the vibroelectronic temperature rises faster in the Rankine–Hugoniot model, due to the translational–rotational temperature already being at the frozen value [see Fig. 15(c)]. This results in an immediate decrease in translational–rotational temperature. This also causes the reactions in NESSRH to initially proceed faster than in NESSST; however, approximately 0.5 mm postshock ionization is reduced and remains below the shock tube relevant model. The discrepancy in temperatures decreases to nearly 0% by 10 mm postshock, however, the discrepancy in molarity remains above 10% for N^{+} and e^{−}. These results highlight the importance of resolving the shock structure when modeling nonequilibrium thermochemistry, with the effects propagating significantly further forward than just the immediate vicinity of the shock.
B. 6.2 km/s shock through 13.3 Pa nitrogen
A lowpressure condition was defined to better evaluate the effect of the mass loss to the boundary layer in a shock tube experiment. The chosen condition is of a 6.2 km/s shock through 13.3 Pa nitrogen gas in a 100mm shock tube, producing a maximal Mirels' length of 73 mm. The results of the simulations are shown in Fig. 16, with the same behavior at the shock evident as the previous condition. This is especially evident for N_{2}^{+} and e^{−}, which has greater than a 200% difference in mole fraction in the first millimeter after the shock. The increased ionization present in NESSRH is due to the larger translational–rotational temperature, which increases the rate of the impact ionization reaction. However, the more significant effect for this condition is the timeofflight effect present due to boundary layer mass loss. This effect becomes apparent by 10 mm postshock and results in a 20% difference in temperature and up to 180% difference in mole fraction by 70 mm postshock [see Figs. 16(c) and 16(d)]. This is an expected phenomenon due to the differing postshock velocity profiles, and this effect was examined in detail by Clarke et al.^{26}
Clarke et al. developed a spatial transformation to transform normal shock results onto a coordinate system relevant for comparison to shock tube results. This methodology used a timeofflight (TOF) approach to equate the postshock location behind a normal shock to the equivalent location behind a shock in a shock tube. This transformation was applied to the results produced by NESSRH to place the Rankine–Hugoniot normal shock results onto a shock tube relevant coordinate system and is denoted as TOF (see Fig. 17). The results are significantly closer in magnitude after the transformation, and this is evident in the comparison of the discrepancy between the TOF transformed results to the viscous shock tube model (NESSST) in Figs. 17(c) and 17(d). The application of the transformation becomes appreciable after approximately 10 mm, with temperature difference being less than 10% and mole fraction within 30%. The difference in temperature decays to an approximately 3% lower value than in NESSST. This is due to the compression effects present in NESSST, resulting in an increase in temperature and pressure behind the shock due to the stagnation of the flow. The effect of the shock structure also becomes more apparent after the transformation is made, as although the vibroelectronic temperature is within 2% after 10 mm, the translational–rotational temperature takes approximately 35 mm to reach less than 5% discrepancy between models. Similarly, the difference in molarity is significantly improved after the transformation; however, the differences due to compression effects and shock structure result in significant differences in the prediction of ions and also nonnegligible differences in the level of dissociation of N_{2}.
These results highlight the possibility of hiding important features of the flow when modeling shock tube experiments with a nonspecialized shock tube solver such as a Rankine–Hugoniotbased solver. Importantly, the thermochemical reaction rates will be incorrectly optimized when matching experimental data using a solver with a flawed hydrodynamic model. When these optimized thermochemical rates are then applied for a vehicle simulation, inaccurate predictions of the vehicle's aerothermodynamics arise. Such reaction rate models using a Rankine–Hugoniotbased solver are the commonly used Park models,^{47,48} while reaction rate models using a blunt body stagnation line solver include the Cruden model^{49} and the Johnston model.^{50} Therefore, when modeling reacting gas shock tube experiments, it is imperative that the shock structure and boundary layers are both appropriately modeled. This, along with further understanding of the twodimensional effects, allows the thermochemistry problem to be appropriately decoupled from the hydrodynamics present in shock tube experiments.
V. COMPARISON TO EXPERIMENTS
The conditions considered to this point have been generic, allowing direct comparison between numerical models to highlight their limitations in simulating the flow phenomena found in shock tubes. Two additional experimental conditions are considered to highlight these discrepancies. The first analyzes a condition relevant for Titan entry with the experiment performed in the NASAEAST facility. The ES6119 test condition is taken from Brandis and Cruden and is a 6.1 km/s shock traveling through 13.3 Pa, 2% CH_{4}, and 98% N_{2} by volume fill gas.^{51} The carbon present in this condition evolves spatially at a slower rate compared to synthetic air, and the low pressure further slows progression to thermodynamic equilibrium. This condition will be analyzed using NESS in shock tube mode, and compared against the original blunt body model used by Brandis and Cruden.^{51} The second condition examines an experiment relevant for low Earth orbit return, performed in the University of Oxford T6 Stalker Tunnel in aluminum shock tube (AST) mode. The T6s491 test condition is a 7.25 km/s shock traveling through 18 Pa, 79.1% N_{2}, and 20.9% O_{2} by volume fill gas.^{52} This condition will be analyzed using NESS in shock tube mode, using Rankine–Hugoniot relaxation and using a blunt body sphere (BB) model. Experimental radiance profiles from both conditions will be compared to modeled radiance profiles using flow solutions coupled with the radiance solver NEQAIR (see Sec. V A).
A. Radiance model
NASA's NEQAIR v15.2^{53} is a line by line code developed to estimate radiance emissions, given a specified profile of temperatures and number densities. This allows an estimate of predicted radiance to be made by coupling a shock tube solver's output of spatially resolved temperature and number density profiles with the NEQAIR program. The model used by the analysis assumed a flux limited nonBoltzmann distribution with local escape factor of 1.0. This corresponds to a line 3 input of N F L 1.0 in the neqair.inp file. Experimentally determined spatial resolution functions (SRF) and instrument line shapes (ILS) convolve the NEQAIR simulations to account for broadening mechanisms contained within the experimental setup, thus enabling direct comparison from experimental to numerical results.
B. NASAEAST experiment for titan entry
The first test case is taken from NASAEAST test series 61, which focused on Titan entry and was examined in the work of Brandis and Cruden^{51} and reanalyzed by Clarke et al. to emphasize the importance of appropriately accounting for boundary layer mass loss in shock tubes.^{26} Reaction rates relevant for Titan entry are of interest due to the presence of CN and C_{2}, and their potential to cause significant increases in radiative heat flux.^{51} The experimental data are taken from shot T6119, a 6.1 km/s shock through a 101.6mmdiameter tube filled with 13.3 Pa, 300 K fill gas composed of 2% CH_{4}, and 98% N_{2} by volume. This test had a flat shock trajectory, minimizing the effects of shock history. The maximal test slug length was calculated to be 77 mm using a Mirels' estimate.^{21} NESSST was used to simulate the test condition and is compared to radiance profiles generated from the original DPLR^{54} simulation of a 3 m sphere. The experimental data used are from regions associated with strong CN radiance, from 370 to 440 nm (CN violet) and 480–680 nm (CN red). The original DPLR^{54} results of the stagnation line of a 3mradius sphere have been kindly provided by Brandis and Cruden. This allows direct evaluation of the radiance estimate produced by NEQAIR 15.2 for both the original DPLR simulation and the NESS simulation, shown in Fig. 18.
Radiance is overpredicted by approximately 50% by both models for the 370–440 nm region, and peak radiance underpredicted by 20% in the 440–680 nm region, indicating an inability of the utilized thermochemistry model to correctly characterize reaction rates. The predicted peak radiance values from the original DPLR simulations (DPLRNEQAIR) are 4% and 8% of the NESSSTNEQAIR values for the 370–440 and 480–680 nm regions, respectively. NESSSTNEQAIR provides the better match for the postshock decay profile for the two wavelength regions compared to the blunt body DPLRNEQAIR simulation, which is made more apparent when the radiance profiles are normalized by peak radiance (see Fig. 19). By 70 mm postshock distance, the predicted peak radiance from the original DPLR simulations is over 100% greater than the predicted NESSSTNEQAIR values. As shown in Sec. IV and in recent work by Clarke et al.,^{26} the time of flight for the shock tube flow will be vastly different to the blunt body sphere solution. This is due to the rapid growth of the boundary layer in the immediate postshock region, slowing the core flow down and subsequently increasing the time of flight of the flow.
This is further illustrated in Fig. 20, where the difference in temperature profiles and CN number densities highlight the discrepancy between models. Although the magnitude of peak number densities and peak temperatures are similar, the variance in spatial distributions result in peak radiance being delayed by approximately 3 mm for the DPLR model. A greater overshoot of the vibroelectronic temperature is observed in the DPLR simulation, likely due to differences in preferential dissociation modeling and grid resolution. The NESS simulation utilized approximately 15 000 points, compared to 800 points along the stagnation line of the DPLR simulation.
This difference in grid resolution is shown in Fig. 21, where the increase in grid mesh in NESS allows the shock structure to be resolved sufficiently with mesh independence. This is possible due to the numerical implementation, where the quasi1D nature of NESS coupled with solving the linearized system of equations with an analytical Jacobian provides the necessary computational efficiency to resolve the shock in thermochemical nonequilibrium. The difference in mesh refinement strategy is also clear, with the DPLR simulation using a monotonically decreasing function for grid point density compared to the targeted increase in grid point density used by NESS to resolve the shock and the resulting complex thermochemistry.
C. T6 low earth orbit return experiments
The second test case is shot T6s491, a 7.25 km/s shock with a relatively flat shock speed profile in 18 Pa synthetic air by volume, taken from an experimental campaign investigating the behavior of synthetic air and nitrogen in shocks ranging from 6 to 8 km/s.^{52} This campaign was conducted in the University of Oxford T6 Stalker Tunnel, operating in the aluminum shock tube (AST) mode.^{3,55} Of particular note is the large diameter of the tube, 225 mm at the viewing window. This produces a smaller boundary layer relative to the size of the core flow, resulting in large maximal Mirels' test slug length and producing radiance profiles with larger signal to noise ratios. Spatial resolution functions (SRF) and instrument line shapes (ILS) were determined by Glenn et al. for each test, allowing the numerical results to be convolved spectrally and spatially for direct comparison to experimental data. The maximal Mirels' test slug length was found to be 454 mm, which compares with an approximated standoff distance of 112 mm for a sphere with a 5 m radius.^{56} These values were utilized by NESS in shock tube mode (NESSST) and NESS in blunt body mode (NESSBB), respectively, with the 5mradius sphere chosen to ensure the estimated standoff distance is greater than the length of the shock tube viewing window.
The resulting radiance profiles are displayed in Fig. 22. The peak radiance values in both regions are overpredicted by all models; however, the peak value for NESSSTNEQAIR and NESSBBNEQAIR is within 2%. This is expected, given that both models only differ in the velocity loss in the postshock region. This minimizes the peak radiance difference in fast reacting flows. In comparison, there is a 10% difference from the NESSRHNEQAIR radiance values to the NESSSTNEQAIR value in the 585–850 nm region, due to the difference in shock structure modeling. Although the absolute magnitude of radiance decays toward a small value, the discrepancy in postshock velocity profiles between the models results in a 20% difference between NESSSTNEQAIR and the other models. One of the primary causes of these differences is visible in Fig. 23(a), where the lack of shock structure in the NESSRH solution produces a sharper decrease in transrotational temperature while the electrovibrational temperature follows a more similar profile to the NESSST and NESSBB results. This results in increased radiance in the 200–400 nm wavelength regions, which is associated with the molecular species, while the differing atomic oxygen number densities [see Fig. 23(b)] causes the greater discrepancy in peak radiance in the 585–850 nm region. This is due to the strong radiance produced by the oxygen 777 nm line, the lack of shock structure in the NESSRH simulation results in a different production rate of atomic oxygen [see inset of Fig. 23(b)]. Figure 23(b) also shows the difference in time of flight between the three models, with the postshock flow stretching the results of NESSRH and compressing the results of NESSBB due to the shock standoff distance being less than the maximal Mirels' length.
VI. CONCLUSION
The unique fluid dynamics present in a shock tube experiment due to boundary layer and shock structure results in a different flow environment compared to normal and blunt body shocks. A unique test bed for reacting gas thermochemistry models has been developed using a quasionedimensional stagnating flow solver with an outflow condition specified by a Mirels' approach. Computational efficiency is achieved through linearization of the system of equations and use of an analytical Jacobian, allowing complex systems containing many chemical species and reactions to be analyzed. This results in the unique capacity of resolving nonequilibrium flow across a shock such that the structure is mesh independent, while retaining a physical method of accounting for shock tube boundary layer effects. This solver is valid for constant shock speed problems, although it is important that any analysis of shock tube data has consideration of the influence of shock trajectory. Therefore, this solver is limited to analyzing shock tube experiments with flat shock speed profiles.
The solver was applied to a variety of test conditions, which directly highlight the importance of using a numerical model accounting for shock structure effects and boundary layer mass loss. Two experimental shock tube test conditions with radiance data were analyzed, with the improved modeling of physical phenomena resulting in a maximal peak radiance difference of 10%, and over 200% change during the decay toward thermochemical equilibrium. The differences between models were most significant for lowerpressure conditions and for slower reacting flows such as those relevant for Titan entry. The comparison to experimental results highlights the necessity of a dedicated shock tube solver for simulation of these flows, and the importance of analyzing and testing thermochemistry models with such a simulation tool. By combining simple approaches accounting for the physical characteristics of shock tube flow, coupled with an elegant and efficient solution method, a unique solver has been developed that can decouple gas thermochemistry from shock tube hydrodynamics. Thus, this solver can provide an improved understanding of the primary reaction mechanisms in highspeed flows relevant for entry, through understanding the sensitivity of simulated shock tube experiments to changes in the chosen thermochemistry model. Therefore, this solver can be directly applied to optimization of reaction rates through analysis of experimental datasets already in existence.
ACKNOWLEDGMENTS
Justin Clarke gratefully acknowledges the Rhodes Trust for funding and supporting his DPhil. Samuel Brody gratefully acknowledges the Marshall Commission for financial support as well as the U.S. Air Force for their support. The authors gratefully acknowledge Brett Cruden and Aaron Brandis for providing their DPLR stagnation line solution of the Titan test case.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
J. Clarke: Conceptualization (equal); Methodology (equal); Writing – original draft (equal). S. Brody: Methodology (supporting); Validation (equal); Writing – original draft (supporting). J. Steer: Validation (supporting); Writing – original draft (supporting). M. McGilvray: Supervision (equal); Writing – original draft (equal). L. Di Mare: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – original draft (supporting).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
NOMENCLATURE
 CD

Contact discontinuity
 D_{s}

Effective diffusion coefficient for species s ( $m2$/s)
 $D\u0302s$

Average vibrational energy per unit mass of species s removed or added to a gas mixture (J/kg)
 $ev,s$

Vibrational energy per unit mass of species s (J/kg)
 $ev,s*$

Vibrational energy per unit mass of species s evaluated at temperature T (J/kg)
 H

Total enthalpy per unit mass of mixture (J/kg)
 h

Mixture enthalpy per unit mass (J/kg)
 h_{s}

Enthalpy per unit mass of species s (J/kg)
 $hs,v$

Vibrational enthalpy per unit mass of species s (J/kg)
 $I\u0302s$

First ionization energy of species s per kgmole (J/kgmol)
 J

Diffusive mass flux (kg/ $m2$/s)
 $JJ$

Jacobian
 $Js$

Species diffusive mass flux (kg/ $m2$/s)
 M_{s}

Molecular weight of species s (kg/kgmol)
 n_{s}

Number of species
 $n\u0307ei$

Molar rate of production per unit volume of species i due to electron impact ionization (kgmol/ $m3$/s)
 P

Pressure (Pa)
 P_{e}

Electron pressure (Pa)
 Q

Collection of reacting gas source and relaxation terms
 q

Vector of viscous state variables
 $qc$

Conductive heat flux (W/ $m2$)
 $qd$

Diffusive enthalpy flux (J/ $m2$/s)
 R

Shock tube radius (m)
 r

Radial coordinate in the cylindrical coordinate system (m)
 $R\xaf$

Specific gas constant of the mixture (J/K/mol)
 T

Translational–rotational temperature (K)
 t

Time (s)
 T_{av}

Twotemperature model average temperature (K)
 T_{v}

Electrovibrational temperature (K)
 t_{flight}

Particle time of flight (s)
 u

Velocity vector (m/s)
 u

u_{r}Radial velocity component (m/s)
 u_{e}

Postshock velocity (m/s)
 u_{s}

Shock speed (m/s)
 v = u_{z}

Axial velocity component (m/s)
 $w\u0307i$

Mass rate production of species s [kg/ $(m)3$/s]
 $x\u2032$

Distance down the shock tube (m)
 y_{s}

Mole fraction of species s
 z

Distance from the contact discontinuity (m)
 $zshock$

Location of shock relative to the contact discontinuity (m)
 α

Constant used for Mirels' outflow (kg/ $m2$/s)
 δ_{b}

Boundary layer thickness (m)
 δ_{m}

Maximal Mirels' length (m)
 θ

Angular coordinate in the cylindrical coordinate system
 κ_{e}

Electronic thermal conductivity (W/m/K)
 κ_{T}

Translational–rotational thermal conductivity (W/m/K)
 κ_{v}

Vibrational thermal conductivity (W/m/K)
 μ

Mixture viscosity (Ns/ $m2$)
 ν_{es}

Effective collision frequency for electrons and species s (1/s)
 ρ

Density (kg/ $m3$)
 ρ_{e}

Partial density of electrons (kg/ $m3$)
 ρ_{s}

Partial density of the species s (kg/ $m3$)
 $\tau $

Stress tensor
 $\u27e8\tau \u27e9$

Total translational–vibrational energy relaxation time
 $\u27e8\tau i\u27e9$

Translational–vibrational energy relaxation time for species i
NOMENCLATURE
 CD

Contact discontinuity
 D_{s}

Effective diffusion coefficient for species s ( $m2$/s)
 $D\u0302s$

Average vibrational energy per unit mass of species s removed or added to a gas mixture (J/kg)
 $ev,s$

Vibrational energy per unit mass of species s (J/kg)
 $ev,s*$

Vibrational energy per unit mass of species s evaluated at temperature T (J/kg)
 H

Total enthalpy per unit mass of mixture (J/kg)
 h

Mixture enthalpy per unit mass (J/kg)
 h_{s}

Enthalpy per unit mass of species s (J/kg)
 $hs,v$

Vibrational enthalpy per unit mass of species s (J/kg)
 $I\u0302s$

First ionization energy of species s per kgmole (J/kgmol)
 J

Diffusive mass flux (kg/ $m2$/s)
 $JJ$

Jacobian
 $Js$

Species diffusive mass flux (kg/ $m2$/s)
 M_{s}

Molecular weight of species s (kg/kgmol)
 n_{s}

Number of species
 $n\u0307ei$

Molar rate of production per unit volume of species i due to electron impact ionization (kgmol/ $m3$/s)
 P

Pressure (Pa)
 P_{e}

Electron pressure (Pa)
 Q

Collection of reacting gas source and relaxation terms
 q

Vector of viscous state variables
 $qc$

Conductive heat flux (W/ $m2$)
 $qd$

Diffusive enthalpy flux (J/ $m2$/s)
 R

Shock tube radius (m)
 r

Radial coordinate in the cylindrical coordinate system (m)
 $R\xaf$

Specific gas constant of the mixture (J/K/mol)
 T

Translational–rotational temperature (K)
 t

Time (s)
 T_{av}

Twotemperature model average temperature (K)
 T_{v}

Electrovibrational temperature (K)
 t_{flight}

Particle time of flight (s)
 u

Velocity vector (m/s)
 u

u_{r}Radial velocity component (m/s)
 u_{e}

Postshock velocity (m/s)
 u_{s}

Shock speed (m/s)
 v = u_{z}

Axial velocity component (m/s)
 $w\u0307i$

Mass rate production of species s [kg/ $(m)3$/s]
 $x\u2032$

Distance down the shock tube (m)
 y_{s}

Mole fraction of species s
 z

Distance from the contact discontinuity (m)
 $zshock$

Location of shock relative to the contact discontinuity (m)
 α

Constant used for Mirels' outflow (kg/ $m2$/s)
 δ_{b}

Boundary layer thickness (m)
 δ_{m}

Maximal Mirels' length (m)
 θ

Angular coordinate in the cylindrical coordinate system
 κ_{e}

Electronic thermal conductivity (W/m/K)
 κ_{T}

Translational–rotational thermal conductivity (W/m/K)
 κ_{v}

Vibrational thermal conductivity (W/m/K)
 μ

Mixture viscosity (Ns/ $m2$)
 ν_{es}

Effective collision frequency for electrons and species s (1/s)
 ρ

Density (kg/ $m3$)
 ρ_{e}

Partial density of electrons (kg/ $m3$)
 ρ_{s}

Partial density of the species s (kg/ $m3$)
 $\tau $

Stress tensor
 $\u27e8\tau \u27e9$

Total translational–vibrational energy relaxation time
 $\u27e8\tau i\u27e9$

Translational–vibrational energy relaxation time for species i