Quasi-one-dimensional non-equilibrium method for shock tube and stagnation line flows

Clarke, J.; Brody, S.; Steer, J.; McGilvray, M.; Di Mare, L.

doi:10.1063/5.0218676

Shock tube experiments are used to investigate non-equilibrium 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, quasi-one-dimensional, non-equilibrium 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 non-linear 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 quasi-one-dimensional 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 post-shock velocity to within 15%. Reacting gas shock tube experiments in synthetic air and synthetic Titan atmospheres were analyzed. Radiance values in the non-equilibrium 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

Non-equilibrium phenomena can significantly affect drag, heat flux, and communication processes¹ 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 inter-molecular interaction parameters. The reaction rate constants and inter-molecular parameters for transport quantities used in numerical predictions of hypersonic flows draw information from a broad range of sources, including experimental ground testing.² 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.³ Various flow non-uniformities are present due to the complex processes occurring during these experiments, including effects from the driver operation,⁴ diaphragm rupture,^5,6 and boundary layer effects.⁷ 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 full-facility simulations are unable to predict these non-uniformities due to experimental complexity and test-to-test variation in both driver operation and diaphragm rupture processes. To date, no full-facility 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.⁸ Thus, the inherent computational expense required to replicate experimental conditions renders a priori full-facility 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 high-temperature 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 non-equilibrium:

Post-shock 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 post-shock 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 post-shock 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 post-shock flow must stagnate against the contact discontinuity when the contact discontinuity and primary shock speeds are matched.²¹ 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.² 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.²⁶ 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 post-shock speed scales with the inverse of the square root of the post-shock distance,²⁸ compared with the linear decrease to stagnation associated with a blunt body and with the constant post-shock 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

FIG. 1.

Diagram showing the flow processes in a shock tube when drawn in the shock frame of reference. Reproduced with permission from Clarke et al.,29 AIAA Paper No. 2023-1797, 2023. Copyright 2023 American Institute of Aeronautics and Astronautics, Inc.

View large Download slide

Diagram showing the flow processes in a shock tube when drawn in the shock frame of reference. Reproduced with permission from Clarke et al.,²⁹ AIAA Paper No. 2023-1797, 2023. Copyright 2023 American Institute of Aeronautics and Astronautics, Inc.

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.¹⁶ 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 non-uniformities. 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.¹⁷ 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 quasi-one-dimensional 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.²¹ Initially, a Park two-temperature 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

The non-equilibrium shock tube solver (NESS) solves the steady system of reacting gas Navier–Stokes equations in quasi-one-dimensional form. Consider the steady compressible and reacting Navier–Stokes equations:

species continuity \nabla \cdot (ρ_{s} u) = - \nabla \cdot J_{s} + {\dot{w}}_{s},

(1)

momentum \nabla \cdot (ρ u \otimes u) = - \nabla p + \nabla \cdot τ,

(2)

vibro - electronic \nabla \cdot (ρ e_{V} u) = - p_{e} \nabla \cdot u + \nabla \cdot ((κ_{v} + κ_{e}) \nabla T_{V} + \sum_{s = 1}^{n_{s}} h_{s, v} J_{s}) + Q,

(3)

total energy \nabla \cdot (ρ H u) = \nabla \cdot (τ \cdot u - q_{c} - q_{d}),

(4)

where Q represents source and relaxation terms defined by

Q = \sum_{s = mol .} ρ_{s} \frac{e_{v, s}^{*} - e_{v, s}}{⟨ τ_{s} ⟩} + 3 ρ_{e} \bar{R} (T - T_{V}) \sum_{s = 2}^{n_{s}} \frac{ν_{e s}}{M_{s}} + \sum_{s = mol .} {\dot{w}}_{s} {\hat{D}}_{s} - \sum_{s = ion .} {\dot{n}}_{e, s} {\hat{I}}_{s} .

(5)

The conductive heat flux is described by Fourier's law

q_{c} = - \sum_{i = 1}^{n_{T}} κ_{i} \nabla T_{i} .

(6)

The diffusive enthalpy flux in a reacting gas mixture is by Fick's law

q_{d} = - \sum_{s = 1}^{n_{s}} h_{s} J_{s} .

(7)

The diffusional mass flux is given by

J_{s} = - ρ D_{s} \nabla y_{s} .

(8)

The viscous stress tensor τ is defined as follows:

τ = μ (\nabla u + {(\nabla u)}^{T} - \frac{2}{3} (\nabla \cdot u)) .

(9)

The definitions of the variables contained in Eq. are found in the study by Gnoffo et al.³⁴ Finally, for reacting gases the mass conservation equations are extended by enforcing conservation of species

\frac{\partial ρ_{s}}{\partial t} + \nabla \cdot (ρ_{s} u) = - \nabla \cdot J_{s} + {\dot{w}}_{s} .

(10)

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 $n_{s} + n_{T} + 3$ equations thus ensuring closure of the problem.

A. Shock tube flow

The flow in a shock tube can be considered as quasi-one-dimensional by assuming the core flow only varies axially, with the radial source term determined by a Mirels' analysis.²¹ This allows a system of quasi-one-dimensional 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).

FIG. 2.

View large Download slide

Shock tube flow can be transformed from the lab frame of reference (a), into the shock frame of reference (b), and then idealized in a quasi-one-dimensional form with a radial velocity outflow condition (c).

This is formalized by making the following assumptions:

\frac{\partial}{\partial θ} = 0,

(11)

u_{θ} = 0,

(12)

u_{r} = u (z),

(13)

u_{z} = v (z),

(14)

ρ u_{r} \sqrt{z_{shock} - z} = α for z < z_{shock} .

(15)

With boundary conditions of

δ_{m} = z_{shock},

(16)

\frac{\partial q_{i}}{\partial z} = 0 at z = 0,

(17)

u_{z} = 0 at z = 0 .

(18)

Equation enforces a Mirels' boundary layer outflow condition, with α constant. This is represented by the arrows on the outside of the tube in Fig. 2(c), indicating the approximate radial velocity profile. This form is equivalent to Eq. (33) in Mirels;²¹ see Appendix C for the derivation. The radius of the core flow will approximately be the tube radius R, assuming a thin boundary layer such that

δ_{b} ≪ R

⁠.

Additionally, the r-momentum 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)].

This novel approach provides a simple yet elegant system of equations accounting for the effect of the boundary layer on the core flow, while retaining the necessary information required to shock and stagnate the flow in a viscous and reacting flow regime. Thus, the simplified system of equations describing the core flow of a shock tube in quasi-one-dimensionality are

species continuity 0 = - \frac{\partial (ρ_{s} v)}{\partial z} - \frac{2 ρ_{s} u}{R} + \frac{\partial}{\partial z} (ρ D_{s} \frac{\partial y_{s}}{\partial z}) + {\dot{w}}_{s},

(19)

z - momentum 0 = - \frac{\partial (ρ v^{2})}{\partial z} - \frac{2 ρ u v}{R} - \frac{\partial p}{\partial z} + \frac{4}{3} \frac{\partial}{\partial z} (μ \frac{\partial v}{\partial z}) - \frac{4}{3 R} \frac{\partial (μ u)}{\partial z} + \frac{2 μ}{R} \frac{\partial u}{\partial z},

(20)

\begin{matrix} vibro - electronic 0 = - \frac{\partial (ρ e_{V} v)}{\partial z} - \frac{2 ρ e_{V} u}{R} - p_{e} \frac{\partial v}{\partial z} - \frac{2 p_{e} u}{R} \\ + \frac{\partial}{\partial z} ((κ_{v} + κ_{e}) \frac{\partial T_{V}}{\partial z} - \sum_{s = 1}^{n_{s}} h_{s, v} ρ D_{s} \frac{\partial y_{s}}{\partial z}) + Q, \end{matrix}

(21)

