1. Introduction
Hypersonic vehicles travel at speeds that are sufficient to generate a strong bow shock. At such high flight speeds, the bow shock raises the temperature of the gas in the shock layer to a level sufficient to cause dissociation and ionisation of molecules and atoms in the gas. In these conditions, the macroscopic properties of the flow depart significantly from those predicted by ideal gas theory. The shock layer is said to be in a state of non-equilibrium when the thermochemical behaviour of the gas evolves on a time scale comparable to the convective time scales of the flow.
Uncertainty in gas non-equilibrium behaviour contributes significantly to overall uncertainty in vehicle surface heat flux, making the design of mass-efficient vehicles a challenge (Palmer, Prabhu & Cruden Reference Palmer, Prabhu and Cruden2014; West et al. Reference West, Brune, Hosder and Johnston2016, Reference West, Johnston and Hosder2017; Carter & Boyd Reference Carter and Boyd2024; Rataczak, Boyd & McMahon Reference Rataczak, Boyd and McMahon2024). Heat transfer via electromagnetic radiation is of particular concern due to its strong dependence on the thermochemical state of the gas. Multi-temperature models are required to accurately describe the distribution of energy within the vibrational, electronic, rotational and translational modes available within the gas, the most widely used being that of Park (Reference Park1989b ). The multi-temperature description of the gas is often combined with a chemistry model which considers that species densities evolve according to a set of temperature-dependent reaction rates, often taking the Arrhenius form:
\begin{equation} k = A T^{n}{\textrm e}^{-E/RT}, \end{equation}
where
$k$
is the reaction rate (most often expressed in units of
$\textrm{cm}^{3}\,mol^{-1}\,s^{-1}$
),
$A$
is a pre-exponential coefficient,
$n$
is a non-dimensional power,
$E$
is the activation energy for the reaction,
$R$
is the universal gas constant and
$T$
is a temperature. In practice, the parameters
$A$
and
$n$
are obtained by tuning against experimental results. Shock tubes are one important source of experimental measurements for this purpose.
The aim of a shock tube experiment is to generate a flow with the same shock layer that would be encountered by a flight vehicle. To do this, it is required that the shock speed, gas composition and gas density match those of the flight case. Shock tubes feature a high-pressure driver region, a low-pressure driven region, a diaphragm separating the two and a station at the downstream end of the facility where the state of the shock layer is observed. The wave processes that occur during a shock tube experiment are represented in figure 1.
The experiment begins with the rupture of the primary diaphragm, allowing the high-pressure driver gas to unsteadily expand into the driven tube. This drives the test gas forward at a speed sufficient to generate a strong normal shock. The trajectory of the shock is monitored during the experiment to obtain information about the actual experimental conditions. This can be achieved by a number of means: as an example, the driven section can be instrumented with many time-of-arrival sensors, typically diaphragm-type piezoelectric pressure transducers, and the speed of the shock can be determined by finite difference. The length of the test slug increases along the length of the tube until it reaches a maximum value governed by the growth of the boundary layer. Upon reaching the observation window, measurements of the radiating shock layer are made by means of optical emission spectroscopy or a similar technique. Analysis of these measurements with a suitable radiation model allows species densities and temperatures to be estimated as a function of distance behind the shock. If constant post-shock properties are assumed, the post-shock distance is linearly related to the particle time of flight; in this way reaction rate constants may be derived.
Figure 1. Position–time diagram of wave processes in a generic shock tube. The ideal trajectory represents the assumption of a constant shock speed equal to the final observed shock speed, whereas the real trajectory represents the true shock profile in the presence of attenuation. These result in different particle times of flight at the instant of observation, potentially leading to the inference of different apparent reaction rate constants. Here s-x and us-x represent the steady and unsteady expansion wave processes that occur at the diaphragm station.
Real shock tubes have been shown to exhibit a number of non-ideal behaviours that complicate the interpretation of experimental data. If these effects are not considered, the reaction rates inferred from experiment will be incorrect. The primary non-ideal effect is growth of a boundary layer in the test slug and driver gas which serves to attenuate the shock and decrease the post-shock velocity by entraining mass from the core flow (Mirels Reference Mirels1963). This effect is most prominent at lower pressures in tubes with smaller cross-sections (Clarke et al. Reference Clarke, Di Mare and McGilvray2023b ). Additionally, the behaviour of the diaphragm and driver gas can introduce complex compression/rarefaction wave patterns that travel downstream and interact with the test slug before it reaches the observation window (White Reference White1958; Bowman Reference Bowman1966; Morgan & Gildfind Reference Morgan and Gildfind2013). The nature of these wave patterns is dependent on many factors such as the driver technology employed, the diaphragm thickness, material and score pattern. Upon reaching the shock front, the compression/rarefaction wave will modulate the shock speed. These phenomena are intrinsically bound to the processes taking place at the start of the experiment, such as details of the diaphragm opening process. Such details do not take place in identical fashion in experiments with identical nominal conditions and cause shock tube experiments to be intrinsically not reproducible.
The consequence of these effects is that real shock tube experiments for the purpose of determining reaction rate constants:
-
(i) feature non-uniform shock speeds and
-
(ii) are difficult to repeat.
The presence of a non-uniform shock speed profile or trajectory gives rise to a complex post-shock velocity profile and thus particle time of flight. Consider figure 2, which shows the time of flight as a function of post-shock distance for a 10 km s−1 shock in 13.3 Pa of synthetic air with and without non-ideal effects. Two cases of linear shock acceleration and deceleration of 500
$\textrm{m s}^{-1}\,m^{-1}$
, but still arriving at the same final speed of 10
$\textrm{km s}^{-1}$
, are presented. The particle time of flight is 45 % lower and 120 % higher than the constant value by the end of the test slug for the decelerating and accelerating cases, respectively. If reaction rates were optimised against real shock tube data (assuming the shock speed is constant throughout the trajectory) using an inviscid one-dimensional (1-D) model for this case, a particle at 30 mm post-shock would erroneously be believed to have been in flight for 42 % less time than in reality, leading to a deduced reaction rate that is 42 % faster than required. In a real shock tube experiment, the particle time of flight is more complex than shown in figure 2 as the shock trajectory is often not well described by a linear trend. Therefore, the use of numerical models featuring the true shock trajectory is required to account for the complex hydrodynamic processes taking place and thus allow accurate determination of thermochemical reaction rates on a test-by-test basis.
Figure 2. Time of flight of particles post-shock for a 10 km s−1 shock in 13.3 Pa of synthetic air. The shock tube diameter is 100 mm in this case. (a) Time of flight. (b) Difference. Reproduced from Clarke et al. (Reference Clarke, Di Mare and McGilvray2023b ).
Several authors have attempted to isolate the effects of non-uniform shock speed on test slug properties. Holbeche (Reference Holbeche1964) observed that temperatures behind an attenuating shock exceeded equilibrium predictions by up to 200 K and derived an expression to predict these discrepancies. Light (Reference Light1973) investigated trends in post-shock electron number density and attributed differences between non-equilibrium theory and experiment to the presence of expansion waves in the test gas. The expansion waves were believed to originate from sharp pressure decrease caused by the electric arc driver technology. Similar behaviour was observed by Cruden (Reference Cruden2012), where electron number densities were found to be significantly in excess of those expected at chemical equilibrium. Brandis et al. (Reference Brandis, Cruden, Prabhu, Bose, McGilvray, Morgan and Morgan2010) performed a cross-facility comparison and proposed that the excess electron number density and radiation observed in the NASA Electric Arc Shock Tube (EAST) facility was due to strong shock deceleration. This was shown to be the case by Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022, Reference Collen, Di Mare, McGilvray and Satchell2023) through the use of a numerical tool developed by Satchell et al. (Reference Satchell, Glenn, Collen, Penty-Geraets, McGilvray and Di Mare2022b ) called LAgrangian Shock Tube Analysis (LASTA 1.0). Satchell’s method used equilibrium thermochemistry to model the test slug based on the experimental shock trajectory alone. This approach assumed that entropy and enthalpy were constant in the Lagrangian reference frame and that the pressure and enthalpy were algebraically related, or that sound wave relaxation within the slug is instantaneous. Additionally, a restriction on the pressure behind the shock wave was imposed such that it could only increase monotonically with the intention of preventing non-physical pressure profiles. This is referred to in the rest of this paper as the pressure-limiter effect. The tool relied on an assumption of thermochemical equilibrium and did not enforce conservation of momentum. These limitations precluded its application to shock tube cases where waves from the driver are significant or where the test slug was predominantly in a non-equilibrium state, i.e. the cases of primary interest for shock tube simulation.
We aim to demonstrate with this work a novel, computationally cheap numerical method that decouples non-ideal fluid mechanics in shock tubes from the object of study, the non-equilibrium thermochemistry. This approach takes account of both forward- and backward-running waves by constraining the model with two boundary conditions derived from experiment: the measured shock speed history and a wall static pressure history at the observation window. The second boundary condition is imposed at the end of the domain and its effect on the test slug is included by reverse time integration. The system of governing equations is cast, discretised and solved in characteristic form. This work offers a first platform for the optimisation of thermochemical reaction rates against shock tube measurements with non-ideal wave processes removed.
The paper begins with a detailed description of model domain, governing equations and method of solution in § 2. The results are then presented including a validation of the basic fluid model in ideal gas (§ 3), the two-temperature model in synthetic air and finally application to a variety of shock tube test cases in air and Titan-entry-relevant mixtures where non-equilibrium and non-uniform shock speed are important considerations (§ 4). A discussion of the key findings is provided in § 5, with attention paid to the sensitivity of the final test slug to the shock trajectory, the chosen reaction rate constants and the combined effect on the reaction rates that may ultimately be deduced. Finally, the major findings of the work are summarised in § 6 and recommendations for future work are made.
2. Methodology
The model idealises the test slug as a 1-D, inviscid core surrounded by a growing boundary layer, as shown in figure 3. The model is conceived in the shock reference frame with
$x$
representing the distance behind the shock,
$x \in [0, \ell ]$
, where
$0$
is the shock front and
$\ell$
is the contact surface. The shock travels left to right at speed
$-U_s$
in the laboratory frame. The test gas travels at speed
$u$
in the laboratory frame and
$u + U_s = u'$
in the shock frame. The 1-D Euler equations are solved at each time point
$t$
along the supplied experimental shock trajectory with a source term to represent mass removal by the boundary layer. The domain size or slug length,
$\ell$
, is allowed to vary along the length of the shock trajectory based on the mass balance between the shock addition and boundary-layer subtraction. The size of the test slug and hence the test time is determined by the location where the mass flow consumed by the boundary layer matches the mass flow processed by the shock. Two boundary conditions are imposed on the model:
-
(i) The primitive variables (
$\rho$
,
$u$
,
$p$
) immediately behind the shock at each time point are set by the shock trajectory. -
(ii) The pressure profile (in space) at the final time point is set from static pressure measurements, e.g. obtained from a transducer.
Figure 3. Diagram of the model domain, shown in the shock-attached reference frame. Here,
$S$
is the core flow diameter,
$\ell$
is the test slug length,
$u$
is the velocity,
$d$
is the tube diameter,
$U_s$
is the shock speed and the subscripts
$s$
and
$w$
refer to the free stream and the wall, respectively.
2.1. Governing equations for the two-temperature plasma
The method solves the 1-D species, momentum, vibrational–electronic energy and total energy conservation equations in characteristic form. The Park two-temperature model (Park Reference Park1989a
) was implemented, although this is not a requirement as any thermochemical model could be utilised with the same approach. This model assumes that the rotational and translational temperatures of heavy particles are in equilibrium at a temperature
$T$
, while the vibrational temperature of molecular species, the electron translational temperature and the electronic excitation of atoms and molecules are in equilibrium at temperature
$T_V$
. The three-dimensional conservation equations are published by Gnoffo, Gupta & Shinn (Reference Gnoffo, Gupta and Shinn1989). For
$\mathcal{N}$
species, the 1-D system is
\begin{equation} \frac {\partial }{\partial t'} \underbrace {\begin{bmatrix} \rho _s\\ \vdots \\ \rho _{\mathcal{N}} \\ \rho u\\ \rho e_V \\ \rho E \end{bmatrix}}_{1} = -\frac {\partial }{\partial x} \underbrace {\begin{bmatrix} \rho _s u\\ \vdots \\ \rho _{\mathcal{N}} u \\ \rho uu + p\\ \rho e_V u + p_e\\ \rho H u \end{bmatrix}}_{2} - \underbrace {\frac {2v}{r} \begin{bmatrix} \rho _s/\rho \\ \vdots \\ \rho _{\mathcal{N}}/\rho \\ u\\ e_V \\ H \\ \end{bmatrix}}_{3} + \underbrace {\begin{bmatrix} \dot {f}_s\\ \vdots \\ \dot {f}_{\mathcal{N}}\\ 0\\ \dot {f}_{V} \\ 0 \\ \end{bmatrix}}_{4}. \end{equation}
In (2.1) the terms labelled 1 are the conserved variables,
$\rho _s$
denotes the density of species
$s$
with the bulk density given by
\begin{equation} \rho = \sum _{s=1}^{\mathcal{N}} \rho _s ,\end{equation}
$u$
is the velocity,
$e_V$
is the mixture vibrational–electronic energy per unit mass and
$E$
is the mixture total energy per unit mass. The terms labelled 2 comprise the flux vector where
$H$
is the mixture total enthalpy per unit mass and
$p$
is the pressure. The pressure obeys Boyle’s law:
\begin{equation} p = \sum _{s=1}^{\mathcal{N}} p_s = R \left ( \sum _{s=2}^{\mathcal{N}} \frac {\rho _{s}T}{M_s} + \frac {\rho _e T_V }{M_e}\right ). \end{equation}
The source terms labelled 3 are included to model the removal of mass from the core flow by the boundary layer. Here,
$v$
is the radial component of velocity and
$r$
is the radius of the tube. The boundary-layer model is presented in greater detail in § 2.2. The source terms labelled 4 account for the effect of thermochemistry. The terms
$\dot {f}_s, \dot {f}_n$
account for the creation or destruction of a particular species due to chemical reactions. The term
$\dot {f}_V$
is comprised of several contributions:
\begin{equation} \begin{aligned} \dot {f}_V &= \underbrace {\sum _{s=mol.} \rho _s \frac {(e_{v,s}^* - e_{v,s})}{\lt \tau _s \gt }}_{5} + \underbrace {2 \rho _e \frac {3}{2} \bar {R} (T - T_V) \sum _{s=2}^{\mathcal{N}} \frac {\upsilon _{es}}{M_s}}_{6} + \underbrace {\sum _{s=mol.} \dot {f}_s \hat {D}_s}_{7} -\underbrace {\sum _{s=mol.} \dot {n}_{e,s} \hat {I}_s}_{8}. \end{aligned} \end{equation}
Term 5 represents the energy exchange or relaxation between trans-rotational and vibro-electronic energy modes, where
$\tau _s$
is the relaxation time. Time
$\tau _s$
is obtained from correlations in the literature and is dependent on the gas mixture. Full details of the relaxation time constant calculation are provided in Appendix C. Parameter
$e_{v,s}$
is calculated using the method in Gnoffo et al. (Reference Gnoffo, Gupta and Shinn1989) which is also detailed in Appendix D. Term 6 represents the energy exchange between modes due to elastic collisions between electrons and heavy particles. Term 7 is the energy lost or gained due to dissociation reactions and, similarly, term 8 is the energy lost or gained due to electron impact ionisation.
It is desirable to manipulate the governing equations such that they are expressed in terms of a set of variables of interest or so-called primitive variables. In order to cast the system with
$T_V$
as a primitive variable, we must express all
$e_V$
terms as functions of
$T_V$
and other primitive variables only. Neglecting chemical source terms, the vibrational–electronic energy conservation equation is
\begin{equation} \frac {\partial (\rho e_V)}{\partial t} + \frac {\partial (\rho e_V u)}{\partial x} + p_e \frac {\partial u}{\partial x} = 0 .\end{equation}
In a problem with many species
$s$
,
$e_V$
may be evaluated as
\begin{align} e_V &= \sum _{s=1}^{\mathcal{N}} \frac {\rho _s}{\rho } e_{V,s}, \end{align}
\begin{align} \rho e_V &= \sum _{s=1}^{\mathcal{N}} \rho _s e_{V,s}. \end{align}
If this definition is substituted into (2.5):
\begin{eqnarray} \begin{aligned} & \frac {\partial }{\partial t} \left (\sum _{s=1}^{\mathcal{N}} \rho _s e_{V,s} \right ) + u \frac {\partial }{\partial x} \left (\sum _{s=1}^{\mathcal{N}} \rho _s e_{V,s} \right ) + \frac {\partial u }{\partial x} \sum _{s=1}^{\mathcal{N}} \rho _s e_{V,s} + p_e \frac {\partial u}{\partial x} = 0 ,\\ & \sum _{s=1}^{\mathcal{N}} \frac {\partial \rho _s}{\partial t} e_{V,s} + \sum _{s=1}^{\mathcal{N}} \rho _s C_{v,V}^{s} \frac {\partial T_V}{\partial t} + u \left ( \sum _{s=1}^{\mathcal{N}} \frac {\partial \rho _s}{\partial x} e_{V,s} + \sum _{s=1}^{\mathcal{N}} \rho _s C_{v,V}^{s} \frac {\partial T_V}{\partial x} \right ) \\ &\quad +\, \frac {\partial u }{\partial x} \sum _{s=1}^{\mathcal{N}} \rho _s e_{V,s} + p_e \frac {\partial u}{\partial x} = 0. \end{aligned} \end{eqnarray}
Subtracting the product of the continuity equation with
$e_{V,s}$
:
\begin{align} \sum _{s=1}^{\mathcal{N}} \rho _s C_{v,V}^{s} \frac {\partial T_V}{\partial t} + u \sum _{s=1}^{\mathcal{N}} \rho _s C_{v,V}^{s} \frac {\partial T_V}{\partial x} + p_e \frac {\partial u}{\partial x} = 0 , \end{align}
\begin{align} \rho C_{v,V} \frac {\partial T_V}{\partial t} + u \rho C_{v,V} \frac {\partial T_V}{\partial x} + p_e \frac {\partial u}{\partial x} = 0 , \end{align}
\begin{align} \frac {\partial T_V}{\partial t} + u \frac {\partial T_V}{\partial x} + \frac {p_e}{\rho C_{v,V}} \frac {\partial u}{\partial x} = 0. \end{align}
The final system in primitive form is then
\begin{align} \frac {\partial }{\partial t'} \begin{bmatrix} \rho _s\\ \vdots \\ \rho _{\mathcal{N}} \\ u\\ T_V \\ p \end{bmatrix} = -\begin{bmatrix} u & \ldots & 0 & \rho _1 & 0 & 0 \\ \vdots & \ddots &\vdots & \vdots & \vdots &\vdots \\ 0 & \ldots &u & \rho _{\mathcal{N}} & 0 & 0\\ 0 &\ldots &0 & u & 0 & 1/\rho \\ 0 &\ldots &0 & \frac {p_e}{\rho C_{v,V}} & u & 0\\ 0 &\ldots &0 & \rho a^{2} & 0 & u \end{bmatrix} \frac {\partial }{\partial x} \begin{bmatrix} \rho _s\\ \vdots \\ \rho _{\mathcal{N}} \\ u\\ T_V \\ p \end{bmatrix} - \rho \frac {2v}{r} \begin{bmatrix} \rho _s/\rho \\ \vdots \\ \rho _{\mathcal{N}}/\rho \\ 0\\ 0 \\ a^2 \\ \end{bmatrix} + \begin{bmatrix} \dot {f}^{\prime}_s\\ \vdots \\ \dot {f}_{\mathcal{N}}^{\prime} \\ 0\\ \dot {f}^{\prime}_{V} \\ 0 \\ \end{bmatrix}.\\[-20pt]\nonumber \end{align}
Here, the electron pressure
$p_e$
appears which may be readily obtained from the vibrational–electronic temperature:
\begin{equation} p_e = \frac {T_V R \rho _e}{M_e}, \end{equation}
where
$M_e$
is the electron molar mass and
$\rho _e$
is the electron density which is itself a primitive variable. The thermochemical source terms are marked with a prime to illustrate that a transformation between conservation and primitive form has occurred. Finally, the system can be expressed in a characteristic basis by solving the eigenvalue problem for the homogeneous system including only the left-hand side and the first term on the right-hand side in (2.12):
\begin{align} \frac {\partial }{\partial t^{\prime}} \begin{bmatrix} w_1\\ \vdots \\ w_{\mathcal{N}} \\ w_{\mathcal{N}+1}\\ w_{\mathcal{N}+2} \\ w_{\mathcal{N}+3} \end{bmatrix} = -\begin{bmatrix} u & \ldots & 0 & 0 & 0 & 0 \\ \vdots & \ddots &\vdots & \vdots & \vdots &\vdots \\ 0 & \ldots &u &0 & 0 & 0\\ 0 &\ldots &0 & u + a& 0 & 0\\ 0 &\ldots &0 & 0 & u & 0\\ 0 &\ldots &0 & 0 & 0 & u - a \end{bmatrix} \frac {\partial }{\partial x} \begin{bmatrix} w_1\\ \vdots \\ w_{\mathcal{N}} \\ w_{\mathcal{N}+1}\\ w_{\mathcal{N}+2} \\ w_{\mathcal{N}+3} \end{bmatrix} - a \frac {2v}{r} \begin{bmatrix} 0\\ \vdots \\ 0\\ 1\\ 0 \\ -1\\ \end{bmatrix} + \begin{bmatrix} \dot {f}^{\prime \prime}_s\\ \vdots \\ \dot {f}^{\prime \prime}_{\mathcal{N}}\\ 0\\ \dot {f}^{\prime \prime}_{V} \\ 0 \\ \end{bmatrix}.\\[-20pt]\nonumber \end{align}
Here,
$w_{s}$
are the characteristic variables. The characteristic variables provide a clear representation of the underlying fluid processes influencing the test slug properties, namely wave processes travelling back and forth through the gas. The speed of sound,
$a$
, is evaluated using the frozen expression from Gnoffo et al. (Reference Gnoffo, Gupta and Shinn1989):
\begin{align} a &= (1 + \beta ) \frac {p}{\rho }, \end{align}
\begin{align} \beta &= \frac {\partial p}{\partial \rho E} = \frac {R}{\rho C_{V,tr}} \sum _{s=2}^{\mathcal{N}} \frac {\rho _s}{M_s}, \\[9pt] \nonumber \end{align}
where
$C_{V,tr}$
is the specific heat at constant volume of the mixture for translational–rotational energy in
$\textrm{J kg}^{-1}\,K^{-1}$
. Details of the solution of the eigenvalue problem are included in Appendix A.
2.2. Boundary layer
The laminar boundary layer is represented through its effects on the core flow. This is included by assuming the radial velocity
$v$
behaves according to a local similarity approximation proposed by Mirels (Reference Mirels1963). Using this approximation, the radial velocity at the edge of the boundary layer can be evaluated as
\begin{equation} v = \frac {\rho _{e,0} u_{e,0}}{4 \rho \sqrt {\ell \ell _{m}}}, \end{equation}
where
$\ell$
is the distance from the shock,
$\rho _{e,0}$
and
$u_{e,0}$
are the immediate post-shock density and velocity in the shock reference frame,
$\rho$
is the density in the test slug evaluated at
$\ell$
and
$\ell _{m}$
is the maximum slug length. The maximum slug length is given by Mirels as
\begin{equation} \ell _{m} = \frac {d^{2}}{16 \beta ^{2}} \left ( \frac {\rho _{e,0}}{\rho _{w,0}} \right )^{2} \frac {u_{e,0}}{U_s - u_{e,0}} \frac {u_{e,0}}{\upsilon _{w,0}},\end{equation}
where
$d$
is the tube diameter,
$\beta$
is a calculated shape factor,
$\rho _{w,0}$
and
$\upsilon _{w,0}$
are the density and kinematic viscosity evaluated at the tube wall temperature and
$U_s$
is the shock speed as before. Mirels’s analysis provides estimates of the shape factor
$\beta$
behind a shock wave:
\begin{equation} \beta = C_{M}^{\alpha } 1.59 \left ( 1 + \frac {1.796 + 0.802 W}{WZ - 1} \right ), \end{equation}
where
\begin{equation} C_M = \frac {\rho _{sh} \mu _{sh}}{\rho _{w,sh} \mu _{w}}. \end{equation}
Here,
$W = U_s/ u_{e,0}$
and
$Z = {(\gamma + 1)}/{(\gamma - 1)}$
. Term
$C_M$
corrects for boundary-layer shape to account for changes in the free-stream properties. Viscosities,
$\mu$
, are determined using Sutherland’s law. These equations are evaluated to obtain the radial velocity profile along the test slug which can then be used to evaluate the source term. It should be noted that Mirels’s analysis was developed for a steady, constant shock speed case, whereas here it is applied for a non-uniform, unsteady shock speed, the physical implication being that it is assumed that the entire boundary-layer shape is instantly resized for changing shock speed. Previous work has shown that this assumption is sufficient for many typical shock tube cases (Satchell et al. Reference Satchell, Glenn, Collen, Penty-Geraets, McGilvray and Di Mare2022b
).
2.3. Characteristics of the two-temperature non-equilibrium system
The characteristics
$w^{-}$
and
$w^{+}$
represent sound waves travelling through the test slug,
$w^{-}$
originating from the shock front and
$w^{+}$
originating from the driver gas. Figure 4(
a) illustrates the directions of propagation for the
$w^{-}$
and
$w^{+}$
waves in a shock tube experiment where time progresses from top to bottom. Accurate numerical modelling of the effect of the
$w^{+}$
characteristic on the test slug properties with forward time integration would require a pressure boundary condition at the contact discontinuity. Unfortunately, experimental information about the position and pressure of the contact surface is difficult to obtain non-intrusively and is typically unavailable; thus, the effect of
$w^{+}$
on the test slug cannot be calculated via forward time integration. However,
$w^{+}$
waves will ultimately propagate through the test slug to the shock front, so
$w^{+}$
is known indirectly through its effect on the future trajectory. By integrating
$w^{+}$
backwards in time, the approach taken here and shown in figure 4(
b), the test gas is within the domain of influence of the shock with respect to
$w^{+}$
and the full test slug history may be calculated using the shock trajectory and final pressure profile as boundary conditions.
Figure 4. Time integration for the
$w^{+}$
characteristic. (a) Forward integration. (b) Backward integration.
Figure 5. Diagrams of the wave processes modelled explicitly in the method of (a) Satchell et al. (Reference Satchell, Glenn, Collen, Penty-Geraets, McGilvray and Di Mare2022b
) and (b) the present work. The diagrams are representative of the processes occurring in a real experiment, but are not drawn to scale. The present work offers improvement in the region where
$u+a$
waves may propagate.
The
$w^{+}$
characteristic waves travel at a finite speed and thus their inclusion via the pressure boundary condition will influence only a subset of the solution domain. This is illustrated by figure 5, an x–t diagram of a shock tube experiment in the shock reference frame. Figure 5(
a) shows the wave processes modelled by Satchell’s method and figure 5(
b) shows the additional wave processes modelled in the present work as a result of the new boundary condition. It can be seen that the influence of the
$u + a$
waves on the slug properties becomes more significant towards the rear of the slug as more waves have had time to pass through this region. Necessarily, the influence of
$u - a$
waves, although still significant, becomes less important towards the rear of the slug. Satchell’s method assumes that the relative flow within the test slug is low subsonic and that the time it takes for sound waves to propagate along the slug length is negligible compared with the duration of the trajectory, or more precisely that within a typical acoustic propagation time variations in shock speed are negligible. This assumption is addressed by the proposed method and improvement is expected in the region where the
$u + a$
waves may propagate.
It should be noted that, although pressure measurements are used as a boundary condition here, it is possible in principle to use other experimental data, such as a centreline temperature measurement. It is also possible to move the boundary condition upstream or downstream, assuming experimental data are available, to the location of greatest interest.
2.4. Numerical scheme
Because the system is solved in Lagrangian form, a frame transformation between the stationary and inertial frames is required. Define the following:
\begin{align} {\textrm d}x &= {\textrm d}y-U_s {\textrm d}t, \end{align}
\begin{align} {\textrm d}\phi &= \dfrac {\partial \phi }{\partial x}{\textrm d}x+ \dfrac {\partial \phi }{\partial t} {\textrm d}t, \end{align}
\begin{align} {\textrm d}\phi &= \dfrac {\partial \phi }{\partial x}\left ( {\textrm d}y-U_s {\textrm d}t\right )+ \dfrac {\partial \phi }{\partial t} {\textrm d}t, \end{align}
\begin{align} {\textrm d}\phi &= \dfrac {\partial \phi }{\partial x} {\textrm d}y + \left ( \dfrac {\partial \phi }{\partial t} -U_s \dfrac {\partial \phi }{\partial x} \right ){\textrm d}t . \end{align}
Time derivative measured by an observer bound to the frame of the shock:
\begin{align} \dfrac {\partial \phi }{\partial t^{\prime}} &= \dfrac {\partial \phi }{\partial t} -U_s \dfrac {\partial \phi }{\partial x}, \end{align}
\begin{align} \dfrac {D}{Dt} &= \dfrac {\partial }{\partial t}+ \left ( u^{\prime}+U_s\right )\dfrac {\partial }{\partial x}= \dfrac {\partial }{\partial t}+ u\dfrac {\partial }{\partial x}. \end{align}
The system of governing equations is solved using a finite-difference method. The forward characteristics are evaluated at the
$i$
th spatial point and
$k$
th time point as
\begin{align} \delta ^t w_{i,k}^{(1 \rightarrow \mathcal{N}, \mathcal{N}+2)}+ \dfrac {\Delta t}{\Delta x} u^{\prime}_{i,k} \delta ^+ w_{i,k}^{(1 \rightarrow \mathcal{N}, \mathcal{N}+2)} + \dot {f}^{\prime \prime} & =0, \end{align}
\begin{align} \delta ^t w_{i,k}^{(\mathcal{N}+1)}+ \dfrac {\Delta t}{\Delta x}\left ( u^{\prime}_{i,k}+a_{i,k} \right )\delta ^+ w_{i,k}^{(\mathcal{N}+1)}+ \dot {f}^{\prime \prime} &=0, \end{align}
where
\begin{align} \delta ^t w_{i,k}^{(1 \rightarrow \mathcal{N}, \mathcal{N}+2)} &= \rho _{i,k}- \rho _{i,k-1}- \dfrac {1}{a_{i,k}^2}\left ( \rho _{i,k}-\rho _{i,k-1}\right ), \end{align}
\begin{align} \delta ^+ w_{i,k}^{(1 \rightarrow \mathcal{N}, \mathcal{N}+2)} &= \rho _{i,k}- \rho _{i-1,k}- \dfrac {1}{a_{i,k}^2}\left ( \rho _{i,k}-\rho _{i-1,k}\right ) , \end{align}
\begin{align} \delta ^t w_{i,k}^{(\mathcal{N}+1)} &= u_{i,k}- u_{i,k-1}+ \dfrac {1}{\rho _{i,k}a_{i,k}}\left ( p_{i,k}-p_{i,k-1}\right ), \end{align}
\begin{align} \delta ^+ w_{i,k}^{(\mathcal{N}+1)} &= u_{i,k}- u_{i-1,k}+ \dfrac {1}{\rho _{i,k}a_{i,k}}\left ( p_{i,k}-p_{i-1,k}\right ). \end{align}
The backward characteristic is evaluated as
\begin{equation} \delta ^t w_{i,k}^{(\mathcal{N}+3)}+ \dfrac {\Delta t}{\Delta x}\left ( u^{\prime}_{i,k}-a_{i,k} \right )\delta ^+ w_{i,k}^{(\mathcal{N}+3)}+ \dot {f}^{\prime \prime} =0, \end{equation}
where
\begin{align} \delta ^t w_{i,k+1}^{(\mathcal{N}+3)} &= u_{i,k+1}- u_{i,k}- \dfrac {1}{\rho _{i,k}a_{i,k}}\left ( p_{i,k+1}-p_{i,k}\right ), \end{align}
\begin{align} \delta ^+ w_{i,k}^{(\mathcal{N}+3)} &= u_{i,k}- u_{i-1,k}- \dfrac {1}{\rho _{i,k}a_{i,k}}\left ( p_{i-1,k}-p_{i,k}\right ). \end{align}
Note that in (2.33), the coefficient of
$w_{i,k}^{(\mathcal{N}+3)}$
is positive and larger than one, which guarantees the stability of the scheme for the
$w^{+}$
family of waves. The grid spacing
$\Delta x$
in (2.27), (2.28) and (2.33) is constructed as
\begin{equation} \Delta x=\ell \Delta \eta ,\end{equation}
where
$\eta$
is a non-dimensional coordinate with value 0 at the shock and 1 at the contact discontinuity. In order for the model to be solvable, an additional relation is needed to determine the slug length
$\ell$
. The additional equation is obtained by imposing that at the contact discontinuity the relative velocity in the shock frame of reference matches the rate of change of the test slug length:
\begin{equation} u^{\prime}= \dfrac {{\textrm d}\ell }{{\textrm d}t}. \end{equation}
This is imposed in the scheme as
\begin{equation} \delta \eta _{k} + \eta _{k-1} - \eta _{k} + \Delta t u_{\mathcal{M}+1,k}^{\prime} = 0. \end{equation}
2.5. Boundary conditions
In this section, the use of experimental boundary conditions is discussed in detail. Firstly, the treatment of the experimentally measured shock trajectory is discussed, specifically the process of converting the experimental measurement of shock speed into a usable boundary condition. Next, the optional use of experimental pressure data is presented.
2.5.1. Shock
The experimental shock speed history is used to determine the fluid properties at the first streamwise cell at each time point. The shock speed trajectory is typically available in the form of
$x{-}t$
data where
$x$
is a vector of timing station locations and
$t$
is the arrival time of the shock at each location. In order to discretise the shock trajectory a function, typically a high-order polynomial or piecewise cubic spline, is fitted to the shock speed vector and average arrival time, i.e. a time between two station arrival times equidistant from each. The function,
$U_s (t)$
, is then evaluated at each simulation time point, allowing the flow properties at the first streamwise cell to be evaluated using the Rankine–Hugoniot relations for a thermally and chemically frozen gas. These are
\begin{align} \frac {\rho _2}{\rho _1} &= \frac {(\gamma + 1) M_{s}^2}{(\gamma -1)M_{s}^2 + 2}, \end{align}
\begin{align} \frac {u_2}{u_1} &= \frac {(\gamma -1)M_{s}^2 + 2}{(\gamma + 1) M_{s}^2}, \end{align}
\begin{align} \frac {p_2}{p_1} &= 2\gamma M_{s}^2 - \frac {\gamma - 1}{\gamma + 1}, \end{align}
where
$1$
and
$2$
represent the pre- and post-shock states,
$M_s = U_s/a_1$
represents the Mach number of the shock and
$\gamma$
is the ratio of specific heats. The post-shock species mass fractions are identical to the pre-shock values. The vibrational temperature
$T_V$
is always initialised as the fill temperature. These flow properties constitute the boundary condition at the shock.
It should be noted that, in many cases, the spatial resolution of experimental shock speed measurements is sufficiently low and/or the variation in shock speed is sufficiently high to mean that the choice of fit has an appreciable effect on the final test slug prediction. The uncertainty in test slug properties introduced because of this is discussed at length in § 5.
2.5.2. Final pressure
The second boundary condition is the streamwise pressure distribution at the final time point. Since the boundary pressure is known as a function of post-shock distance,
$p_b(x)$
, the residual can be evaluated for each Newton iteration at each streamwise point as
$r(x_i) = p(x_i) - p_b(x_i)$
, where
$p(x_i)$
is the current solution. The other primitive variables are not set explicitly by this boundary condition.
The domain length
$\eta$
is an unknown variable and thus changes between iterations. Therefore, there is not a direct correspondence between
$p(x)$
and
$p_b(x)$
and it is necessary to obtain
$p_b(x_i)$
via interpolation. In most applications, data for the pressure boundary condition are available as a function of time, i.e.
$p_b(t)$
. In these cases the interpolation is performed using the time of flight for each streamwise point,
$\varTheta = x_i /u_i$
.
It should be noted that if a wall static pressure transducer measurement is used as the simulation boundary condition there is an implicit assumption that the pressure history at the wall is identical to the pressure history at the tube centreline. This is not strictly correct because curvature of the shock can lead to a different time of arrival and pressure field distribution between the wall and centreline. The degree of curvature should be analysed on a test-to-test basis to assess the validity of the assumption using a relation such as that of De Boer (Reference De Boer1963). For the test conditions analysed in the present work, shock curvature is not expected to be of the order of 1 mm and therefore not expected to introduce significant error.
2.6. Method of solution
The method described in the previous section generates a large system of algebraic equations in
$(\mathcal{N}+3) \times (\mathcal{M} + 1) \times \mathcal{T}$
variables in as many unknowns. This system is highly nonlinear because of the reaction source terms and because of the overall mass balance for the determination of the test slug length which keeps all the equations tightly coupled. A solution can be efficiently obtained by Newton iterations and by exploiting the structure of the Jacobian matrix of the system. The system is solved using an exact Jacobian where the inverse is approximated by an incomplete factorisation method. The Jacobian is of the form
$\unicode{x1D645} \in \mathbb{R}^{(\mathcal{N} + 3) \mathcal{\times } (\mathcal{N} + 3) \mathcal{\times } \mathcal{T \times }(\mathcal{M} + 1)}$
, where
$\mathcal{N}$
is the number of species,
$\mathcal{T}$
is the number of time points and
$\mathcal{M}$
is the number of spatial points. A Schur complement is added to the Jacobian matrix that contains terms related to the determination of the non-dimensional test slug length,
$\eta$
. The Schur complement is of the form
$\unicode{x1D64E} \in \mathbb{R}^{(\mathcal{N} + 3) \mathcal{\times } \mathcal{T \times }(\mathcal{M} + 1)}$
. Given the numerical scheme presented in § 2.4, derivatives for adjacent spatial and temporal points must be stored and thus the Jacobian takes a nested tridiagonal form. This structure is represented below with
$\unicode{x1D642}$
denoting the tridiagonal matrix for the full spatial grid at a single time point:
\begin{align} \unicode{x1D642} &= \begin{bmatrix} \unicode{x1D642}_1 & \unicode{x1D642}_2 & 0 & \ldots & 0 \\ \unicode{x1D642}_{1} & \unicode{x1D642}_2 & \unicode{x1D642}_3 & \ddots & \vdots \\ 0 & \unicode{x1D642}_2 & \unicode{x1D642}_3 & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \unicode{x1D642}_{\mathcal{T}} \\ 0 & \ldots & 0 & \unicode{x1D642}_{\mathcal{T}-1} & \unicode{x1D642}_{\mathcal{T}} \end{bmatrix}, \end{align}
\begin{align} \unicode{x1D642}_k &= \begin{bmatrix} \unicode{x1D642}_{1,k} & 0 & 0 & \ldots & 0 \\ \unicode{x1D642}_{1,k} & \unicode{x1D642}_{2,k} & 0 & \ddots & \vdots \\ 0 & \unicode{x1D642}_{2,k} & \unicode{x1D642}_{3,k} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \ldots & 0 & \unicode{x1D642}_{\mathcal{M},k} & \unicode{x1D642}_{\mathcal{M}+1,k} \end{bmatrix}. \end{align}
Note that while the
$k+1$
entries are stored as they are required for the reverse time integration approach, the
$i+1$
entries are not required by the scheme and are thus left unpopulated. Here
$\unicode{x1D642}_{i,k}$
denotes the Jacobian and Schur complement entries for the
$i$
th spatial point and
$k$
th time point:
\begin{align} \unicode{x1D642}_{i,k} = \left [ \begin{array}{ccc|c} \frac {\partial F_1}{\partial q_1} & \ldots & \frac {\partial F_1}{\partial q_{\mathcal{N}+3}} & 0 \\ \vdots & \ddots & \vdots & 0 \\ \frac {\partial F_{\mathcal{N}+3}}{\partial q_1} & \ldots & \frac {\partial F_{\mathcal{N}+3}}{\partial q_{\mathcal{N}+3}} & 0 \\[12pt] \hline\\[-5pt] \frac {\partial F_\eta }{\partial q_1} & \ldots & \frac {\partial F_\eta }{\partial q_{\mathcal{N}+3}} & 1 \end{array} \right ]. \end{align}
Here,
$F$
and
$q$
represent the corrections to the characteristic variables and the vector and primitive variables, respectively. The structure of the linear system is then
\begin{equation} \unicode{x1D645} \unicode{x1D64C} = \unicode{x1D64D} ,\end{equation}
where
$\unicode{x1D64C}$
and
$\unicode{x1D64D}$
are the correction and residual vectors. Vector
$\unicode{x1D64D}$
is stored in characteristic form, whereas
$\unicode{x1D64C}$
is stored in primitive form. If the system is partitioned as follows:
\begin{equation} \unicode{x1D645} = \begin{bmatrix} \dfrac {\partial \unicode{x1D64D}_B}{\partial \unicode{x1D64C}_B} & \dfrac {\partial \unicode{x1D64D}_B}{\partial \unicode{x1D64C}_T}\\[10pt] \dfrac {\partial \unicode{x1D64D}_T}{\partial \unicode{x1D64C}_B} & 1\\ \end{bmatrix} = \begin{bmatrix} \unicode{x1D645}_{BB} & \unicode{x1D645}_{BT}\\ \unicode{x1D645}_{TB} & 1 \end{bmatrix}, \; \unicode{x1D64C} = \ \begin{bmatrix} \unicode{x1D64C}_B\\ \unicode{x1D64C}_T \end{bmatrix}, \; \unicode{x1D64D} = \ \begin{bmatrix} \unicode{x1D64D}_B\\ \unicode{x1D64D}_T \end{bmatrix}, \end{equation}
where
\begin{equation} \unicode{x1D64C}_B, \unicode{x1D64D}_B \in \mathbb{R}^{\mathcal{T} \times (\mathcal{N} + 3) \times (\mathcal{M} + 1)}, \; \unicode{x1D64C}_T, \unicode{x1D64D}_T \in \mathbb{R}^{\mathcal{T \times }(\mathcal{M} + 1)}, \end{equation}
then
$\unicode{x1D64C}$
can be obtained via
\begin{align} \unicode{x1D64C}_B \to { \unicode{x1D64C}{\kern0.5pt}'}_B & = \unicode{x1D64C}_B + \unicode{x1D645}_{BB}^{-1} ( \unicode{x1D64D}_B - \unicode{x1D645}_{BB} \unicode{x1D64C}_B - \unicode{x1D63C}_{BT} \unicode{x1D64C}_{T}), \end{align}
\begin{align} \unicode{x1D64C}_T \to { \unicode{x1D64C}{\kern0.5pt}'}_T &= \unicode{x1D64C}_T + \omega ( \unicode{x1D64D}_T - \unicode{x1D645}_{TB} \unicode{x1D64C}_B - \unicode{x1D64C}_{T}), \end{align}
where
$\tilde { \unicode{x1D645}}^{-1}$
is an approximation of
$J^{-1}$
constructed via an incomplete factorisation method, with
$\omega \lt 1$
a relaxation factor. The solution vector is obtained via Newton iterations with
\begin{equation} \unicode{x1D64C} \to \unicode{x1D64C} - \tilde { \unicode{x1D645}}^{-1} \unicode{x1D64D}. \end{equation}
3. Validation
The methodology described in the previous section has been implemented as part of a software package called LAgrangian Shock Tube Analysis 2.0 (LASTA 2.0). Results obtained using LASTA 2.0 are presented here. The results presented were selected to achieve the following aims:
-
(i) Validate the implementation of the quasi-1-D Euler equations with boundary-layer source terms by predicting the properties of an ideal test gas in the presence of shock speed non-uniformity.
-
(ii) Validate the implementation of the two-temperature non-equilibrium model.
These aims have been addressed by firstly comparing predictions made using the method with results from a viscous, axisymmetric Navier–Stokes solver for several shock trajectories in argon. Then, the two-temperature model implementation was validated by comparison to an established quasi-1-D method. All simulations were run on a workstation with an Intel(R) Xeon(R) W-2265 CPU clocked at 3.50 GHz and 64.0 GB of RAM. Comparisons with experiment are presented in the next section.
3.1. Fluid model
The method was validated against a series of numerical simulations by Satchell et al. (Reference Satchell, Di Mare and McGilvray2022a ) using the framework for overset simulation of shock tubes (FROSST). These validation cases use an ideal gas model and thus remove the complication of thermochemistry. The FROSST method is an axisymmetric, explicit, structured finite-volume method for solution of the Navier–Stokes equations for ideal gas simulation of shock tube flows. The code is written using an overset formulation that allows a highly refined grid to be attached to both shock and contact discontinuity with the intention of reducing computational cost when compared to typical axisymmetric shock tube solvers. The FROSST method has the capability to model the effect of primary diaphragm opening, allowing for the production of non-uniform shock trajectories. It has been validated against a number of historical and canonical shock tube cases.
Satchell et al. (Reference Satchell, Di Mare and McGilvray2022a
) devised three ideal gas shock tube test cases with the same final shock speed but differing initial shock speeds by varying the driver behaviour, thus allowing the effect of non-uniform shock speed on the test slug properties to be evaluated. The trajectories are shown in figure 6, consisting of a decelerating, accelerating and relatively constant shock speed case. All cases result in a shock speed of
$2100$
$\textrm{m s}^{-1}$
at a distance
$8$
m from the primary diaphragm. Pure argon was used as the test gas at an initial pressure and temperature of
$66.7$
Pa and
$300$
K, respectively. A tube diameter of
$100$
mm was selected. The FROSST simulations were viscous.
Figure 6. Shock speed as a function of distance from the primary diaphragm for the three test cases.
For the new method, the domain was discretised into 150 streamwise points behind the shock with 10 % grid clustering towards the shock front, that is,
$x(i)=x(i-1) + {\textrm d}x(i)$
, where
${\textrm d}x(i) = 1.05 {\textrm d}x(i-1)$
with
${\textrm d}x(1)=0.005$
m. For each simulation, the full shock trajectory and the pressure profile at
$9$
m from the primary diaphragm from the FROSST test cases were used as boundary conditions. The supplied shock trajectory was discretised into 51 time points. The pressure profile at any location along the length of the tube may in principle be used as the boundary condition. However,
$9$
m was selected so that the profiles at
$8$
m, the location where the test case shock trajectories converge on
$2100$
$\textrm{m s}^{-1}$
, were within the domain of influence of but not explicitly set by the new boundary condition. This arrangement is analogous to interrogating a location slightly upstream of a pressure transducer in an experimental facility, which is ultimately the intended use case when modelling shock layer radiation.
The full test slug solution for the decelerating (2800
$\textrm{m s}^{-1}$
) case is shown in figure 7(
a), where the
$z$
axis is static pressure. The boundary conditions imposed on the solver can be clearly seen; the decreasing pressure with time at 0 m post-shock is the shock boundary condition and occurs as a result of the decelerating shock trajectory. The post-shock pressure profile at the final time point is the pressure boundary condition and is taken from the FROSST simulation. The paths of particles and characteristic waves in the test slug have also been shown. These were determined by selecting uniformly distributed seed points at the shock and pressure boundary conditions and then numerically integrating the positions based on the test slug solution. The paths are seen to cluster together towards the end of the test slug, showing the effect of compression due to growth of the boundary layer. Figure 7(
b) shows, for the same particle paths, the local temperature as a function of time of flight. Here,
$t = 0$
at the time the particle enters the test slug through the shock. This clearly shows that individual particles experience different fluid properties throughout their histories. This fact has implications for modelling of non-equilibrium behaviour, which are discussed in more detail later in the paper.
Figure 7. Full test slug solution for the decelerating (2800
$\textrm{m s}^{-1}$
) case. The surface colouring is proportional to the static pressure. (a) Static pressure contour. (b) Particle histories.
Figure 8. (a–f) Profiles of pressure and temperature obtained at three axial locations along the shock tube for the two-dimensional axisymmetric simulations and the present work.
Comparisons of static pressure and temperature between the new method and the FROSST solver are shown in figure 8. Profiles at
$8$
,
$8.5$
and
$9$
m from the primary diaphragm are shown. The pressure profiles exactly match at
$9$
m from the diaphragm, the location of the boundary condition. Pressure profiles predicted by the new method at
$8$
and
$8.5$
m also show excellent agreement with FROSST where the difference is less than 1 % for all cases.
Profiles of temperature, which are not set explicitly as a boundary condition in the new method, agree with FROSST within 4 % for the first 200 μs. The temperature predicted by the method for the decelerating case begins to diverge towards the rear of the test slug due to the arrival of a cold jet of driver gas in the FROSST simulation. This is expected, since the three-dimensional nature of the contact discontinuity cannot be fully described by a quasi-1-D method. The difference is in part due to boundary-layer implementation. Contour plots of temperature in the FROSST simulations are shown in Satchell et al. (Reference Satchell, Di Mare and McGilvray2022a ), illustrating the arrival of the jet.
These comparisons demonstrated that the new method provides satisfying predictions of the flow in the core flow of a shock tube in the presence of shock speed non-uniformity at low computational cost. They have shown that additional experimental information, namely a wall static pressure history, can be included via a reverse time-integration method.
3.2. Two-temperature non-equilibrium thermochemistry model
To validate the implementation of the non-equilibrium thermochemistry model, the method was compared with the non-equilibrium shock tube solver or NESS (Clarke et al. Reference Clarke, Collen, McGilvray and Di Mare2023a ). Both NESS and the new method utilise an in-house thermochemistry library known as OCEAN. The OCEAN library provides functions to evaluate thermodynamic properties, transport properties and reaction and relaxation rates in non-equilibrium two-temperature plasmas. The parameters of interest and their exact Jacobians with respect to a set of working variables are returned. In this section, the NESS solver and test case are briefly outlined and results are discussed.
The NESS solves the compressible, steady Navier–Stokes equations in quasi-1-D form to establish the flow state behind the normal shock wave in a shock tube. The NESS utilises a Park two-temperature approach to model the thermochemical relaxation behind the shock wave. The NESS is the first tool to combine the solution of the viscous Navier–Stokes equations with appropriate source terms for mass loss from the core flow to the boundary layer; however, NESS is used here in its most basic form with no boundary-layer source terms and initialised with the Rankine–Hugoniot post-shock properties. The NESS has been validated against the POSHAX 3 solver (Gollan Reference Gollan2007) and a variety of experimental data (Clarke et al. Reference Clarke, Brody, Steer, McGilvray and Di Mare2024).
The validation case chosen was a 4.72
$\textrm{km s}^{-1}$
shock in a 133.3 Pa synthetic air mixture (80.2 %
$\textrm{N}_{2}$
/19.8 %
$\textrm{O}_{2}$
by volume). Both solvers were initialised with the Rankine–Hugoniot post-shock properties. The NESS was run using Rankine–Hugoniot post-shock properties for the purposes of validation only. The chemical reaction scheme of Marrone & Treanor (Reference Marrone and Treanor1963) was used. The new method was run with 200 streamwise grid points compared with 4000 for the NESS simulation. The shock speed was constant. No boundary-layer effects were included for either solver.
Results from the validation simulations are shown in figure 9. Excellent agreement was generally achieved between the methods and differences closer to the shock are attributed to the viscous effects modelled by NESS that are not included in the present work.
Figure 9. (a–d) Comparison of non-equilibrium fluid properties between the method of Clarke et al. (Reference Clarke, Brody, Steer, McGilvray and Di Mare2024) and the present work for a 4.27
$\textrm{km s}^{-1}$
shock in 133.3 Pa of 80.2 %
$\textrm{N}_{2}$
/19.8 %
$\textrm{O}_{2}$
(v/v).
3.3. Grid study
A grid refinement study was also performed to assess the effect of the number of post-shock cells on the non-equilibrium region. Here, the simulation settings were maintained from the previous section but the shock speed and fill pressure were 5.3
$\textrm{km s}^{-1}$
and 26.6 Pa, respectively. The number of streamwise grid points were increased in subsequent simulations. The resulting vibrational–electronic temperature profiles and percentage difference from the highest-grid-density case are shown in figure 10. The maximum temperature difference between the cases and the finest mesh (n = 500) was 307, 23, 7 and 3 K. This corresponds to a convergence of within 1.5 % for the 200-point case, which is deemed as a satisfactory compromise between computational time and accuracy and is used as a minimum grid size for the remainder of the document.
Figure 10. Dependence of the predicted vibrational-electronic temperature on the number of streamwise points chosen: (a)
$T_V$
; (b) percentage difference in
$T_V$
.
The slug length should be taken into consideration when determining an appropriate grid size. The grid for the new method is built in non-dimensional units and the physical size of the grid changes depending on the slug length at a given time step. The slug length for this refinement study was 406 mm; thus a minimum grid density of 0.5 points per millimetre is required. In practice, the minimum grid size is dependent on the thermochemical model and grid independence should be confirmed on a case-by-case basis.
A refinement study considering the number of time steps chosen was also performed. Here, it was necessary to interrogate a non-uniform shock trajectory so that time-dependent differences in the solution appeared. A fictitious shock trajectory was devised for this purpose. The trajectory decelerates by 318
$\textrm{m s}^{-1}$
over 5.1 m and follows a cubic polynomial trend, ending at 5.3
$\textrm{km s}^{-1}$
. The simulations were run with 220 streamwise grid points and the number of timewise points was varied between 15 and 95. Boundary-layer effects were included in the simulations so that an assessment of the assumption of instantaneously boundary-layer reshaping discussed in § 2.2 could be made. It would be expected that, if this was a poor assumption, changing the number of time steps would have a large effect on the solution because the incorrect boundary-layer shape had been used.
Figure 11. Dependence of the predicted vibrational–electronic temperature on the number of time steps chosen: (a)
$T_V$
; (b) percentage difference in
$T_V$
.
The predicted vibrational–electronic temperatures for four timewise grid sizes (15, 30, 60 and 95) are shown in figure 11(
a). The percentage difference between each result and the prediction for the highest grid density (
$t=95$
) is shown in figure 11(
b). The
$t=60$
and
$t=90$
cases converge to within 0.5 % which is deemed acceptable for subsequent simulations. This also demonstrates that, at least for this test case, the instantaneous boundary-layer reshaping assumption is acceptable.
4. Application
Given satisfactory validation of the underlying fluid mechanics and two-temperature thermochemistry models, the new method is now used to model non-equilibrium thermochemistry in real shock tubes. Direct comparisons between the numerical model and optical emission spectroscopy measurements have been made. This required the use of a radiation model, which is discussed first. Results from simulations of test conditions relevant to lunar return and Titan entry are presented from two experimental facilities.
4.1. Facilities overview
Optical emission spectroscopy data from the NASA EAST and Oxford T6 Aluminium Shock Tube (AST) facilities are modelled in this section. For the purposes of this study, the key differences between the two are:
-
(i) The internal diameter: 101.6 mm for the EAST cases and 225 mm for the T6 AST cases.
-
(ii) The driver technology: an electric arc for EAST and a free piston for T6.
These differences are highlighted as they may be expected to change the relative importance of certain phenomena taking place in the test gas. With regard to the diameter, as shown in (2.17) and (2.18), the radial velocity at the boundary-layer edge in the Mirels model is inversely proportional to the tube diameter and thus viscous effects are expected to be more important in the EAST cases. The differing driver technologies are likely to introduce different upstream wave patterns into the test gas, leading to different test slug properties. For example, the free piston continues to compress the driver gas after diaphragm rupture, potentially lessening the magnitude of reflected expansion (
$u+a$
) waves reaching the test slug when compared with the arc driver where there is no further compression.
Aside from these differences, the operation of the two facilities for the test cases presented was nearly identical. The reader is directed to Collen et al. (Reference Collen, Doherty, Subiah, Sopek, Gildfind, Penty, Gollan and Morgan2021) and Cruden (Reference Cruden2014) for a comprehensive overview of the operation of T6 and EAST, respectively.
4.2. Radiation simulations
Radiation simulations in the present work were run with the Nonequilibrium Radiative Transport and Spectra (NEQAIR) program, developed by Park (Reference Park1985). NEQAIR is a method for simulating radiative transport and benchmark approach for comparison with shock tube data. The full capabilities of the solver are summarised by Cruden & Brandis (Reference Cruden and Brandis2014). NEQAIR takes flow-field properties (temperature, species number densities) as inputs and outputs spectral radiation characteristics, including emission and absorption spectra along a given line of sight.
All NEQAIR simulations in this work were run with the recommended default settings. NEQAIR version 15.2.2 was used. For the non-Boltzmann simulations, the local escape flux limited method with a characteristic length of 1 cm was selected as recommended by Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022). The shock tube mode was selected with species number densities and position imported directly from the fluid simulations. Boundary conditions were set with no initial radiance and blackbody at the final Line of sight point temperature so that incident radiant heat flux was returned. The simulation wavelength range was set depending on the experimental data with which to be compared; the spectral range has been reported alongside the results. All bands were included in the simulations. Automatic grid spacing was selected with 10 points per line and a range of 10. Spatial scanning with convolution functions provided with the experimental data was performed. NEQAIR input files used for radiation simulations are available on request.
We note the improved non-Boltzmann modelling implemented in NEQAIR 15.2.2 and its importance for interpreting the results of the present work in comparison with previous work. In particular, these changes greatly improve the CN modelling when compared with experiment and change the predicted non-Boltzmann radiation in the 760–785 nm region when compared with NEQAIR 15.0 used in Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022, Reference Collen, Di Mare, McGilvray and Satchell2023).
4.3. A 9.4
$\textrm{km s}^{-1}$
shock in 26.6 Pa of
$\textrm{N}_2$
/
$\textrm{O}_2$
for lunar return
This section concerns modelling of a shock tube experiment at lunar-return-relevant conditions. The experiment was performed in the Oxford T6 Stalker Tunnel AST. The diaphragm rupture conditions and piston dynamics were selected such that a non-uniform shock trajectory was produced; the experimental shock speed measurements are shown as points in figure 12. The effects of shock trajectory for this test were studied by Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022) using Satchell’s method and predictions of the post-shock equilibrium radiance in the 760–785 nm region were observed when compared with an assumption of constant shock speed. However, Satchell’s method was limited by the assumption of thermochemical equilibrium and monotonically increasing pressure. Therefore, predictions of non-equilibrium behaviour were not made and predictions of the post-shock pressure were non-physical. This test case is then an opportunity to benchmark the new method against the current state-of-the-art approach.
The fit that was applied to the experimental shock speeds for simulations with the new method is shown in figure 12. The pressure history from a transducer located 9.58 m from the primary diaphragm is also shown, this measurement being used to constrain the pressure at the end of the domain. A 2 μs moving average filter (shown) was applied to the raw experimental trace before it was used in the simulation to improve solution stability. An 11-species gas model (
$\textrm{e}^-$
, N, O,
$\textrm{N}^+$
,
$\textrm{O}^+$
,
$\textrm{N}_2$
, NO,
$\textrm{O}_2$
,
$\textrm{NO}^+$
,
$\textrm{N}_2^+$
,
$\textrm{O}_2^+$
) was used. The reaction scheme was that of Park (Reference Park1993). The Millikan–White approach to translational–vibrational relaxation times is utilised (Millikan & White Reference Millikan and White1963). The final grid consisted of 800 streamwise points and 62 time points, taking 0.947 CPU hours to solve.
Figure 12. Shock speed as a function of distance from the primary diaphragm for the lunar return test case, T6s132. The experimental pressure history used as a boundary condition is also shown with the fit applied overlaid. This transducer was a PCB113B27 model located just downstream of the optical emission spectrometer at 9.58 m.
4.3.1. On the choice of state population method
The radiance emitted in the 760–785 nm region is dominated by the 777 nm O triplet (
$\textrm{3p}^{5}P_{1-3}$
states) and the
$\textrm{N}_2$
first positive system (
$\textrm{B}^{3}$
$\varPi$
g–
$\textrm{A}^{3}$
$\varSigma$
u+). The choice of state population method in NEQAIR (i.e. non-Boltzmann versus Boltzmann) changes the population densities of these states which in turn has a significant effect on the predicted radiance in this region. We now examine the radiance predictions for each case to determine the most appropriate choice for use with the new method in the context of interrogating non-ideal flow effects.
We firstly consider a test case with the new method where the shock trajectory is as shown in figure 12, the final pressure is constrained using measurements from the facility at the observation location as in figure 12 and boundary-layer effects are included. Figure 13 shows a comparison between the non-Boltzmann population method with settings as described in § 4.2, a Boltzmann population method and the experiment. All other NEQAIR settings and inputs were held constant. The state energy distribution within O atoms was likely Boltzmann after 18 mm as an improved agreement with experiment is offered until 60 mm, where the results converge. Given the temperatures are identical between cases (figure 13
b), this is attributed to the reduced population of the
$\textrm{3p}^{5}P_{1-3}$
states, shown in figure 13
c).
The Boltzmann results show reduced agreement against experiment when compared with the non-Boltzmann in the non-equilibrium region (i.e. between 0 and 18 mm post-shock). However, we acknowledge that disagreement in the non-equilibrium region is generally to be expected for many reasons. These include but are not limited to:
-
(i) The assumption of Rankine–Hugoniot post-shock properties, rather than modelling the shock with finite thickness and with diffusion terms as done by Clarke et al. (Reference Clarke, Brody, Steer, McGilvray and Di Mare2024); particularly at high shock speed, the Rankine–Hugoniot method will result in high translational temperatures (
$\gt 40\,000$
K) that are non-physical. -
(ii) The non-equilibrium reaction and relaxation rates used in the fluid and radiation simulations; these are fundamentally the subject of study and are associated with high uncertainty, as shown in figure 5.
-
(i) The two-temperature, 11-species gas model; the reaction processes taking place in the gas are state-dependent and thus the model used here is an approximation.
The contribution of the present work is to quantify uncertainty in the shock trajectory as a source of error. The method described herein serves as a platform to optimise non-equilibrium models against shock tube experiment with the non-uniform experimental effects accounted for. We therefore proceed with an assumption of Boltzmann state distributions in the following air radiation simulations to demonstrate the capabilities of the method in the equilibrium region where there is highest confidence in the fluid properties and leave tuning of the non-equilibrium and non-Boltzmann rate parameters against experiment for future work.
Figure 13. Comparison of (a) radiance, (b) temperature and (c) O
$\textrm{3p}^{5}P_{1-3}$
number density when using a non-Boltzmann and Boltzmann state population method in NEQAIR. The experiment is shown in black and the CEA equilibrium values for the final shock speed of 9.40
$\textrm{km s}^{-1}$
are shown as dashed lines. The shaded region represents perturbations of
$\pm$
100
$\textrm{m s}^{-1}$
to the nominal shock trajectory.
4.3.2. Effect of the boundary conditions
To show how each feature of the method contributes to changes in predictions of radiance, three cases are presented:
-
(i) Constant shock speed of 9.40
$\textrm{km s}^{-1}$
, unconstrained pressure, boundary-layer effects neglected. -
(ii) Shock trajectory from figure 12, unconstrained pressure, boundary-layer effects included.
-
(iii) Shock trajectory from figure 12, pressure constrained using measurements from the facility at the observation location, boundary-layer effects included.
Figure 14(a) shows the predicted radiance in the 760–785 nm region using NEQAIR 15.2.2, the experimental result and the difference between the simulations and experiment. Figure 14(b,c) shows the predicted temperatures and O
$\textrm{3p}^{5}P_{1-3}$
state number densities (corresponding to the upper levels of the 777 nm triplet). The CEA equilibrium values for the final shock speed have been included. Significant differences in temperature and number density, and hence radiance, were observed between the simulation approaches in the equilibrium region. Inclusion of the decelerating shock trajectory as a boundary condition gives rise to an increasing radiance trend in the equilibrium region when compared with the constant-shock-speed case, a result of the gradient in temperature and enthalpy, which is consistent with the findings of Collen et al. (Reference Collen, Di Mare, McGilvray and Satchell2023).
Figure 14. Comparison of (a) radiance, (b) temperature and (c) O
$\textrm{3p}^{5}P_{1-3}$
number density with the new method run with various combinations of boundary conditions. The experiment is shown in black and the CEA equilibrium values for the final shock speed of 9.40
$\textrm{km s}^{-1}$
are shown as dashed lines. The shaded region represents perturbations of
$\pm$
100
$\textrm{m s}^{-1}$
to the nominal shock trajectory.
Although qualitative agreement with experiment is improved when shock trajectory effects are included, the magnitude of radiance is increasingly over-predicted towards the rear of the slug. The radiance at 100 mm post-shock is over-predicted by 0.35
$\textrm{W sr}^{-1}\,cm^{-2}$
at 100 mm post-shock when the pressure is not constrained, compared with an under-prediction of 0.12
$\textrm{W sr}^{-1}\,cm^{-2}$
for the constant-shock-speed case. This is in part caused because, as discussed previously,
$u+a$
waves from the driver have not been appropriately modelled.
The static pressures at the end of the simulation domain (9.5 m from the primary diaphragm) are shown in figure 15 and provide an illustration of the shortcomings of using only the shock trajectory to model the test slug. When the slug properties are determined using the shock trajectory alone, the predicted pressures disagree with experiment. Differences in pressure give rise to different thermochemical states and thus different radiance predictions. Thus, when the simulation pressure is constrained to the experimental values, improved agreement with the experimental radiance was observed whilst preserving the increasing trend due to the shock trajectory.
Figure 15. Pressure predicted for the lunar return case. Experimental pressure measurements were obtained using a PCB113B27 static pressure transducer located 9.58 m from the primary diaphragm.
Despite constraining the model with two experimental measurements, appreciable differences from the experimental radiance remain. Many factors are at play, including the radiation model (discussed in the next section) and thermochemistry model, but the quality of the experimental measurements used as boundary conditions is of most immediate relevance to the present work. The shaded region in figure 14(a) shows the solution bounds when the nominal shock speed is perturbed by
$\pm$
100
$\textrm{m s}^{-1}$
. Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022) report shock speed uncertainties of between 50 and 500
$\textrm{m s}^{-1}$
for this test condition and so 100
$\textrm{m s}^{-1}$
perturbations were selected as a representative, order-of-magnitude estimate. The effect of shock speed uncertainty is discussed in greater detail in § 5. The shock speed perturbations were found, for this condition, to strongly affect the radiance predictions in the equilibrium region, up to
$\pm$
18 % at 100 mm post-shock. Therefore, without improved shock speed measurements, for example using the microwave reflectometry technique developed by Dufrene, Holden & Ringuette (Reference Dufrene, Holden and Ringuette2015), it cannot be said with certainty that the true shock trajectory was modelled.
Some comments concerning the quality of the experimental pressure measurements have already been made, namely that measurements obtained at the wall are not expected to be truly representative of the core flow. We add that the rise time of the PCB113B27 pressure sensor used in this case is quoted as
$\leqslant 1$
μs or potentially greater than 30 % of the test time. We also note that information from the experimental pressure boundary condition may only propagate through a subset of the solution domain. Figure 16 shows pressure distribution across the full domain for this test case. Here, the shock trajectory boundary condition and pressure boundary condition are clearly seen to dictate the slug properties at the final time step, as well as the domain of influence for the pressure boundary condition. The figure illustrates that
$u+a$
waves are not experimentally constrained for the majority of the solution domain.
Figure 16. Pressure distribution through the full test slug for the lunar return case.
4.3.3. Comparison with Satchell’s method
Figure 17 shows a comparison between Satchell’s method and the present work. For this case, Satchell’s method yields a generally higher radiance prediction than the present work when the pressure is constrained to the experiment. This occurs because of the reduced temperature and O
$\textrm{3p}^{5}P_{1-3}$
number density associated with the decreased pressure. This behaviour can be recovered to a certain extent applying an equilibrium boundary condition at the shock and constraining the final pressure to that predicted by Satchell’s method. In this case, the temperature, O number density and thus radiance predictions are essentially identical to those of Satchell’s method. The slight discrepancies are attributed to different approaches to trajectory discretisation between the present work and Satchell’s method.
Figure 17. Comparison of (a) radiance, (b) temperature and (c) O number density with the new method run with various combinations of boundary conditions and Satchell’s method. The experiment is shown in black and the CEA equilibrium values for the final shock speed of 9.40
$\textrm{km s}^{-1}$
are shown as dashed lines. The shaded region represents perturbations of
$\pm$
100
$\textrm{m s}^{-1}$
to the nominal shock trajectory. No spatial convolution was applied for Satchell’s method.
Pressure profiles for the experiment and the simulation results presented in figure 17 are shown in figure 18. Again, the shortcomings of modelling only the shock trajectory are evident as Satchell’s method over-predicts the experimentally observed static pressure.
Figure 18. Pressure predicted for the lunar return case. Experimental pressure measurements were obtained using a PCB113B27 static pressure transducer located 9.58 m from the primary diaphragm.
4.4. A 4.8
$\textrm{km s}^{-1}$
shock in 85 Pa of
$\textrm{N}_2$
/
$\textrm{CH}_4$
for Titan entry
A study of the radiation in a Titan-relevant atmosphere is presented in this section. The post-shock gas in these mixtures is often in a state of non-equilibrium due to the high concentration of
$\textrm{N}_2$
, the presence of
$\textrm{CH}_4$
and the resulting formation of CN. The radiation emitted by such mixtures has been the subject of study in shock tubes for many years to provide support and post-flight analysis of the Huygens entry and more recently the design of the Dragonfly heat shield.
Recent shock tube studies of this mixture have highlighted the need to include non-uniform shock speed effects in the modelling approach (Tibère-Inglesse et al. Reference Tibère-Inglesse, Johnston, Brandis and Cruden2023). Sopek et al. (Reference Sopek, Glenn, Clarke, Di Mare, Collen and McGilvray2024) made efforts to do this by modelling a series of experiments performed in the T6 facility using both Satchell’s method and the time-of-flight correction proposed by Clarke et al. (Reference Clarke, Di Mare and McGilvray2023b ) for viscous shock tubes. Although shock speed variation was clearly present in the experiments analysed, the predictions of Satchell’s method were limited by the assumption of equilibrium thermochemistry. The test gas was in a strong state of non-equilibrium in most cases for the majority of the test slug and thus the numerical results greatly under-predicted the experiment. Such experiments thus constitute a worthy test case for the new method described in this paper.
A test case from a recent experimental campaign in the EAST facility was selected to benchmark the new method in a strongly non-equilibrium environment. The test case was nominally 4.82
$\textrm{km s}^{-1}$
in 85 Pa of 97.8 %
$\textrm{N}_2$
/2.2 %
$\textrm{CH}_4$
by volume. Figure 19 shows the experimentally measured shock speeds and fitted curve used in the simulations. A 19-species gas model (
$\textrm{e}^-$
, N, C, H,
$\textrm{N}^+$
,
$\textrm{C}^+$
,
$\textrm{H}^+$
,
$\textrm{N}_2$
, CN,
$\textrm{CH}_4$
,
$\textrm{CH}_3$
,
$\textrm{CH}_2$
, CH,
$\textrm{C}_2$
,
$\textrm{H}_2$
, NH, HCN,
$\textrm{N}_2^+$
,
$\textrm{CN}^+$
) was used for the Titan entry simulations with the Gokcen (Reference Gokcen2007) reaction scheme. No pressure boundary condition was applied for this test case as a pressure history of sufficient quality located at the optical emission spectrometer was not available. This should not substantially affect the findings in this section as the effect of shock trajectory is likely to be dominant due to the lower shock velocity.
Figure 19. Shock speed as a function of distance from the primary diaphragm for the Titan entry test case, E65-45.
The effect of non-uniformity in the experimental shock trajectory can be isolated by comparing the results of two simulations: one with the experimental trajectory as a boundary condition and one using a constant shock speed equal to the final shock speed. This comparison is made in figure 20 for the E65-45 test case, selected as it features the strongest shock speed non-uniformity of the Titan-relevant cases analysed. Figure 20 shows the radiance, the temperatures and the CN number density. Although the two cases feature essentially identical temperatures at the final time point, the experimental case features higher temperatures earlier in the trajectory due to the increased shock speed, leading to a reduction in CN number density that is reflected in the radiance profile.
Predictions made using Satchell’s method are also shown. Given the assumption of equilibrium thermochemistry, this method is clearly not directly applicable to shock tube cases with strong thermochemical non-equilibrium such as this one.
Figure 20. Comparison of (a) radiance, (b) temperature and (c) CN number density for E65-45 using the new method with the experimental shock trajectory and a constant value.
5. Discussion
Given that the boundary conditions for the numerical simulations are derived for experiment, they have an uncertainty associated with them. Experimental sources of uncertainty that propagate to the numerical simulation include the shock trajectory (see Appendix B), the static pressure measurements and the radiation measurements (Cruden Reference Cruden2014). The thermochemical model contains uncertainties in the reaction rate constants and relaxation parameters (Palmer et al. Reference Palmer, Prabhu and Cruden2014).
An important advantage of the method demonstrated in the present work is the ability to probe, in isolation, the effect of the experimental shock trajectory and the thermochemical model on the predicted final radiance profile; such an analysis is presented in this section for a single Titan entry test case published by Sopek et al. (Reference Sopek, Glenn, Clarke, Di Mare, Collen and McGilvray2024), namely T6s234. The effect of shock trajectory uncertainty on the final test slug state is determined by simulating a test case using trajectories determined by the nominal, upper and lower bounds of experimental shock speed uncertainty. Then, the effect of uncertainty in the reaction rate constants on the test slug properties is determined by simulating the same case with the upper and lower bounds for two rates of interest. These results are finally used to demonstrate how trajectory effects may be misinterpreted as a need to alter reaction rate constants.
5.1. Trajectory study
The effects of uncertainty in the experimentally measured shock trajectory are discussed in this section. Shock speed uncertainty was determined using the method in Appendix B, which is identical to the method in Appendix G of James (Reference James2018). Once the experimental uncertainties are known, a fit can be applied to them in the same way as was done for the nominal shock trajectory. Figure 21 shows the fits used for the analysis in this section.
Figure 21. Shock speed as a function of distance from the primary diaphragm for the three test cases.
As with many impulse facilities, the density of time of arrival measurements was higher closer to the observation window for T6s234. This was accompanied by an increase in the standard deviation of the shock speed measurements. This is because, as discussed in Appendix B, the uncertainty in
$\Delta x$
and
$\Delta t$
does not depend on the sensor spacing (
$\Delta x$
) and thus the uncertainty in the shock speed measurement (
$\delta V_s$
) is inversely proportional to
$\Delta x$
and
$\Delta t$
. Put simply, decreasing sensor spacing inherently increases shock speed uncertainty. It is common practice to use the mean of the final
$n$
measurements as the nominal shock speed in this scenario to decrease the effective uncertainty. This can often differ appreciably from the final shock speed measurement and can affect the predicted test slug properties as shown in § 6 of Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022). Ideally, a method that incorporates uncertainty history would be used, which is done here.
The resulting radiance predictions for the 240–440 nm region are shown in figure 22(a) with figure 22(b) showing the difference between the three simulations and the experimental radiance and figure 22(c) showing the point of origin in the tube for each point in the test slug; the horizontal axis corresponds to the horizontal axis in figure 21, allowing for the uncertainty in the shock speed to be correlated with post-shock distance. The peak upper and lower radiance values bound the experimental peak radiance, with the error with respect to the experiment becoming essentially the same at 15 mm post-shock.
Figure 22. (a–c) Predictions of 240–440 nm radiance for T6s234 using the nominal, lower and upper uncertainty bounds for shock speed history from figure 21.
These results show where efforts should be directed for improvement of the various physical models involved in the predictions. Differences between the simulation and experiment cannot be attributed to uncertainty in the shock trajectory after 15 mm as the experimental radiance is no longer bounded at this point. This suggests that it is valid to seek improved agreement in this region by modification of other parameters, such as the reaction rate constants, alone. Similarly, if improvement in the region between 0 and 15 mm is sought, more precise measurement of the shock speed is required to decrease uncertainty.
The
$\textrm{N}_2$
and CN number densities associated with the simulations are shown in figure 23. It can be seen that uncertainty in the shock trajectory results in uncertainty of up to
$\pm$
3.6 % and over
$\pm$
100 % in the
$\textrm{N}_2$
and CN number densities, respectively.
Figure 23. (a,c) The CN and (b,d)
$\textrm{N}_2$
species number densities predicted for the three shock trajectories shown in figure 21.
5.2. Rate constant study
Uncertainty in the numerical predictions of radiance also arises from uncertainty in the gas model, in particular the reaction rate constants themselves. Gokcen (Reference Gokcen2007) provides uncertainty bounds in table 3 for the reactions used in the present work as a variable
$F$
. The variable
$F$
is defined such that the range of values of
$k_f$
is bounded by multiplication and division of
$k_f$
by a factor
$F$
. In this section the effect of uncertainty in two reaction rates, namely C +
$\textrm{N}_2$
$\longleftrightarrow$
CN + N and
$\textrm{N}_2$
+ M
$\longleftrightarrow$
2N + M, on the predicted post-shock radiance for the nominal shock trajectory is considered;
$F$
values of 3.0 and 2.0 are quoted for these rates, respectively. These reactions were selected as they comprise the main CN formation mechanism at the relevant conditions and thus most strongly affect the radiance observed. Simulations were performed by perturbing only one rate at a time to study the effects in isolation. All other simulation parameters were held constant.
Figure 24 shows the numerical radiance predictions for the rate uncertainty study. Modification of both rates has a large effect on the predicted radiance. This demonstrates why tuning of rate constants to match experimental data is an attractive choice to improve agreement between numerical models and experimental radiance measurements. The results show that, at least for this particular test case, uncertainty due to the shock trajectory contributes equally to the value of peak non-equilibrium radiance, with uncertainties in reaction rate constants becoming larger drivers further from the shock.
Figure 24. Effective uncertainty bounds on radiance predictions using the new method for the shock trajectory, C +
$\textrm{N}_2$
$\longleftrightarrow$
CN + N reaction rate and
$\textrm{N}_2$
+ M
$\longleftrightarrow$
2N + M reaction rate.
5.3. On modelling the shock tube flow as quasi-one-dimensional
Although the quasi-1-D framework adopted in the present work enables efficient reconstruction of the test slug thermochemical state, it necessarily omits several flow features known to arise in real shock tube environments. In practice, the boundary layer develops in a fully three-dimensional manner and may intrude into the nominally uniform core flow, thereby altering the local pressure, temperature and species concentrations in ways that are not captured by a 1-D representation. Moreover, the diaphragm rupture process and geometric discontinuities at the diaphragm station can generate transverse shock structures, such as non-planar fronts and cross-sectional shocklets, that introduce spatial fluctuations in thermodynamic and chemical quantities immediately behind the primary shock. These effects lie outside the scope of the present formulation but may influence the reconstructed fields, and thus warrant careful consideration when assessing the fidelity and applicability of quasi-1-D shock tube models. In this section an assessment of the expected importance of these effects is provided for the two test conditions analysed in § 4.
With respect to the relative importance of transverse disturbances produced by the driver, Bowman (Reference Bowman1966) provides a criterion (equation (7.3) in that reference) for determining whether the flow in a given shock tube experiment will be dominated by unsteady transverse waves, quasi-steady boundary-layer effects or viscosity. The criterion is as follows:
\begin{equation} 0.6 \lt \frac {L}{p_1 r^2} \lt 2.3, \end{equation}
where
$L$
is the distance from the primary diaphragm in mm,
$p_1$
is the fill pressure in mm Hg and
$L$
is the distance from the primary diaphragm to the location in question in mm. For the air (§ 4.3) and Titan (§§ 4.4 and 5) cases analysed in the present work, the criterion was found to be 3.76, 4.80 and 7.60 at the observation window, respectively, placing all cases inside Bowman’s boundary-layer regime. Therefore, transverse waves from the driver and/or highly three-dimensional boundary-layer growth that might occur in a viscosity-dominated regime are not expected to be significant, and it is justified to model the cases as quasi-1-D.
Given that these experiments lie in the boundary-layer regime, the expected thickness of the boundary layer should be estimated to determine if the boundary layer intrudes appreciably into the core flow. Hooker (Reference Hooker1961) provides the following expression for evaluating the thickness of a laminar boundary layer in a shock tube flow which is compatible with the Mirels model used in the present work:
\begin{equation} \delta = \beta \left (\frac {\mu _w \ell _m}{\rho _{w,0} (U_s - u_{e,0}) }\right )^{1/2}, \end{equation}
where
$\delta$
is the boundary-layer thickness and
$\beta$
,
$\mu$
,
$\rho _{w,0}$
,
$U_s$
,
$u_{e,0}$
and
$\ell _m$
are as defined in § 2.2. Evaluating this expression for the air (§ 4.3) and Titan (§§ 4.4 and 5) cases returns expected boundary-layer thicknesses of 0.062 , 0.117 and 0.13 mm, respectively. These values correspond to less than 0.1 % of the tube diameter, and thus it is appropriate to neglect the boundary layer in these cases. Additionally, this analysis is expected to overestimate the boundary-layer thickness as it assumes fully developed flow, i.e. the test slug length is maximised.
5.4. Summary
The results presented in this section show that, at least for one test case, uncertainty in measurements of shock speed obtained via finite difference contributes appreciably to uncertainty in non-equilibrium peak radiance. However, uncertainty in reaction rate constants, particularly for the nitrogen dissociation reaction, becomes the dominant contributor away from the shock. This promises that given a combination of sufficiently precise measurement of the shock speed and use of the method described in this paper, uncertainty in test slug properties due to the shock trajectory could essentially be eliminated.
Despite the improved agreement offered by the new method against experiment when compared with previous attempts, clearly further improvements could be made. The most promising avenue would be to recast this method to solve the two- or three-dimensional viscous Navier–Stokes equations. In addition to capturing the multi-dimensional effects present in shock tubes such as shock curvature, such an approach would resolve the reacting boundary layer and remove the need to rely on the Mirels model. The Rankine–Hugoniot post-shock boundary condition, which is not realised in practice due to diffusion, could be removed as the inclusion of viscous terms would inherently capture the shock thickness (Clarke et al. (Reference Clarke, Brody, Steer, McGilvray and Di Mare2024) report the importance of this in detail). Finally, having obtained a full two-dimensional axisymmetric or three-dimensional flow field the full capability of the NEQAIR program to simulate non-local, non-Boltzmann radiative transport could be deployed. Absorption of radiation by the boundary layer has been postulated as a cause of discrepancy between simulation and experiment by several authors (e.g. Jelloian et al. Reference Jelloian, Minesi, Spearrin, Tibere-Inglesse, MacDonald and Cruden2023); however, a definitive link has proved elusive due to the lack of a numerical method that could isolate the many non-ideal effects present in typical low-pressure shock tube flows. Such an approach will ultimately enable the excitation and de-excitation pathways most responsible for discrepancies with experiment to be isolated, thereby identifying those reaction and relaxation coefficients that may be candidates for modification or optimisation against high-fidelity shock tube data.
6. Conclusion
A novel, rigorous approach to shock tube modelling has been validated against established numerical tools and applied to several test cases. It has been demonstrated that, via a reverse time integration method, it is possible to extract valuable information about wave trains originating from the driver gas by utilising available pressure data as a boundary condition. This method removes the need for non-physical constraints on pressure that are present in previous attempts, notably Satchell et al. (Reference Satchell, Glenn, Collen, Penty-Geraets, McGilvray and Di Mare2022b ), which allowed for the implementation of non-equilibrium thermochemistry. This approach allowed the effects of shock speed variation on highly non-equilibrium shock tube tests to be studied for the first time, extending the work of Collen et al. (Reference Collen, Satchell, Di Mare and McGilvray2022). It was shown that, for a Titan entry experiment, the method provided improved predictions of radiance in the far post-shock region because of the reduced CN number density. If the shock speed were modelled as constant, this effect could only be achieved by erroneously modifying reaction and relaxation rate constants. Finally, the method was used to produce uncertainty bounds on the post-shock radiance profile for a Titan test due to uncertainty in the shock trajectory and uncertainty in the underlying reaction rate constants. Uncertainty in experimental shock speed measurements was found to contribute appreciably close to the shock for a particular test case, but was greatly outweighed by uncertainty in reaction rate constants further from the shock.
In light of the results presented herein, the following recommendations for future work are made:
-
(i) Efforts to reduce uncertainty in experimental shock speed measurements should be undertaken, with the implementation of higher-fidelity techniques such as that of Dufrene et al. (Reference Dufrene, Holden and Ringuette2015) being particularly appropriate for use with the new method.
-
(ii) High-fidelity pressure transducers with short response times should be installed at the optical measurement location in shock tubes so that they can later be used to reconstruct the forward-running waves.
-
(iii) Optimisation algorithms for determination of thermochemical reaction and relaxation rates against shock tube data should take into account shock speed non-uniformity by utilising a method such as the one presented as the underlying fluid model.
Acknowledgements
The authors would like to thank M. Satchell for his time and support, S. Brody for his contributions to the OCEAN thermochemistry library, and T. Drezet for his development of the incomplete factorisation method.
Funding
The authors acknowledge the European Space Agency (ESA) for financial support of the first author and for supporting this research.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Structure of the eigenvectors
Let the subscripts
$n+1$
,
$n+2$
,
$n+3$
represent the vibrational–electronic energy equation, the momentum equation and the total energy equations, respectively. Let subscripts
$1-n$
represent the continuity equations. Assume eigenvalues of
$\lambda _{1-n,n+1} = u$
,
$\lambda _{n+2,n+3} = u\pm a$
. To find the left eigenvector
$\boldsymbol{l}_{\boldsymbol{n+1}}$
for
$\lambda _{n+1}$
we must solve
$\boldsymbol{l}_{\boldsymbol{n+1}} (A - \lambda I)= 0$
:
\begin{equation} \boldsymbol{l}_{\boldsymbol{n+1}} (A - \lambda I)=\begin{bmatrix} \ell _1\kern-3.5pt & \ldots\kern-3.5pt & \ell _n & \ell _{n+1} & \ell _{n+2} & \ell _{n+3} \end{bmatrix} \; \begin{bmatrix} 0 & \ldots & 0 & 0 & \rho _1 & 0 \\ \vdots & \ddots &\vdots & \vdots & \vdots &\vdots \\ 0 & \ldots &0 & 0 & \rho _n & 0\\[3pt] 0 &\ldots &0 & 0 & \frac {p_e}{\rho C_{v,V}} & 0\\[3pt] 0 &\ldots &0 & 0 & 0 & 1/\rho \\ 0 &\ldots &0 & 0 & \rho a^{2} & 0 \end{bmatrix} = 0. \end{equation}
The constraint resulting from the
$n+1$
column of
$A$
is
\begin{equation} l_1 \rho _1 + \ldots + l_n \rho _n + l_{n+1} \frac {p_e}{\rho C_{v,V}} + l_{n+3} \rho a^{2} = 0 .\end{equation}
The eigenvector
$\boldsymbol{l}_{\boldsymbol{n+1}}$
is then
\begin{equation} \boldsymbol{l}_{\boldsymbol{n+1}} = \begin{bmatrix} 0 & \ldots & 0 & 1 & 0 & -\frac {p_e}{\rho ^{2} a^{2} C_{v,V}} \end{bmatrix}. \end{equation}
The resulting matrix of left eigenvectors:
\begin{equation} L = \begin{bmatrix} 1 & \ldots & 0 & 0 & 0 & -\frac {\rho _1}{\rho a^{2}} \\ \vdots & \ddots &\vdots & \vdots & \vdots &\vdots \\ 0 & \ldots & 1 & 0 & 0 & -\frac {\rho _n}{\rho a^{2}}\\[3pt] 0 & \ldots & 0 & 1 & 0 & -\frac {p_e}{\rho ^{2} a^{2} C_{v,V}}\\[3pt] 0 &\ldots &0 & 0 & \rho a & 1\\ 0 &\ldots &0& 0 & -\rho a & 1 \end{bmatrix}. \end{equation}
As implemented in the program:
\begin{equation} L = \begin{bmatrix} 1 & \ldots & 0 & 0 & 0 & -\frac {\rho _1}{\rho a^{2}} \\ \vdots & \ddots &\vdots & \vdots & \vdots &\vdots \\ 0 & \ldots & 1 & 0 & 0 & -\frac {\rho _n}{\rho a^{2}}\\[3pt] 0 & \ldots & 0 & 1 & 0 & -\frac {p_e}{\rho ^{2} a^{2} C_{v,V}}\\[3pt] 0 &\ldots &0 & 0 & 1 & \frac {1}{\rho a}\\[3pt] 0 &\ldots &0& 0 & 1 & - \frac {1}{\rho a} \end{bmatrix}. \end{equation}
Appendix B. Method of determining time of arrival
The approach to determination of shock arrival time and uncertainty varies between research groups. This includes the measurement technique used (Dufrene et al. Reference Dufrene, Holden and Ringuette2015; James et al. Reference James, Cox, Komonen, Barltrop, Wikner and McIntyre2021) and the post-processing approach (Brandis Reference Brandis2009; James Reference James2019).
The method used for determination of shock arrival time and uncertainty for all T6 test cases presented in the present work is summarised in Appendix G of James (Reference James2018). Given the fundamental importance of the experimental shock speed measurement to the LASTA 2.0 simulation approach, this method is summarised here. The shock speed is simply calculated as
\begin{equation} V_s = \frac {x_2 - x_1}{t_2 - t_1} = \frac {\Delta x}{\Delta t}, \end{equation}
where
$x_1$
and
$x_2$
are the axial locations of two adjacent sensors and
$t_2$
and
$t_1$
are the shock arrival times at the respective sensors. The uncertainty in shock speed is calculated using
\begin{equation} \delta V_s = V_s \sqrt {\left (\frac {\delta \Delta x}{\Delta x} \right )^{2} + \left (\frac {\delta \Delta t}{\Delta t} \right )^{2}}, \end{equation}
where
$\delta \Delta x$
represents the uncertainty in
$\Delta x$
and
$\delta \Delta t$
represents the uncertainty in
$\Delta t$
. The uncertainty in
$\Delta x$
is comprised of:
-
(i) Uncertainty in the measurement of the sensor locations.
-
(ii) Uncertainty in the response of the pressure sensor due to its finite size.
-
(iii) Uncertainty in the curvature of the shock.
These three contributions are assigned a single uncertainty of
$\pm 2.0 \times 10^{-3}$
m for each sensor location. Given that the distance uncertainties at each sensor location are independent, the total distance uncertainty (
$\delta \Delta x$
) is given by
\begin{equation} \delta \Delta x = \sqrt {\delta x_{1}^{2} + \delta x_{2}^{2}}. \end{equation}
Similarly, one source of time uncertainty and one source of error are considered for
$\delta \Delta x$
:
-
(i) Uncertainty in the shock time of arrival due to the post-processing approach.
-
(ii) Sampling rate error from the data acquisition system.
The time of arrival is determined manually to avoid the often large uncertainty and unreliability associated with an automated process. During post-processing, the user is presented with a graphical interface showing the time history for each sensor and they are prompted to select the data points immediately before and after they believe the shock arrived. This results in a range of possible arrival times, from which the midpoint is treated as the nominal value and half the width is treated as the uncertainty. This is illustrated in figure 25( a). The tangent method of determining shock arrival is also shown.
The sampling rate error is accounted for by including the size of a single sample (
$\delta t_{sr}$
) since the shock may arrive at any point during the sample. This must be considered for both sensors. The T6 data acquisition samples at 2.5 MHz so
$\delta t_{sr} = 0.4$
μs. The final time uncertainty is then given by
\begin{equation} \delta \Delta t = \sqrt {\delta t_{1}^{2} + \delta t_{2}^{2} + 2\delta t_{sr}^{2}}. \end{equation}
Figure 25. (a) Window method. (b) Tangent method.
Appendix C. Approach to determination of the relaxation time constant for an interacting pair
Values for the characteristic temperature for each molecular species found using the harmonic vibrational frequency found in NIST (2022) are used to generate values usable in the Millikan and White equation for translational–vibrational relaxation time (Millikan & White Reference Millikan and White1963), using the form
\begin{equation} p\tau _{MW,ij}=\exp \left [A_{\textit{ij}}(T^{-1/3}-b_{\textit{ij}})-18.42 \right ], \end{equation}
where
$A_{\textit{ij}}$
and
$b_{\textit{ij}}$
are specified for each interacting pair
${\textit{ij}}$
by
\begin{align} A_{\textit{ij}} &=0.00116\mu _{\textit{ij}}^{1/2}\theta _i^{4/3}, \end{align}
\begin{align} b_{\textit{ij}} &=0.015\mu _{\textit{ij}}^{1/4}, \end{align}
where
$\mu$
is the reduced molecular mass between the two species and
$\theta _i$
is the characteristic temperature of the dissociating species. These approximate values are only used where values for specific pairs are not available in the literature. The model described by equation (56) in Gnoffo et al. (Reference Gnoffo, Gupta and Shinn1989) is utilised to modify the translational–vibrational relaxation time at high temperatures, which requires a value for the effective cross-section for vibrational relaxation
$\sigma _v$
. The form of
$\sigma _v$
(in
${\textrm m}^2$
) is taken from Park (Reference Park1993):
\begin{equation} \sigma _v=3\times 10^{-21} \left (\frac {50\,000}{T}\right )^2. \end{equation}
Where available, electron-neutral cross-sections are taken from Gnoffo et al. (Reference Gnoffo, Gupta and Shinn1989) for use in the Appleton and Bray correlation (Appleton Reference Appleton1966), used to estimate the relaxation time between translational and electronic energy modes. The remaining species are considered as trace species and thus are assumed to have a negligible impact on the energy transfer.
Appendix D. Approach to determination of
$e_{V,s}$
Energy
$e_{V,s}$
is the vibrational–electronic energy of each species. Energy
$e_{V,s}$
may be evaluated as
\begin{equation} e_{V,s} = \int _{T_{ref}}^{T_V} C_{v,V}^{s}\, {\textrm d}T^{\prime} .\end{equation}
Here,
$C_{v,V}^{s}$
(
$\textrm{J kg}^{-1}\,K^{-1}$
) is the specific heat at constant volume for species
$s$
for vibrational–electronic energy,
$T_{ref}$
is a curve-fit reference temperature, typically 298 K, and
$T_V$
is the vibrational–electronic temperature. The vibrational–electronic specific heat at constant volume for species
$s$
can be evaluated by utilising the curve fit for specific heat at constant pressure evaluated at temperature
$T_V$
and subtracting both the contribution for the translational and rotational enthalpy evaluated at
$T_V$
:
\begin{align} C_{v}^{s} (T_V) &= C_{p}^{s} (T_V) - \frac {\bar {R}}{M_s}, \end{align}
\begin{align} C_{v,V}^{s} (T_V) &= C_{v}^{s} (T_V) - C_{v,t}^{s} - C_{v,r}^{s}, \end{align}
where
$C_{v,v}^{s}$
,
$C_{v,t}^{s}$
and
$C_{v,r}^{s}$
are the specific heats at constant volume for species
$s$
for vibrational, translational and rotational enthalpy, respectively (
$\textrm{J kg}^{-1}\,K^{-1}$
),
$C_{p}^{s}$
is the specific heat at constant pressure for species
$s$
,
$\bar {R}$
is the universal gas constant (
$8314.3$
$\textrm{J mol}^{-1}\,K^{-1}$
) and
$M_s$
is the molecular weight of species
$s$
. Assuming that the translational and rotational modes are fully excited,
$C_{v,t}^{s}$
and
$C_{v,r}^{s}$
can be evaluated simply as
\begin{equation} C_{v,t}^{s} = \frac {3 \bar {R}}{2 M_s} ,\quad C_{v,r}^{s} = \frac {\bar {R}}{M_s}. \end{equation}
Term
$C_{p}^{s}$
can be evaluated from curve fits/tabulated values from the literature. In the present work we used the NASA JANAF polynomials (Gordon & McBride Reference Gordon and McBride1999) which are of the form
\begin{equation} C_{p}^{s} (T) = \frac {\bar {R}}{M_s} \sum _{k=1}^{5} A_{k}^{s} T^{k-1}. \end{equation}
The specific heat at constant volume for species
$s$
for vibrational–electronic energy is then
\begin{equation} C_{v,V}^{s} (T) = \frac {\bar {R}}{M_s} \left (\sum _{k=1}^{5} A_{k}^{s} T^{k-1} - \frac {7}{2}\right ) .\end{equation}
The vibrational–electronic energy of a species:
\begin{equation} e_{V,s} = \int _{T_{ref}}^{T_V} \frac {\bar {R}}{M_s} \left (\sum _{k=1}^{5} A_{k}^{s} T^{k-1} - \frac {7}{2}\right ) {\textrm d}T^{\prime} .\end{equation}
And the total vibrational–electronic energy per unit mass:
\begin{equation} e_V = \sum _{s=i}^{n} \left (\frac {\rho _s}{\rho } \int _{T_{ref}}^{T_V} \frac {\bar {R}}{M_s} \sum _{k=1}^{5} A_{k}^{s} T^{k-1} - \frac {7}{2}\right ) {\textrm d}T^{\prime} .\end{equation}
Open access