\begin{matrix} total energy 0 = - \frac{\partial (ρ H v)}{\partial z} - \frac{2 ρ H u}{R} + \frac{2 μ v}{R} \frac{\partial u}{\partial z} - \frac{4 μ u}{3 R} \frac{\partial v}{\partial z} + \frac{\partial}{\partial z} (μ (\frac{4}{3} v \frac{\partial v}{\partial z} + u \frac{\partial u}{\partial z} - \frac{4}{3 R} u v)) \\ + \frac{\partial}{\partial z} (\sum_{i = 1}^{n_{T}} κ_{i} \frac{\partial T_{i}}{\partial z} - \sum_{s = 1}^{n_{s}} h_{s} ρ D_{s} \frac{\partial y_{s}}{\partial z}), \end{matrix}

(22)

Mirels equation 0 = ρ u \sqrt{z_{shock} - z} - α for z < z_{shock},

(23)

shock location 0 = g_{1} (q) = 0.5 P (0) - P (z_{shock}),

(24)

Mirels length 0 = g_{2} (z_{shock}) = δ_{m} - z_{shock} .

(25)

This system of equations can be simply described in vector form as

[\begin{matrix} F (q, α, z_{shock}) \\ g_{1} (q) \\ g_{2} (z_{shock}) \end{matrix}] = 0,

(26)

where q is the state vector, with the system completed by appending Eqs. and using the additional variables of

z_{shock}

and α. δ_m is used as a constant input, determined using the maximal Mirels' length.²¹ This form can be extended or adapted to enforce other relevant conditions by including additional equations.

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 $z_{shock}$ ⁠. This allows the exact Jacobian J to be found.

B. Discretization and solution method

The system of equations is discretized using upwind schemes for first-order derivatives, and central difference for second-order 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 + \frac{1}{2}$ ⁠) (see Fig. 3).

FIG. 3.

View large Download slide

Discretization stencil used by NESS, note cell-centered temperatures, pressures, and radial velocities, while axial velocity is defined at the cell interface.

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.

Collating Eqs. in the notation of Eq. , we can build a block matrix of the problem over the spatial domain

0 = [\begin{matrix} G (\tilde{q}, α, z_{k_{shock}}) \\ g_{1} (\tilde{q}) \\ g_{2} (z_{k_{shock}}) \end{matrix}] = [\begin{matrix} F_{1} (q^{1}, q^{2}, α, z_{k_{shock}}) \\ ⋮ \\ F_{k - 1} (q^{k - 2}, q^{k - 1}, q^{k}, α, z_{k_{shock}}) \\ F_{k} (q^{k - 1}, q^{k}, q^{k + 1}, α, z_{k_{shock}}) \\ F_{k + 1} (q^{k}, q^{k + 1}, q^{k + 2}, α, z_{k_{shock}}) \\ ⋮ \\ F_{N} (q^{N - 1}, q^{N}, α, z_{k_{shock}}) \\ g_{1} (\tilde{q}) \\ g_{2} (z_{k_{shock}}) \end{matrix}],

(27)

where

\tilde{q}

is the block vector of state variables.

\tilde{q} = [\begin{matrix} q^{1} \\ ⋮ \\ q^{k - 1} \\ q^{k} \\ q^{k + 1} \\ ⋮ \\ q^{N} \end{matrix}] .

(28)

Through calculation of the analytical derivatives, the full Jacobian is found and the problem can be solved using Newton–Raphson iterations

[\begin{matrix} {\tilde{q}}_{n + 1} \\ α_{n + 1} \\ z_{k_{shock}, n + 1} \end{matrix}] = [\begin{matrix} {\tilde{q}}_{n} \\ α_{n} \\ z_{k_{shock}, n} \end{matrix}] - J_{J}^{- 1} ({\tilde{q}}_{n}, α_{n}, z_{k_{shock}, n}) [\begin{matrix} {\tilde{q}}_{n} \\ α_{n} \\ z_{k_{shock}, n} \end{matrix}] .

(29)

To compute the inverse of the Jacobian J efficiently, we first break

J_{J}

into the block matrix form of

J_{J} = [\begin{matrix} \begin{matrix} \frac{\partial F_{1}}{\partial q_{1}} & \frac{\partial F_{1}}{\partial q_{2}} & 0 & \dots & 0 & \frac{\partial F_{1}}{\partial α} & \frac{\partial F_{1}}{\partial z_{k_{shock}}} \\ \frac{\partial F_{2}}{\partial q_{1}} & \frac{\partial F_{2}}{\partial q_{2}} & \frac{\partial F_{2}}{\partial q_{3}} & 0 & \dots & 0 & \frac{\partial F_{2}}{\partial α} & \frac{\partial F_{2}}{\partial z_{k_{shock}}} \\ ⋱ & ⋱ & ⋱ & ⋮ & ⋮ \\ 0 & \frac{\partial F_{k}}{\partial q_{k-1}} & \frac{\partial F_{k}}{\partial q_{k}} & \frac{\partial F_{k}}{\partial q_{k+1}} & \dots & 0 & \frac{\partial F_{k}}{\partial α} & \frac{\partial F_{k}}{\partial z_{k_{shock}}} \\ ⋱ & ⋱ & ⋱ & ⋮ & ⋮ \\ ⋮ & 0 & \frac{\partial F_{N-1}}{\partial q_{N-2}} & \frac{\partial F_{N-1}}{\partial q_{N-1}} & \frac{\partial F_{N-1}}{\partial q_{N}} & \frac{\partial F_{N-1}}{\partial α} & \frac{\partial F_{N-1}}{\partial z_{k_{shock}}} \\ 0 & \dots & 0 & \frac{\partial F_{N}}{\partial q_{N-1}} & \frac{\partial F_{N}}{\partial q_{N}} & \frac{\partial F_{N}}{\partial α} & \frac{\partial F_{N}}{\partial z_{k_{shock}}} \\ \frac{\partial g}{\partial q_{1}} & \dots & \frac{\partial g}{\partial q_{k}} & \dots & \frac{\partial g}{\partial q_{N-1}} & \frac{\partial g}{\partial q_{N}} & 0 & 0 \\ 0 & \dots & \dots & 0 & 0 & - 1 \end{matrix} \end{matrix}],

(30)

where

\frac{\partial F_{k}}{\partial q_{j}}

denotes the Jacobian of the vector function

F_{k}

⁠, associated with cell k with respect to the state variables of cell j. This can be rewritten in the following form for simplification:

J = [\begin{matrix} A & B \\ C & D \end{matrix}] .

(31)

The matrix A is block tridiagonal and can be efficiently inverted with a block Thomas algorithm. The Schur complement can then be taken as

M / A = D - C A^{- 1} B

allowing evaluation of the Jacobian inverse as

J^{- 1} = [\begin{matrix} A^{- 1} + A^{- 1} B {(M / A)}^{- 1} C A^{- 1} & - A^{- 1} B {(M / A)}^{- 1} \\ - {(M / A)}^{- 1} C A^{- 1} & {(M / A)}^{- 1} \end{matrix}] .

(32)

We also note the relative ease of computing the inverse of the 2 × 2 Schur complement M/A, thus reducing the total computational cost of calculating the inverse of the Jacobian.

The grid is recursively refined by determining the difference in-state 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,³⁵ 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 post-shock velocity profile
- Stagnation of the flow against the driver gas
Appropriate implementation of the non-equilibrium 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 i7-1165G7 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.
- Post-shock relaxation of a single-temperature model with reacting chemistry.
- Stagnation line evaluation of a cylinder in a nitrogen flow using a two-temperature, non-reacting gas model
A non-reacting 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 in-state 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.³⁶ 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).

FIG. 4.

Comparison of analytical results of the axial velocity (a) and pressure rise (b) evaluated using Eqs. (13) and (14) from Morduchow and Libby,36 and NESS.

View large Download slide

Comparison of analytical results of the axial velocity (a) and pressure rise (b) evaluated using Eqs. (13) and (14) from Morduchow and Libby,³⁶ and NESS.

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 post-shock 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.³⁷ 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.

FIG. 5.

Discrepancy between the analytical results for the axial velocity (a), and pressure rise (b) evaluated using Eqs. (13) and (14) from Morduchow and Libby,36 and NESS evaluated with an increasing number of grid points n.

View large Download slide

Discrepancy between the analytical results for the axial velocity (a), and pressure rise (b) evaluated using Eqs. (13) and (14) from Morduchow and Libby,³⁶ and NESS evaluated with an increasing number of grid points n.

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. First-order 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.

FIG. 6.

View large Download slide

Logarithm of the maximal discrepancy between the pressure profile of Morduchow and Libby and NESS with increased grid resolution.

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 in-house thermochemistry library developed at the University of Oxford (OCEAN). The library provides functions to evaluate thermodynamic properties, transport properties, reaction, and relaxation rates in non-equilibrium two-temperature 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,³⁸ comparing the kinetic schemes of Marrone et al.³⁹ and Gupta et al.⁴⁰ This was implemented in the University of Queensland's POSHAX3 algorithm,¹⁹ which is a one-dimensional Rankine–Hugoniot normal shock solver, which applies a chemical relaxation scheme behind a frozen post-shock state. The flow considered is a Mach 12.28 shock in synthetic air at 300 K and 133.32 Pa, using a single-temperature model. The single-temperature 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,³⁸ 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 post-shock. This is expected, as the same kinetic reaction scheme from Gupta et al. is used,⁴⁰ 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.

FIG. 7.

View large Download slide

Validation of NESS against POSHAX3 for a Mach 12.28 shock. (a) Temperature and density distribution and (b) mole concentrations.

2. Translational–vibrational relaxation

The translational–vibrational relaxation of the thermochemistry used by NESS was validated using a non-reacting flow of nitrogen past a one meter infinite cylinder originally considered by Giordano et al.⁴¹ A steady inflow of chemically inert, Mach 6.5 nitrogen at 500 Pa was evaluated with NESS using a two-temperature 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⁴² and Millikan and White⁴³ are shown in Fig. 8 alongside the results from Giordano using the same parameters.⁴¹ 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.³⁸

FIG. 8.

View large Download slide

Temperature profiles of the Mach 6.5, 500 Pa nitrogen stagnation line problem developed by Giordano. (a) Profile using Millikan and White constants and (b) profile using Blackman constants.

C. Non-reacting 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 post-shock 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²⁷ 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 full-facility simulation using the FROSST code.⁸ A 50-mm-radius 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 post-shock 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 post-shock distances is explained by the influence of shock speed variation. Satchell et al.³² and Collen et al.²⁰ 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.¹⁷ 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.

FIG. 9.

View large Download slide

Shock trajectory of the FROSST simulation.

FIG. 10.

View large Download slide

Temperature profile of the 3 km/s shock in 66.6 Pa argon.

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 post-shock velocity profile is especially important for reacting non-equilibrium flows.²⁶ 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%.

FIG. 11.

View large Download slide

Pressure profiles of the 3 km/s shock in 66.6 Pa argon.

FIG. 12.

View large Download slide

Axial velocity profiles of the 3 km/s shock in 66.6 Pa argon.

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⁴⁴ 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 non-dimensional length given by $1 - \sqrt{(z_{shock} - z) / δ_{m}}$ ⁠.²¹ 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.

FIG. 13.

View large Download slide

Radial velocity profile of NESS (a) and the comparative Mirels' analysis (b) of NESS with FROSST for a 3 km/s shock in 66.6 Pa argon.

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 non-equilibrium is considered, a 6.2 km/s shock in 133.3 Pa pure nitrogen, in a 100-mm-diameter shock tube.⁴⁵ This case was examined by Kim and Boyd⁴⁶ using a variety of well-detailed temperature models and assuming relaxation from the conditions immediately post a chemically and vibrationally frozen normal shock using the Rankine–Hugoniot equations. The two-temperature results of Kim and Boyd are compared to NESS-RH 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.,³⁴ 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 non-viscous 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.

FIG. 14.

$Comparison of temperatures (a) and mole fractions (b) using post-normal shock relaxation to the results of Kim and Boyd.46 6.2 km/s shock through 133.3 Pa nitrogen. Reproduced with permission from Kim and Boyd, AIAA Paper No. 2013-3150, 2013. Copyright 2013 American Institute of Aeronautics and Astronautics, Inc.$

View large Download slide

Comparison of temperatures (a) and mole fractions (b) using post-normal shock relaxation to the results of Kim and Boyd.⁴⁶ 6.2 km/s shock through 133.3 Pa nitrogen. Reproduced with permission from Kim and Boyd, AIAA Paper No. 2013-3150, 2013. Copyright 2013 American Institute of Aeronautics and Astronautics, Inc.

IV. IMPACT OF SHOCK STRUCTURE AND BOUNDARY LAYER EFFECTS

The same test case used in Sec. III D is used in conjunction with a low-pressure 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²¹ and the distance for reactions to occur will be increased.²⁶ 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 100-mm-diameter shock tube filled with pure nitrogen and identical thermochemistry models. NESS-ST resolves the flow through the shock and accounts for mass loss to the boundary layer, while NESS-RH 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⁴⁶ in Sec. III D, i.e., a 6.2 km/s shock through 133.3 Pa nitrogen gas in a 100-mm 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 post-shock pressure is high, especially when the maximal Mirels' length is large.²⁶ Due to the resolution of shock thickness found in NESS-ST, 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 vibro-electronic temperature is not observed, resulting in a spatially different combination of the two temperatures.

FIG. 15.

$Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 133.3 Pa nitrogen evaluated using a Rankine–Hugoniot solver (NESS-RH) and a viscous shock tube model (NESS-ST).$

View large Download slide

Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 133.3 Pa nitrogen evaluated using a Rankine–Hugoniot solver (NESS-RH) and a viscous shock tube model (NESS-ST).

Furthermore, detail regarding the discrepancies between the solutions is displayed in Figs. 15(c) and 15(d). Immediately behind the shock, the vibro-electronic 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 NESS-RH to initially proceed faster than in NESS-ST; however, approximately 0.5 mm post-shock ionization is reduced and remains below the shock tube relevant model. The discrepancy in temperatures decreases to nearly 0% by 10 mm post-shock, however, the discrepancy in molarity remains above 10% for N⁺ and e⁻. These results highlight the importance of resolving the shock structure when modeling non-equilibrium 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 low-pressure 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 100-mm 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₂⁺ and e⁻, which has greater than a 200% difference in mole fraction in the first millimeter after the shock. The increased ionization present in NESS-RH 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 time-of-flight effect present due to boundary layer mass loss. This effect becomes apparent by 10 mm post-shock and results in a 20% difference in temperature and up to 180% difference in mole fraction by 70 mm post-shock [see Figs. 16(c) and 16(d)]. This is an expected phenomenon due to the differing post-shock velocity profiles, and this effect was examined in detail by Clarke et al.²⁶

FIG. 16.

$Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 13.3 Pa nitrogen evaluated using a Rankine–Hugoniot solver (NESS-RH) and a viscous shock tube model (NESS-ST).$

View large Download slide

Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 13.3 Pa nitrogen evaluated using a Rankine–Hugoniot solver (NESS-RH) and a viscous shock tube model (NESS-ST).

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 time-of-flight (TOF) approach to equate the post-shock 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 NESS-RH 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 (NESS-ST) 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 NESS-ST. This is due to the compression effects present in NESS-ST, 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 vibro-electronic 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 non-negligible differences in the level of dissociation of N₂.

FIG. 17.

$Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 13.3 Pa nitrogen evaluated using TOF transformed results from the Rankine–Hugoniot solver (NESS-RH) and the viscous shock tube model (NESS-ST).$

View large Download slide

Temperature profiles (a) and mole fractions (b) with their respective differences (c) and (d) for the 6.2 km/s shock in 13.3 Pa nitrogen evaluated using TOF transformed results from the Rankine–Hugoniot solver (NESS-RH) and the viscous shock tube model (NESS-ST).

These results highlight the possibility of hiding important features of the flow when modeling shock tube experiments with a non-specialized shock tube solver such as a Rankine–Hugoniot-based 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–Hugoniot-based solver are the commonly used Park models,^47,48 while reaction rate models using a blunt body stagnation line solver include the Cruden model⁴⁹ and the Johnston model.⁵⁰ 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 two-dimensional 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 NASA-EAST facility. The ES61-19 test condition is taken from Brandis and Cruden and is a 6.1 km/s shock traveling through 13.3 Pa, 2% CH₄, and 98% N₂ by volume fill gas.⁵¹ 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.⁵¹ 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₂, and 20.9% O₂ by volume fill gas.⁵² 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⁵³ 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 non-Boltzmann 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. NASA-EAST experiment for titan entry

The first test case is taken from NASA-EAST test series 61, which focused on Titan entry and was examined in the work of Brandis and Cruden⁵¹ and reanalyzed by Clarke et al. to emphasize the importance of appropriately accounting for boundary layer mass loss in shock tubes.²⁶ Reaction rates relevant for Titan entry are of interest due to the presence of CN and C₂, and their potential to cause significant increases in radiative heat flux.⁵¹ The experimental data are taken from shot T61-19, a 6.1 km/s shock through a 101.6-mm-diameter tube filled with 13.3 Pa, 300 K fill gas composed of 2% CH₄, and 98% N₂ 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.²¹ NESS-ST was used to simulate the test condition and is compared to radiance profiles generated from the original DPLR⁵⁴ 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⁵⁴ results of the stagnation line of a 3-m-radius 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.

FIG. 18.

View large Download slide

Comparison of experimental radiance between (a) 370–440 and (b) 480–680 nm from NASA-EAST test T61-19 to NEQAIR estimates from NESS-ST and the original DPLR blunt body simulation. Experimental and blunt body results reproduced with permission from Brandis and Cruden,⁵¹ AIAA Paper No. 2017-4535, 2017. Not subject to copyright in the United States.

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 (DPLR-NEQAIR) are 4% and 8% of the NESS-ST-NEQAIR values for the 370–440 and 480–680 nm regions, respectively. NESS-ST-NEQAIR provides the better match for the post-shock decay profile for the two wavelength regions compared to the blunt body DPLR-NEQAIR simulation, which is made more apparent when the radiance profiles are normalized by peak radiance (see Fig. 19). By 70 mm post-shock distance, the predicted peak radiance from the original DPLR simulations is over 100% greater than the predicted NESS-ST-NEQAIR values. As shown in Sec. IV and in recent work by Clarke et al.,²⁶ 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 post-shock region, slowing the core flow down and subsequently increasing the time of flight of the flow.

FIG. 19.

View large Download slide

Comparison of experimental radiance between (a) 370–440 and (b) 480–680 nm from NASA-EAST test T61-19 to NEQAIR estimates from NESS-ST and the original DPLR blunt body simulation, normalized by peak radiance. Experimental and blunt body results reproduced with permission from Brandis and Cruden,⁵¹ AIAA Paper No. 2017-4535, 2017. Not subject to copyright in the United States.

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 vibro-electronic 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.

FIG. 20.

View large Download slide

Comparison of temperature (a) and CN number density (b) from simulations of NASA-EAST test T61-19 by NESS-ST and the original DPLR blunt body simulation. Blunt body results reproduced with permission from Brandis and Cruden,⁵¹ AIAA Paper No. 2017-4535, 2017. Not subject to copyright in the United States.

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 quasi-1D 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 non-equilibrium. 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.

FIG. 21.

View large Download slide

Comparison of grid resolution of NESS-ST and the original DPLR blunt body simulation of NASA-EAST test T61-19. Blunt body results reproduced with permission from Brandis and Cruden,⁵¹ AIAA Paper No. 2017-4535, 2017. Not subject to copyright in the United States.

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.⁵² 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 stand-off distance of 112 mm for a sphere with a 5 m radius.⁵⁶ These values were utilized by NESS in shock tube mode (NESS-ST) and NESS in blunt body mode (NESS-BB), respectively, with the 5-m-radius sphere chosen to ensure the estimated stand-off 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 NESS-ST-NEQAIR and NESS-BB-NEQAIR is within 2%. This is expected, given that both models only differ in the velocity loss in the post-shock region. This minimizes the peak radiance difference in fast reacting flows. In comparison, there is a 10% difference from the NESS-RH-NEQAIR radiance values to the NESS-ST-NEQAIR 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 post-shock velocity profiles between the models results in a 20% difference between NESS-ST-NEQAIR 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 NESS-RH solution produces a sharper decrease in transrotational temperature while the electro-vibrational temperature follows a more similar profile to the NESS-ST and NESS-BB 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 NESS-RH 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 post-shock flow stretching the results of NESS-RH and compressing the results of NESS-BB due to the shock stand-off distance being less than the maximal Mirels' length.

FIG. 22.

View large Download slide

Comparison of experimental radiance between 200 and 440 nm (a) and 585–850 nm (b) of T6s491 to NEQAIR estimates from NESS-ST, NESS-RH, and NESS-BB simulations.

FIG. 23.

View large Download slide

Comparison of temperatures (a) and atomic oxygen number densities (b) from simulations of T6s491 by NESS-ST, NESS-RH, and NESS-BB.

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 quasi-one-dimensional 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 non-equilibrium 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 lower-pressure 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 high-speed 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 (⁠ $m^{2}$ /s)
${\hat{D}}_{s}$: Average vibrational energy per unit mass of species s removed or added to a gas mixture (J/kg)
$e_{v, s}$: Vibrational energy per unit mass of species s (J/kg)
$e_{v, 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)
$h_{s, v}$: Vibrational enthalpy per unit mass of species s (J/kg)
${\hat{I}}_{s}$: First ionization energy of species s per kg-mole (J/kg-mol)
J: Diffusive mass flux (kg/ $m^{2}$ /s)
$J_{J}$: Jacobian
$J_{s}$: Species diffusive mass flux (kg/ $m^{2}$ /s)
M_s: Molecular weight of species s (kg/kg-mol)
n_s: Number of species
${\dot{n}}_{e i}$: Molar rate of production per unit volume of species i due to electron impact ionization (kg-mol/ $m^{3}$ /s)
P: Pressure (Pa)
P_e: Electron pressure (Pa)
Q: Collection of reacting gas source and relaxation terms
q: Vector of viscous state variables
$q_{c}$: Conductive heat flux (W/ $m^{2}$ ⁠)
$q_{d}$: Diffusive enthalpy flux (J/ $m^{2}$ /s)
R: Shock tube radius (m)
r: Radial coordinate in the cylindrical coordinate system (m)
$\bar{R}$: Specific gas constant of the mixture (J/K/mol)
T: Translational–rotational temperature (K)
t: Time (s)
T_av: Two-temperature model average temperature (K)
T_v: Electro-vibrational temperature (K)
t_flight: Particle time of flight (s)
u: Velocity vector (m/s)
u: u_rRadial velocity component (m/s)
u_e: Post-shock velocity (m/s)
u_s: Shock speed (m/s)
v = u_z: Axial velocity component (m/s)
${\dot{w}}_{i}$: Mass rate production of species s [kg/ ${(m)}^{3}$ /s]
$x'$: Distance down the shock tube (m)
y_s: Mole fraction of species s
z: Distance from the contact discontinuity (m)
$z_{shock}$: Location of shock relative to the contact discontinuity (m)
α: Constant used for Mirels' outflow (kg/ $m^{2}$ /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 (N-s/ $m^{2}$ ⁠)
ν_es: Effective collision frequency for electrons and species s (1/s)
ρ: Density (kg/ $m^{3}$ ⁠)
ρ_e: Partial density of electrons (kg/ $m^{3}$ ⁠)
ρ_s: Partial density of the species s (kg/ $m^{3}$ ⁠)
$τ$: Stress tensor
$⟨ τ ⟩$: Total translational–vibrational energy relaxation time
$⟨ τ_{i} ⟩$: Translational–vibrational energy relaxation time for species i

NOMENCLATURE

CD: Contact discontinuity
D_s: Effective diffusion coefficient for species s (⁠ $m^{2}$ /s)
${\hat{D}}_{s}$: Average vibrational energy per unit mass of species s removed or added to a gas mixture (J/kg)
$e_{v, s}$: Vibrational energy per unit mass of species s (J/kg)
$e_{v, 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)
$h_{s, v}$: Vibrational enthalpy per unit mass of species s (J/kg)
${\hat{I}}_{s}$: First ionization energy of species s per kg-mole (J/kg-mol)
J: Diffusive mass flux (kg/ $m^{2}$ /s)
$J_{J}$: Jacobian
$J_{s}$: Species diffusive mass flux (kg/ $m^{2}$ /s)
M_s: Molecular weight of species s (kg/kg-mol)
n_s: Number of species
${\dot{n}}_{e i}$: Molar rate of production per unit volume of species i due to electron impact ionization (kg-mol/ $m^{3}$ /s)
P: Pressure (Pa)
P_e: Electron pressure (Pa)
Q: Collection of reacting gas source and relaxation terms
q: Vector of viscous state variables
$q_{c}$: Conductive heat flux (W/ $m^{2}$ ⁠)
$q_{d}$: Diffusive enthalpy flux (J/ $m^{2}$ /s)
R: Shock tube radius (m)
r: Radial coordinate in the cylindrical coordinate system (m)
$\bar{R}$: Specific gas constant of the mixture (J/K/mol)
T: Translational–rotational temperature (K)
t: Time (s)
T_av: Two-temperature model average temperature (K)
T_v: Electro-vibrational temperature (K)
t_flight: Particle time of flight (s)
u: Velocity vector (m/s)
u: u_rRadial velocity component (m/s)
u_e: Post-shock velocity (m/s)
u_s: Shock speed (m/s)
v = u_z: Axial velocity component (m/s)
${\dot{w}}_{i}$: Mass rate production of species s [kg/ ${(m)}^{3}$ /s]
$x'$: Distance down the shock tube (m)
y_s: Mole fraction of species s
z: Distance from the contact discontinuity (m)
$z_{shock}$: Location of shock relative to the contact discontinuity (m)
α: Constant used for Mirels' outflow (kg/ $m^{2}$ /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 (N-s/ $m^{2}$ ⁠)
ν_es: Effective collision frequency for electrons and species s (1/s)
ρ: Density (kg/ $m^{3}$ ⁠)
ρ_e: Partial density of electrons (kg/ $m^{3}$ ⁠)
ρ_s: Partial density of the species s (kg/ $m^{3}$ ⁠)
$τ$: Stress tensor
$⟨ τ ⟩$: Total translational–vibrational energy relaxation time
$⟨ τ_{i} ⟩$: Translational–vibrational energy relaxation time for species i

APPENDIX A: DISCRETIZED SYSTEM OF EQUATIONS

By discretizing Eqs. , we obtain the following system of equations for a given cell k:

\begin{matrix} species continuity 0 = (- ρ_{s, k + 1} v_{k} + ρ_{s, k} v_{k - 1}) - \frac{2 ρ_{s, k} u_{k}}{R} (z_{k} - z_{k - 1}) + \frac{2 ρ_{k + \frac{1}{2}} D_{s, k + \frac{1}{2}} (y_{s, k + 1} - y_{s, k})}{z_{k + 1} - z_{k - 1}} - \frac{2 ρ_{k - \frac{1}{2}} D_{s, k - \frac{1}{2}} (y_{s, k} - y_{s, k - 1})}{z_{k} - z_{k - 2}} + {\dot{w}}_{s, k}, \end{matrix}

(A1)

\begin{matrix} z - momentum 0 = (- ρ_{k + 1} v_{k + \frac{1}{2}}^{2} + ρ_{k} v_{k - \frac{1}{2}}^{2}) - \frac{(ρ_{k} u_{k} + ρ_{k + 1} u_{k + 1}) v_{k} (z_{k + 1} - z_{k - 1})}{2 R} - P_{k + 1} + P_{k} + \frac{4}{3} \frac{μ_{k + 1} (v_{k + 1} - v_{k})}{z_{k + 1} - z_{k}} - \frac{4}{3} \frac{μ_{k} (v_{k} - v_{k - 1})}{z_{k} - z_{k - 1}} - \frac{4}{3} (μ_{k + 1} u_{k + 1} - μ_{k} u_{k}) + \frac{2 μ_{k}}{R} (u_{k + 1} - u_{k}), \end{matrix}

(A2)

\begin{matrix} vibro - electronic 0 = (- ρ_{k + 1} e_{k + 1} v_{k} + ρ_{k} e_{k} v_{k - 1}) - \frac{ρ_{k} e_{k} u_{k}}{R} (z_{k} - z_{k - 1}) - p_{e, k} (v_{k} - v_{k - 1}) - \frac{2 p_{e, k} u_{k}}{R} \\ + 2 ((κ_{v, k + \frac{1}{2}} + κ_{e, k + \frac{1}{2}}) \frac{T_{V, k + 1} - T_{V, k}}{z_{k + 1} - z_{k - 1}} - (κ_{v, k - \frac{1}{2}} + κ_{e, k - \frac{1}{2}}) \frac{T_{V, k} - T_{V, k - 1}}{z_{k} - z_{k - 2}}) \\ - \sum_{s = 1}^{n_{s}} [\frac{2 h_{s, v, k + \frac{1}{2}} ρ_{k + \frac{1}{2}} D_{s, k + \frac{1}{2}} (y_{s, k + 1} - y_{s, k})}{z_{k + 1} - z_{k - 1}} - \frac{2 h_{s, v, k - \frac{1}{2}} ρ_{k - \frac{1}{2}} D_{s, k - \frac{1}{2}} (y_{s, k} - y_{s, k - 1})}{z_{k} - z_{k - 2}}] + Q_{k}, \end{matrix}

(A3)

\begin{array}{l} total energy 0 = (- ρ_{k + 1} (h_{k + 1} + \frac{1}{2} v_{k}^{2}) v_{k} + ρ_{k} (h_{k} + \frac{1}{2} v_{k - 1}^{2}) v_{k - 1}) - 2 ρ_{k} (h_{k} + \frac{1}{2} v_{k - \frac{1}{2}}^{2}) u_{k} \frac{(z_{k} - z_{k - 1})}{R} + \frac{2 μ_{k} v_{k}}{R} (u_{k + 1} - u_{k}) \\ - \frac{2 μ_{k} u_{k}}{3 R} (v_{k + 1} - v_{k - 1}) + \frac{2 μ_{k + 1} u_{k + 1} (u_{k + 1} - u_{k})}{z_{k + 1} - z_{k - 1}} - \frac{2 μ_{k} u_{k} (u_{k} - u_{k - 1})}{z_{k} - z_{k - 2}} + \frac{4}{3} μ_{k + 1} v_{k} \frac{v_{k + 1} - v_{k}}{z_{k + 1} - z_{k}} - \frac{4}{3} μ_{k} v_{k - 1} \frac{v_{k} - v_{k - 1}}{z_{k} - z_{k - 1}} - \frac{4}{3 R} (μ_{k + 1} u_{k + 1} v_{k} - μ_{k} u_{k} v_{k - 1}) - \sum_{s = 1}^{n_{s}} [\frac{2 h_{s, k + \frac{1}{2}} ρ_{k + \frac{1}{2}} D_{s, k + \frac{1}{2}} (y_{s, k + 1} - y_{s, k})}{z_{k + 1} - z_{k - 1}} - \frac{2 h_{s, k - \frac{1}{2}} ρ_{k - \frac{1}{2}} D_{s, k - \frac{1}{2}} (y_{s, k} - y_{s, k - 1})}{z_{k} - z_{k - 2}}] + \sum_{i = 1}^{n_{T}} 2 (κ_{i, k + \frac{1}{2}} \frac{T_{i, k + 1} - T_{i, k}}{z_{k + 1} - z_{k - 1}} - κ_{i, k - \frac{1}{2}} \frac{T_{i, k} - T_{i, k - 1}}{z_{k} - z_{k - 2}}), \end{array}

(A4)

Mirels equation 0 = ρ_{k} u_{k} \sqrt{z_{k_{shock}} - z_{k + \frac{1}{2}}} - α for z < z_{k_{shock}},

(A5)

shock location 0 = g_{1} (\tilde{q}) = 0.5 P_{1} - P_{k_{shock}},

(A6)

Mirels length 0 = g_{2} (z_{k_{shock}}) = δ_{m} - z_{k_{shock}} .

(A7)

APPENDIX B: PROBLEMS DESCRIBED BY DIFFERING OUTFLOW MODELS

NESS can be adapted to evaluate a relaxing Rankine–Hugoniot shock problem, a normal shock problem, and a blunt body stagnation line problem in addition to the shock tube problem. By changing Eqs. and , a blunt body stagnation line problem can be approximately evaluated using a relationship for the radial velocity profile, which results in a linear decrease in the density–axial velocity product²⁶

blunt body equation 0 = ρ u - α for z < z_{shock},

(B1)

stand - off length 0 = Δ - z_{shock},

(B2)

where Δ is the shock stand-off distance. Alternatively, a normal shock problem can be evaluated by specifying that the radial velocity is zero everywhere in the domain, i.e.,

u (z) = 0, \forall z > 0

⁠, and removing the stagnation condition

v (0) = 0

⁠. It can be further adapted to the relaxing Rankine–Hugoniot shock problem by specifying the initial conditions as the immediate post-shock values given by the Rankine–Hugoniot equations, and assuming T_V = 300 K for the two-temperature model.

APPENDIX C: MIRELS RADIAL VELOCITY PROFILE

We start the derivation using Eq. (33) from Mirels.²¹ This specifies the axial density–velocity profile at a certain post-shock distance l,

ρ u_{z} = {(ρ u_{z})}_{0} (1 - \sqrt{\frac{l}{l_{m}}}) .

(C1)

Consider a portion of a shock tube with length

Δ z

and assume negligible boundary layer thickness compared to the tube diameter (see Fig. 24). Then, using a mass balance at length l₁, we have

ρ_{1} u_{z 1}

entering the cell face, and

ρ_{2} u_{z 2}

leaving the cell face at l₂. The difference in mass must be removed by

ρ u_{r}

over the length of the cell. This can be formalized as

(ρ_{2} u_{z 2} - ρ_{1} u_{z 1}) π r^{2} = 2 π r Δ z ρ u_{r} .

(C2)

FIG. 24.

View large Download slide

Portion of shock tube used for deriving Mirels' radial velocity profile.

Substituting in Eq. :

\frac{{(ρ u_{z})}_{0}}{\sqrt{l_{m}}} (\sqrt{l_{2}} - \sqrt{l_{1}}) π r^{2} = 2 π r Δ z ρ u_{r},

(C3)

ρ u_{r} = \frac{r {(ρ u_{z})}_{0} (\sqrt{l_{1} + Δ z} - \sqrt{l_{1}})}{2 Δ z \sqrt{l_{m}}} .

(C4)

Taking the limit as

Δ z

goes to zero

ρ u_{r} = lim_{Δ z \to 0} \frac{r {(ρ u_{z})}_{0} (\sqrt{l_{1} + Δ z} - \sqrt{l_{1}})}{2 Δ z \sqrt{l_{m}}},

(C5)

ρ u_{r} = lim_{Δ z \to 0} \frac{r {(ρ u_{z})}_{0} (\sqrt{l_{1} + Δ z} - \sqrt{l_{1}})}{2 Δ z \sqrt{l_{m}}} \frac{(\sqrt{l_{1} + Δ z} + \sqrt{l_{1}})}{(\sqrt{l_{1} + Δ z} + \sqrt{l_{1}})},

(C6)

ρ u_{r} = lim_{Δ z \to 0} \frac{r {(ρ u_{z})}_{0} (l_{1} + Δ z - l_{1})}{2 Δ z \sqrt{l_{m}} (\sqrt{l_{1} + Δ z} + \sqrt{l_{1}})},

(C7)

ρ u_{r} = lim_{Δ z \to 0} \frac{r {(ρ u_{z})}_{0}}{2 \sqrt{l_{m}} (\sqrt{l_{1} + Δ z} + \sqrt{l_{1}})}

(C8)

\Rightarrow ρ u_{r} = \frac{r {(ρ u_{z})}_{0}}{4 \sqrt{l_{m}} \sqrt{l}} \Rightarrow ρ u_{r} \sqrt{l} = α .

(C9)

Here, α is given by

α = \frac{r {(ρ u_{z})}_{0}}{4 \sqrt{l_{m}}} .

(C10)

REFERENCES

1.

S.

Ramjatan

,

A.

Lani

,

S.

Boccelli

,

B.

Van Hove

,

O.

Karatekin

,

T.

Magin

, and

J.

Thoemel

, “

Blackout analysis of Mars entry missions

,”

J. Fluid Mech.

904

,

A26

(

2020

).

https://doi.org/10.1017/jfm.2020.714

Google Scholar

Crossref

2.

C.

Park

,

Nonequilibrium Hypersonic Aerothermodynamics

(

Wiley

,

New York

,

1990

).

Google Scholar

3.

P.

Collen

,

L. J.

Doherty

,

S. D.

Subiah

,

T.

Sopek

,

I.

Jahn

,

D.

Gildfind

,

R.

Penty Geraets

,

R.

Gollan

,

C.

Hambidge

,

R.

Morgan

, and

M.

McGilvray

, “

Development and commissioning of the T6 Stalker tunnel

,”

Exp. Fluids

62

,

225

(

2021

).

https://doi.org/10.1007/s00348-021-03298-1

Google Scholar

Crossref

4.

M.

McGilvray

,

A. G.

Dann

, and

P. A.

Jacobs

, “

Modelling the complete operation of a free-piston shock tunnel for a low enthalpy condition

,”

Shock Waves

23

,

399

–

406

(

2013

).

https://doi.org/10.1007/s00193-013-0437-8

Google Scholar

Crossref

5.

C. J.

Simpson

,

T. R.

Chandler

, and

K. B.

Bridgman

, “

Effect on shock trajectory of the opening time of diaphragms in a shock tube

,”

Phys. Fluids

10

,

1894

–

1896

(

1967

).

https://doi.org/10.1063/1.1762384

Google Scholar

Crossref

6.

E. M.

Rothkopf

and

W.

Low

, “

Diaphragm opening process in shock tubes

,”

Phys. Fluids

17

,

1169

–

1173

(

1974

).

https://doi.org/10.1063/1.1694860

Google Scholar

Crossref

7.

H.

Mirels

, “

Attenuation in a shock tube due to unsteady-boundary-layer action

,”

Report No. NACA-TR-1333

,

1957

.

Google Scholar

8.

M.

Satchell

,

P.

Collen

,

M.

McGilvray

, and

L.

Di Mare

, “

Numerical simulation of shock tubes using shock tracking in an overset formulation

,” AIAA Paper No. 2020-2722,

2020

.

https://doi.org/10.2514/6.2020-2722

Google Scholar

9.

K.

Bensassi

and

A. M.

Brandis

, “

Time accurate simulation of nonequilibrium flow inside the NASA Ames electric arc shock tube

,” AIAA Paper No. 2019-0798,

2019

.

https://doi.org/10.2514/6.2019-0798

Google Scholar

10.

D.

Chandel

,

I.

Nompelis

, and

G. V.

Candler

, “

Computations of high enthalpy shock propagation in electric arc shock tube (EAST) at NASA Ames

,” AIAA Paper No. 2018-1722,

2018

.

Google Scholar

11.

D. V.

Kotov

,

H. C.

Yee

,

M.

Panesi

,

D. K.

Prabhu

, and

A. A.

Wray

, “

Computational challenges for simulations related to the NASA electric arc shock tube (EAST) experiments

,”

J. Comput. Phys.

269

,

215

–

233

(

2014

).

https://doi.org/10.1016/j.jcp.2014.03.021

Google Scholar

Crossref

12.

E.

Freedman

and

J. W.

Daiber

, “

Decomposition rate of nitric oxide between 3000 and 4300° K

,”

J. Chem. Phys.

34

,

1271

–

1278

(

1961

).

https://doi.org/10.1063/1.1731731

Google Scholar

Crossref

13.

M. G.

Dunn

and

J. A.

Lordi

, “

Measurement of O₂⁺ + e⁻ dissociative recombination in expanding oxygen flows

,”

AIAA J.

8

,

614

–

618

(

1970

).

https://doi.org/10.2514/3.5730

Google Scholar

Crossref

14.

S.

Byron

,

Interferometric Measurement in a Shock Tube of Dissociation Rates for Air and Its Component Gases

(

Cornell University

,

New York

,

1959

).

Google Scholar

15.

C.

Park

, “

Two-temperature interpretation of dissociation rate data for N₂ and O₂

,” AIAA Paper No. 88-0458,

1988

.

https://doi.org/10.2514/6.1988-458

Google Scholar

16.

M.

Satchell

,

M.

McGilvray

, and

L.

Di Mare

, “

Analytical method of evaluating nonuniformities in shock tube flows: Theory and development

,”

AIAA J.

60

,

654

–

615

(

2021

).

https://doi.org/10.2514/1.J060990

Google Scholar

Crossref

17.

J.

Steer

,

J.

Clarke

,

M.

McGilvray

, and

L.

Di Mare

, “

LASTA 2.0: Validation of a reverse time integration method

,” AIAA Paper No. 2024-0447,

2024

.

https://doi.org/10.2514/6.2024-0447

Google Scholar

18.

A. B.

Glenn

,

P. L.

Collen

, and

M.

McGilvray

, “

Experimental non-equilibrium radiation measurements for low-earth orbit return

,” AIAA Paper No. 2022-2154,

2022

.

https://doi.org/10.2514/6.2022-2154

Google Scholar

19.

D. F.

Potter

, “

Modelling of radiating shock layers for atmospheric entry at Earth and Mars

,” Ph.D. thesis (

University of Queensland

,

2011

).

Google Scholar

20.

P. L.

Collen

,

M.

Satchell

,

L.

Di Mare

, and

M.

McGilvray

, “

The influence of shock speed variation on radiation and thermochemistry experiments in shock tubes

,”

J. Fluid Mech.

948

,

A51

(

2022

).

https://doi.org/10.1017/jfm.2022.727

Google Scholar

Crossref

21.

H.

Mirels

, “

Test time in low-pressure shock tubes

,”

Phys. Fluids

6

,

1201

–

1214

(

1963

).

https://doi.org/10.1063/1.1706887

Google Scholar

Crossref

22.

B. A.

Cruden

,

A. M.

Brandis

, and

D. K.

Prabhu

, “

Compositional dependence of radiance in CO₂/N₂/Ar systems

,” AIAA Paper No. 2013-2502,

2013

.

https://doi.org/10.2514/6.2013-2502

Google Scholar

23.

A.

Brandis

,

C.

Johnston

,

B.

Cruden

, and

D.

Prabhu

, “

Investigation of nonequilibrium radiation for Mars entry

,” AIAA Paper No. 2013-1055,

2013

.

https://doi.org/10.2514/6.2013-1055

Google Scholar

24.

A. M.

Brandis

,

C. O.

Johnston

,

B. A.

Cruden

,

D.

Prabhu

, and

D.

Bose

, “

Uncertainty analysis and validation of radiation measurements for Earth reentry

,”

J. Thermophys. Heat Transfer

29

,

209

–

221

(

2015

).

https://doi.org/10.2514/1.T4000

Google Scholar

Crossref

25.

B. A.

Cruden

and

A. M.

Brandis

, “

Measurement of radiative nonequilibrium for air shocks between 7 and 9 km/s

,”

J. Thermophys. Heat Transfer

34

,

154

–

180

(

2020

).

https://doi.org/10.2514/1.T5735

Google Scholar

Crossref

26.

J.

Clarke

,

L.

Di Mare

, and

M.

McGilvray

, “

Spatial transformations for reacting gas shock tube experiments

,”

AIAA J.

61

,

3365

(

2023

).

https://doi.org/10.2514/1.J062604

Google Scholar

Crossref

27.

R. E.

Duff

, “

Shock-tube performance at low initial pressure

,”

Phys. Fluids

2

,

207

–

216

(

1959

).

https://doi.org/10.1063/1.1705910

Google Scholar

Crossref

28.

H.

Mirels

, “

Flow nonuniformity in shock tubes operating at maximum test times

,”

Phys. Fluids

9

,

1907

–

1912

(

1966

).

https://doi.org/10.1063/1.1761542

Google Scholar

Crossref

29.

J.

Clarke

,

P. L.

Collen

,

M.

McGilvray

, and

L.

di Mare

, “

Numerical simulation of a shock tube in thermochemical non-equilibrium

,” AIAA Paper No. 2023-1797,

2023

.

https://doi.org/10.2514/6.2023-1797

Google Scholar

30.

M.

Warshay

, “

Effects of boundary layer buildup in shock tubes upon chemical rate measurements

,”

Report No. NASA-TN-4795

(

NASA

,

Cleveland, OH

,

1968

).

Google Scholar

31.

C. J.

Doolan

and

P. A.

Jacobs

, “

Modeling mass entrainment in a quasi-one-dimensional shock tube code

,”

AIAA J.

34

,

1291

–

1293

(

1996

).

https://doi.org/10.2514/3.13227

Google Scholar

Crossref

32.

M.

Satchell

,

A.

Glenn

,

P.

Collen

,

R.

Penty-Geraets

,

M.

McGilvray

, and

L.

di Mare

, “

Analytical method of evaluating nonuniformities in shock tube flows: Application

,”

AIAA J.

60

,

669

–

676

(

2022

).

https://doi.org/10.2514/1.J060991

Google Scholar

Crossref

33.

M.

Satchell

,

L.

di Mare

, and

M.

McGilvray

, “

Flow nonuniformities behind accelerating and decelerating shock waves in shock tubes

,”

AIAA J.

60

,

1537

–

1548

(

2022

).

https://doi.org/10.2514/1.J060375

Google Scholar

Crossref

34.

P. A.

Gnoffo

,

R. N.

Gupta

, and

J. L.

Shinn

, “

Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium

,”

Report No. NASA-TP-2867

(

NASA

,

1989

).

Google Scholar

35.

S.

Brody

,

K. S.

Lau

,

J.

Clarke

,

M.

McGilvary

, and

L. D.

Mare

, “

Numerical simulation of transpiration cooling on stagnation line in thermochemical non-equilibrium

,” AIAA Paper No. 2024-0648,

2024

.

Google Scholar

36.

M.

Morduchow

and

P. A.

Libby

, “

On a complete solution of the one-dimensional flow equations of a viscous, heat-conducting, compressible gas

,”

J. Aeronaut. Sci.

16

,

674

–

684

(

1949

).

https://doi.org/10.2514/8.11882

Google Scholar

Crossref

37.

S.

Gordon

and

B. J.

McBride

, “

Computer program for calculation of complex chemical equilibrium

,”

Report No. NASA-RP-1311

,

1994

.

Google Scholar

38.

R. J.

Gollan

, “

The computational modelling of high-temperature gas effects with application to hypersonic flows

,” Ph.D. thesis (

University of Queensland

,

2008

).

Google Scholar

39.

P. V.

Marrone

and

C. E.

Treanor

, “

Chemical relaxation with preferential dissociation from excited vibrational levels

,”

Phys. Fluids

6

,

1215

–

1221

(

1963

).

https://doi.org/10.1063/1.1706888

Google Scholar

Crossref

40.

R. N.

Gupta

,

J. M.

Yos

,

R. A.

Thompson

, and

K.-P.

Lee

, “

A review of reaction rates and thermodynamic and transport properties for an 11 species air model for chemical and thermal nonequilibrium calculations to 30 000 K

,”

Report No. NASA-RP-1232

(

NASA

,

1990

).

Google Scholar

41.

D.

Giordano

,

V.

Bellucci

,

G.

Colonna

,

M.

Capitelli

,

I.

Armenise

, and

C.

Bruno

, “

Vibrationally relaxing flow of N₂ past an infinite cylinder

,”

J. Thermophys. Heat Transfer

11

,

27

–

35

(

1997

).

https://doi.org/10.2514/3.854

Google Scholar

Crossref

42.

V.

Blackman

, “

Vibrational relaxation in oxygen and nitrogen

,”

J. Fluid Mech.

1

,

61

–

85

(

1956

).

https://doi.org/10.1017/S0022112056000056

Google Scholar

Crossref

43.

R. C.

Millikan

and

D. R.

White

, “

Systematics of vibrational relaxation

,”

J. Chem. Phys.

39

,

3209

–

3213

(

1963

).

https://doi.org/10.1063/1.1734182

Google Scholar

Crossref

44.

P. C.

De Boer

, “

Curvature of shock fronts in shock tubes

,”

Phys. Fluids

6

,

962

–

971

(

1963

).

https://doi.org/10.1063/1.1706852

Google Scholar

Crossref

45.

S. P.

Sharma

and

W.

Gillespie

, “

Nonequilibrium and equilibrium shock front radiation measurements

,” AIAA Paper No. 90-0139,

1990

.

https://doi.org/10.2514/6.1990-139

Google Scholar

46.

J. G.

Kim

and

I. D.

Boyd

, “

Modeling of strong nonequilibrium in nitrogen shock waves

,” AIAA Paper No. 2013-3150,

2013

.

https://doi.org/10.2514/6.2013-3150

Google Scholar

47.

C.

Park

, “

Review of chemical-kinetic problems of future NASA missions, I: Earth entries

,”

J. Thermophys. Heat Transfer

7

,

385

–

398

(

1993

).

https://doi.org/10.2514/3.431

Google Scholar

Crossref

48.

C.

Park

,

J. T.

Howe

,

R. L.

Jaffe

, and

G. V.

Candler

, “

Review of chemical-kinetic problems of future NASA missions, II: Mars entries

,”

J. Thermophys. Heat Transfer

8

,

9

–

23

(

1994

).

https://doi.org/10.2514/3.496

Google Scholar

Crossref

49.

B. A.

Cruden

and

A. M.

Brandis

, “

Measurement of radiative non-equilibrium for air shocks between 7–9 km/s

,” AIAA Paper No. 2017-4535,

2017

.

https://doi.org/10.2514/6.2017-4535

Google Scholar

50.

C. O.

Johnston

and

A. M.

Brandis

, “

Modeling of nonequilibrium CO fourth-positive and CN violet emission in CO₂-N₂ gases

,”

J. Quant. Spectrosc. Radiat. Transfer

149

,

303

–

317

(

2014

).

https://doi.org/10.1016/j.jqsrt.2014.08.025

Google Scholar

Crossref

51.

A. M.

Brandis

and

B. A.

Cruden

, “

Titan atmospheric entry radiative heating

,” AIAA Paper No. 2017-4534,

2017

.

Google Scholar

52.

A. B.

Glenn

,

P. L.

Collen

, and

M.

McGilvray

, “

Radiation measurements of shockwaves in synthetic air and pure nitrogen

,” AIAA Paper No. 2024-0226,

2024

.

https://doi.org/10.2514/6.2024-0226

Google Scholar

53.

E. E.

Whiting

,

P.

Chul

,

Y.

Liu

,

O.

Arnold

, and

A.

Paterson

, “

NEQAIR96, nonequilibrium and equilibrium radiative transport and spectra program: User's manual

,”

Report No. NASA-RP-1389

(

NASA

,

1996

).

Google Scholar

54.

M. J.

Wright

,

T.

White

, and

N.

Mangini

, “

Data parallel line relaxation (DPLR) code user manual Acadia-Version 4.01.1

,”

Report No. NASA-TM-2009-215388

(

NASA

,

2009

).

Google Scholar

55.

P. L.

Collen

,

L.

Di Mare

,

M.

McGilvray

, and

M.

Satchell

, “

Analysis of shock deceleration effects on radiation experiments in the NASA electric arc shock tube

,” AIAA Paper No. 2022-0267,

2022

.

https://doi.org/10.2514/6.2022-0267

Google Scholar

56.

C.-Y.

Wen

and

H. G.

Hornung

, “

Non-equilibrium dissociating flow over spheres

,”

J. Fluid Mech.

299

,

389

–

405

(

1995

).

https://doi.org/10.1017/S0022112095003545

Google Scholar

Crossref

Published open access through an agreement withJISC Collections

Quasi-one-dimensional non-equilibrium method for shock tube and stagnation line flows

I. INTRODUCTION

II. NUMERICAL FORMULATION

A. Shock tube flow

B. Discretization and solution method

III. VALIDATION

A. Shock structure

B. Thermochemistry

1. Chemistry

2. Translational–vibrational relaxation

C. Non-reacting shock tube test

D. Reacting shock tube case

IV. IMPACT OF SHOCK STRUCTURE AND BOUNDARY LAYER EFFECTS

A. 6.2 km/s shock through 133.3 Pa nitrogen

B. 6.2 km/s shock through 13.3 Pa nitrogen

V. COMPARISON TO EXPERIMENTS

A. Radiance model

B. NASA-EAST experiment for titan entry

C. T6 low earth orbit return experiments

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

NOMENCLATURE

NOMENCLATURE

APPENDIX A: DISCRETIZED SYSTEM OF EQUATIONS

APPENDIX B: PROBLEMS DESCRIBED BY DIFFERING OUTFLOW MODELS

APPENDIX C: MIRELS RADIAL VELOCITY PROFILE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Quasi-one-dimensional non-equilibrium method for shock tube and stagnation line flows

I. INTRODUCTION

II. NUMERICAL FORMULATION

A. Shock tube flow

B. Discretization and solution method

III. VALIDATION

A. Shock structure

B. Thermochemistry

1. Chemistry

2. Translational–vibrational relaxation

C. Non-reacting shock tube test

D. Reacting shock tube case

IV. IMPACT OF SHOCK STRUCTURE AND BOUNDARY LAYER EFFECTS

A. 6.2 km/s shock through 133.3 Pa nitrogen

B. 6.2 km/s shock through 13.3 Pa nitrogen

V. COMPARISON TO EXPERIMENTS

A. Radiance model

B. NASA-EAST experiment for titan entry

C. T6 low earth orbit return experiments

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

NOMENCLATURE

NOMENCLATURE

APPENDIX A: DISCRETIZED SYSTEM OF EQUATIONS

APPENDIX B: PROBLEMS DESCRIBED BY DIFFERING OUTFLOW MODELS

APPENDIX C: MIRELS RADIAL VELOCITY PROFILE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only