Near-continuum, hypersonic oxygen flow over a double cone simulated by direct simulation Monte Carlo informed from quantum chemistry

Abstract

A large-scale, fully resolved direct simulation Monte Carlo (DSMC) computation of a non-equilibrium, reactive flow of pure oxygen over a double cone is presented. Under the simulated near-continuum conditions, the computational demands are shown to be significant because of the wide range of length scales that must be resolved. Therefore, robust grid adaption capabilities and efficient parallelization of the Stochastic PArallel Rarefied-gas Time-accurate Analyzer (SPARTA) code that is utilized in this work are essential. The thermochemical and transport collision models were selected for efficiency and simplicity. First-principles data, obtained from the highly accurate direct molecular simulation method, were used to inform the collision models’ parameters. Importantly, because SPARTA implements molecular collision models using collision-specific energies, the resulting macroscopic relaxation rates were evaluated a posteriori via zero-dimensional heat bath simulations. The comparisons of surface properties, namely heat flux and pressure, show very close agreement with previous computational fluid dynamics (CFD) results. Differences with the measurements were found to be similar to the CFD simulations. The unresolved discrepancy with the measurements could be due to inconsistent free stream conditions with the actual experimental data or missing physical phenomena altogether, for example atomic and molecular oxygen electronically excited states, three-dimensional effects, or more complex gas–surface interactions. As shown in this work, the advantages of obtaining a DSMC particle solution for these flows reside in the method's ability to be directly informed from first principles and to seamlessly describe internal energy non-equilibrium for all modes. With the advent of exascale computing and beyond, particle methods will be an increasingly important tool to verify the validity of physical assumptions in reduced-order models via fully resolved, experimental-scale simulations, down to the level of molecular-level distributions.

1. Introduction

Shock waves interactions with a boundary layer (SBLI) represent a key feature of hypersonic flows. The physical phenomena that characterize SBLI are highly complex because of various nonlinear couplings between the fluid mechanics and real gas effects, such as vibrational energy excitation and chemical reactivity. The physical description of these flows is further complicated by the fact that characteristic fluid time scales approach those of internal energy relaxation and chemical reactions, thus resulting in local thermochemical non-equilibrium. Failure to properly describe the finite rate nature of these processes leads to incorrect numerical predictions (Holloway, Hanquist & Boyd 2020).

An understanding of SBLI is also crucial from a technological point of view for the design of hypersonic vehicles. These interactions may, in fact, have a dominant role in thermal protection system heating because they produce local heat transfer peaks of extreme severity. An illustrious example of the importance of SBLI is the X-15 airplane in the 1960s (Anderson 2006, pp. 395–396).

For these reasons, the simulation of flows dominated by SBLI has been traditionally viewed as a benchmark for analytic and, more importantly, computational methods. In the case of simulations, because of the wide range of scales involved in such flows, it was observed that adequate spatial resolution was needed to obtain a reliable numerical prediction. For continuum methods (computational fluid dynamics (CFD)), this results in a requirement for a rigorous grid convergence analysis (Nompelis, Candler & Holden 2003). For particle methods (direct simulation Monte Carlo (DSMC)), a judicious selection of numerical parameters, primarily number of particles in each grid cell and mean collision separation of the order of the local mean free path, must be made (Moss & Bird 2005).

A classic set of experimental data was obtained on a $25^{\circ }/55^{\circ }$ double cone test object in the Calspan-University at Buffalo Research Center (CUBRC) Large Energy National Shock (LENS) tunnel (Holden 1999). The particular geometry was selected to provide stable flows (Holden & Wadhams 2003), although this has been partially disputed in later studies (Tumuklu, Levin & Theofilis 2018a; Tumuklu, Theofilis & Levin 2018b). The low-pressure free stream conditions resulted in separated, laminar flows due to the relatively low Reynolds numbers. Measurements for flows of pure nitrogen, pure oxygen, and air were obtained. The strong expansion in the nozzle of the facility led to a significant degree of non-equilibrium in the free stream conditions, which had to be determined with nozzle simulations (Nompelis et al. 2003). This was a significant cause of uncertainty in the measurements obtained on the test object, but heat transfer rates and pressure were measured on the double cone with a quoted uncertainty of approximately 5 % (Holden 1999).

Due to the variety of conditions, a vast number of numerical studies have focused on the simulation of the CUBRC double cone experiments. For example, a review of CFD predictions for select test cases is detailed by Knight et al. (2012). In general, both surface heat flux and pressure distributions were found to be in relatively close agreement with the experimental data. Also, the separation region was well predicted for most experimental cases. Various investigators further improved the numerical predictions by refining the simulation conditions, for example, by accounting for partial vibrational energy accommodation at the wall (Nompelis et al. 2003), non-uniformity of the free stream (Nompelis et al. 2003), better characterization of the non-equilibrium in the free stream (Candler et al. 2002; Ninni et al. 2022) and updated thermochemical models (Gao et al. 2021). However, the agreement between experiments and computations was shown to worsen as the free stream enthalpy was increased (Nompelis et al. 2010).

Despite the overall good success of CFD at simulating the CUBRC double cone experiments, several quantitative discrepancies for various runs were not resolved. This prompted uncertainty quantification analyses to verify that the stated experimental conditions were consistent with the measurements (Ray et al. 2020). Also, the numerical predictions for many runs appeared to be rather insensitive to the thermochemical models and the small variations were not enough to account for the differences with experiments (Kianvashrad & Knight 2019; Holloway, Chaudhry & Boyd 2022). Finally, the experimental run time (around 1 ms) was brought into question as insufficient to establish a steady-state flow over the double cone (Tumuklu et al. 2018a).

At variance with CFD, application of DSMC to the simulation of the CUBRC runs has proven much more challenging due to the computational costs associated with the application of DSMC to simulate near-continuum flows. In fact, the free stream conditions in all CUBRC experimental tests are characterized by densities of roughly $1\ {\rm g}\ {\rm m}^{-3}$ and translational temperatures that vary between less than 100 K and a few hundred kelvin. Therefore, the molecular mean free path in the free stream is around $O(10^{-5})$ m. Consequently, the Knudsen number based on the double cone diameters is in the range between approximately $10^{-5}$ and $10^{-4}$.

Early attempts at DSMC simulation of double cone SBLI features resulted in significant differences between the numerical predictions for the separation and the measurements (Wang et al. 2001; Gimelshein et al. 2002b). However, it was soon recognized that additional spatial resolution was needed to correctly predict the separation of the flow on the first cone. The study of Moss & Bird (2005) illustrated this by a rigorous analysis of the dependence of various fluid mechanic features of the SBLI on the DSMC numerical parameters, particularly the mean collision separation between colliding particles. Similar conclusions were shown by Roy et al. (2003), who compared DSMC and CFD, but used reduced free stream densities compared with the typical experimental values. Notably, the work of Moss & Bird (2005) was based on the simulation of Run 7, which had the lowest free stream pressure of all the CUBRC experimental runs and consisted of pure nitrogen with negligible chemical activity in the shock layer. Run 7 was also selected by Plimpton et al. (2019) to demonstrate the Stochastic PArallel Rarefied-gas Time-accurate Analyzer (SPARTA) DSMC code capabilities to fully resolve the complex fluid dynamic features of hypersonic double cone flows.

The considerable advancements in computational capability in the last 20 years are enabling the application of DSMC to the near-continuum regime. In the case of double cone flows, very large-scale computations were recently conducted by Tumuklu et al. (2018b) to investigate the effect of Reynolds number on the SBLI structure and its stability. In the work of Tumuklu et al. (2018b), up to $10\times 10^9$ particles were used and the simulation was time integrated for several hundred thousand time steps on grids containing between $300\times 10^6$ and $600\times 10^6$ collision cells. A minimum of 5–6 particles per collision cell was maintained. Since the focus of the study was to investigate fluid dynamic instabilities, the selected free stream conditions were those of Run 7, i.e. characterized by strong vibrational non-equilibrium in the shock layer, but no significant chemical reactivity. Only the free stream static pressure was increased to obtain larger Reynolds numbers. On the other hand, Gimelshein & Wysong (2019) simulated high-enthalpy chemically reacting air flows over the $25^{\circ }/55^{\circ }$, double cone. The computations, similar to the work of Tumuklu et al. (2018b), were conducted with the SMILE DSMC code (Ivanov, Markelov & Gimelshein 1998) equipped with sophisticated DSMC models. Approximately $500\times 10^6$ particles were used, but no further information is provided with respect to grid size and number of particles per collision cells. However, Gimelshein & Wysong (2019) state that sensitivity studies were conducted to obtain a solution independent of the DSMC parameters. A comparison between the low and high fidelity DSMC models revealed a relative insensitivity of the flows to the thermochemical models, and thus, differences with the experimental data could not be resolved.

Quite recently, the direct molecular simulation (DMS) method (Schwartzentruber, Grover & Valentini 2018) was applied to two-dimensional, reactive flows in thermochemical non-equilibrium. The DMS method can be viewed as an ab initio variant of DSMC because it replaces simplified collision models with actual molecular dynamics trajectories integrated on an ab initio potential energy surface (PES). Due to the cost of trajectory integration, the DMS method is significantly more computationally expensive than standard DSMC. Initially, it was applied to canonical flows around a cylinder with a $O(10^{-2})$ free stream Knudsen number (Grover & Valentini 2021; Valentini et al. 2021; Grover et al. 2022b). Later, a Mach 21 flow of dissociating nitrogen over a blunt wedge was simulated with DMS. In this case, the free stream Knudsen number based on the wedge nose radius was $O(10^{-3})$ (Grover et al. 2022a).

It is clear from the discussion so far that the application of DSMC to near-continuum flows still represents a significant computational challenge due to the complexity resulting from the wide range of length scales that must be resolved. In terms of fidelity, DSMC is generally positioned between the highly accurate molecular dynamics (MD) and DMS methods, and engineering CFD solvers. Because DSMC relies on fewer assumptions than CFD in terms of thermochemical models and vicinity to local thermodynamic equilibrium, and it naturally incorporates thermal fluctuations, it remains extremely valuable to compare the DSMC statistical solutions with those obtained from CFD at conditions where both methods are fully applicable. Furthermore, the DSMC collision models tend to be more easily informed from ab initio DMS or quasiclassical trajectory calculation data at the level of collision cross-sections (Grover & Valentini 2021).

Therefore, the objective of this work is to simulate a flow in thermochemical non-equilibrium over a $25^{\circ }/55^{\circ }$ double cone with significant chemical reactivity in the shock layer. The free stream static pressure leads to a near-continuum flow field, thus enabling well-posed comparisons with CFD. Furthermore, the free stream enthalpy is sufficient to cause significant chemical reactivity in the shock layer. The DSMC models were selected for computational efficiency and, therefore, can be viewed as lower fidelity compared with more advanced formulations (Gimelshein & Wysong 2019). However, where possible, the parameters were optimized to reproduce previous computational results for transport phenomena, internal energy relaxation, and chemical dissociation obtained from highly accurate, ab initio PES.

The article is organized as follows. Section 2 contains a description of the computational method and the details of the thermochemical and transport models utilized in the computations. A brief summary of the code employed is also provided. The flow conditions and geometry are listed in § 3. The results are summarized in § 4 and the conclusions in § 5.

2. Simulation method

The DSMC method was originally proposed by Bird (1994) for the simulation of rarefied and transitional flows. Although initially restricted to rarefied conditions because of its computational cost, DSMC has recently found applications in areas of fluid dynamics near the continuum regime (Gallis et al. 2017), and particularly for problems dominated by thermochemical non-equilibrium. In this work, we use the open-source SPARTA DSMC code developed by Plimpton et al. (2019), as it was designed with the goal of executing efficiently on the largest high-performance computing platforms. The code version utilized for all simulations presented here is 9 September 2019.

The DSMC method is an inherently time accurate method, as the simulation time step must be of the order of the smallest mean collision time in the flow. For this reason, cells in the flow field must have a size of the order of the local mean free path. Unlike deterministic MD, only a subset of real particles in each control volume is simulated. Every simulated particle can be viewed as representing a large number of gas molecules and atoms. This is controlled by assigning a so-called particle weight, which can be made cell-specific, as in the case of axisymmetric flows. Similarly, every collision can be viewed as representing a large number of near-identical molecular interactions in each collision cell. The no time counter (NTC) algorithm (Bird 1994) is used to stochastically select collision partners. The collisions occurring in each cell at every time step, in turn, drive the time evolution of velocity and energy distribution functions. In this light, DSMC can be viewed as a numerical solution of the Boltzmann equation in that the collision operator is modelled stochastically via a large number of molecular interactions, both elastic and inelastic. Although this was only rigorously proved for monoatomic (structureless) gases by Wagner (1992), there is overwhelming anecdotal evidence on polyatomic, chemically reactive systems (Plimpton et al. 2019).

In DSMC, a mapping from initial (reactant) to final (product) molecular states due to a collision must be performed. Analytic functions of a set of precollision quantities, e.g. total collision energy and impact parameter, are generally used and are commonly referred to as collision models. Numerous collision models of various sophistication (Bird 1994; Boyd & Schwartzentruber 2017) have been proposed since the inception of the method, particularly for internal energy relaxation (both rotational and vibrational modes) and chemical reactivity.

2.1. Viscosity cross-section model

The variable hard-sphere (VHS) cross-section model (Bird 1994) is used for particle collisions. The model parameters, namely reference diameter and viscosity index, were selected based on the recent ab initio DMS data of Valentini et al. (2023) and are listed in table 1. The resulting shear viscosity predictions for molecular oxygen are shown in figure 1, where the VHS and DMS data align very well. Similar agreement is found for atomic oxygen.

Table 1. The DMS-informed VHS collision-specific reference temperature and diameter, and viscosity index for oxygen.

Figure 1. The DMS-based viscosity curves used to parametrize the VHS model.

2.2. Rotational and vibrational relaxation models

The DSMC models for internal energy relaxation are often formulated based on a collision number $Z_i$ (Bird 1994), with $i=\{r,v\}$ for rotation and vibration, respectively. Due to the probabilistic nature of DSMC, particle collisions are either considered elastic (i.e. no translational collision energy transfer with internal energy modes) or inelastic (e.g. translation–vibration or translation–rotation energy transfers). Therefore, the collision number can be viewed as the average number of collisions a molecule must undergo before energy is exchanged between internal and external modes. Under this interpretation, the reciprocal of the collision number $1/Z_i$ can be approximately viewed as the probability for such energy exchanges.

At the continuum level, the transfer between the average internal energy $\overline {E_i}(t)$ and the average equilibrium energy $\overline {E}^*_i(T)$ is typically modelled with a time relaxation differential equation,

(2.1)\begin{equation} \frac{{\rm d} E_i}{{\rm d}t} = \frac{\bar{E}^*_i(T)-\overline{E_i}(t)}{\tau_i}. \end{equation}

For $i=r$, (2.1) is typically referred to as Jeans’ equation, whereas the Landau–Teller equation is for $i=v$.

A connection between the microscopic information of $Z_i$ and the macroscopic characteristic relaxation time that appears in (2.1) is recovered by using the relation

(2.2)\begin{equation} \tau_i = C_f Z_i \tau_c, \end{equation}

where $\tau _c$ is the mean collision time of the gas. In general, the correction factor $C_f\ne 1$, and it depends on both the functional form of $Z_i$ and the implementation of the relaxation process (e.g. permitting or prohibiting double relaxation), as discussed by other authors (Lumpkin, Haas & Boyd 1991; Haas et al. 1994; Gimelshein, Gimelshein & Levin 2002a; Zhang & Schwartzentruber 2013). Equivalently, the probability of relaxing the energy mode must not be simply taken as $1/Z_i$, if one wants to precisely impose macroscopic relaxation rates in DSMC, albeit in a statistical manner.

Relaxation number models generally account for the influence of the gas state on the relaxation process via a functional dependence on either the so-called cell temperature (i.e. a suitable cell average of particle energies, both translational and internal) or the pair-specific collision energy. In both cases, the correction factor must be determined if a known macroscopic relaxation time $\tau _i$ is to be imposed. The task of obtaining $C_f$ may be a significant challenge except for the simplest cases (e.g. constant collision numbers), particularly for collision-specific relaxation numbers.

The SPARTA code simply implements $1/Z_i$ as the probability of relaxing mode $i$ using a pair-selection scheme that allows for double relaxation and the Borgnakke–Larsen method (Borgnakke & Larsen 1975) to distribute the total collision energy between the collision partners. Therefore, imposing a known relaxation rate, for example the Millikan–White correlation for translation–vibration energy transfer (Millikan & White 1963), is not possible a priori without the determination of $C_f$. For this work, however, the estimation of effective continuum relaxation times resulting from the collision-specific models utilized in SPARTA is done a posteriori using isothermal heat bath simulations. Where possible, model parameters were further calibrated to obtain relaxation rates in accordance with ab initio DMS data.

For rotational relaxation, we used a constant rotational relaxation number $Z_r = 4$ (Bird 1994) as it was shown to be relatively independent of temperature in previous MD studies (Valentini, Zhang & Schwartzentruber 2012; Kosyanchuk & Yakunchikov 2021). For vibration–translation energy exchange, the vibrational relaxation number $Z_v$ was modelled according to Bird (1994)

(2.3)\begin{equation} Z_v(T) = \frac{B_1}{T^{\omega}} \exp\left( \frac{B_2}{T^{{1}/{3}}}\right), \end{equation}

where $T$ is generally defined as the cell-averaged translational temperature, $\omega$ is the VHS viscosity index, and $B_1$ and $B_2$ are adjustable parameters. Because SPARTA is not based on cell-averaged temperatures, the value for $T$ appearing in (2.3) is computed as

(2.4)\begin{equation} T = \frac{\frac{1}{2} m_r c^2_r + E_v}{k_B \left( \frac{5}{2} - \omega + \overline{\zeta_v}\right)}, \end{equation}

where $E_v$ is the vibrational energy in the inelastic collision. Then, $E_v$ is summed to the relative kinetic energy of the two colliding particles $1/2 m_r c^2_r$ (where $m_r$ is the reduced mass and $c_r$ the relative speed). The energy term is rescaled with Boltzmann's constant $k_B$ and the degrees of freedom $5/2 - \omega + \overline {\zeta _v}$ (Bird 1994) involved in the collision. As discussed by Eckert & Gallis (2022), (2.4) coupled with the Borgnakke–Larsen model is shown to satisfy detailed balance.

Here, we calibrated $B_1$ and $B_2$ by estimating the vibrational relaxation rate with a set of zero-dimensional, isothermal relaxation calculations that are shown in figure 2(a). The rotational and vibrational temperatures of the system $T_{r,v} = \overline {E_{r,v}}/k_B$ were initialized to 250 K, while the translational temperature was controlled by resampling the centre-of-mass velocities from a Maxwell–Boltzmann distribution. The e-folding method was then used to compute $\tau _v$ using (2.1). To estimate $\tau _v$ for O + O$_2$ collisions, the system contained 1 % of molecular oxygen by mole and, therefore, the relaxation of O$_2$ was only due to collisions with O atoms. A similar approach was used by Valentini et al. (2015) to estimate N + N$_2$ relaxation rates using the DMS method. The resulting $p \tau _v$ data are shown in figure 2(b). For O$_2$ + O$_2$ collisions, the vibrational relaxation model used here agrees well with the recent measurements of Streicher, Krish & Hanson (2020), the Millikan–White correlation (Millikan & White 1963), and the computational results of Grover, Torres & Schwartzentruber (2019). For O + O$_2$, the model parameters in (2.3) were chosen to reproduce the recent ab initio computations of Grover et al. (2019). Table 2 contains $B_1$ and $B_2$ used in the simulations presented in § 4.

Figure 2. (a) Isothermal vibrational energy relaxation for O$_2$ + O$_2$ collisions at various equilibrium temperatures $T_{eq}$. (b) Vibrational relaxation times for O$_2$ + O$_2$ and O + O$_2$ collisions and comparisons with DMS and experimental data.

Table 2. The DMS-informed vibrational relaxation number parameters.

2.3. Dissociation model

The total collision energy (TCE) model (Bird 1994) is used in this work. The TCE reaction cross-section, specialized to the VHS model, is given by Bird (1994)

(2.5)$$\begin{gather} \sigma_{diss} = \sigma_{VHS} C_1 (E_c - E_a)^{C_2} \left(1 - \frac{E_a}{E_c}\right)^{\bar{\zeta}+3/2-\omega}, \quad\text{if}\ E_c>E_a, \end{gather}$$

(2.6)$$\begin{gather}\sigma_{diss} = 0, \quad \text{if}\ E_c \le E_a, \end{gather}$$

where $\sigma _{VHS}$ represents the VHS cross-section (a power-law of the relative speed), $E_c$ is the collision energy contained in the $\bar {\zeta }$ (average) internal degrees of freedom, and $\omega$ is the VHS viscosity index. Application of (2.5) then results in a rate coefficient of the form

(2.7)\begin{equation} {k}_{diss} = A T^{\eta} \exp{\frac{-E_a}{k_B T}}. \end{equation}

As shown by Bird (1994), $C_2$ is then simply given by

(2.8)\begin{equation} C_2 = \eta - 1 + \omega. \end{equation}

A more complex relation is also found between $C_1$ and $\text {k}_{diss}$ (Bird 1994, p. 127). Furthermore, in the SPARTA implementation, $E_c$ only contains the energy contribution from rotation, as the contribution from internal modes is arbitrary in the TCE model (Bird 1994, p. 128).

We decided to optimize $C_1$ and $C_2$ (table 3) to obtain dissociation rates that match previously published DMS data (Grover et al. 2019) obtained from first-principles PESs. Once again, this was verified with zero-dimensional, isothermal chemical reactors a posteriori, such as those shown in figure 3(a) for O$_2$ + O$_2$ dissociation. The reaction rates obtained with these DMS-informed TCE parameters are shown in figure 3(b). No recombination is modelled in the flow field computations.

Table 3. The DMS-informed TCE model parameters.

Figure 3. (a) Isothermal dissociation in a zero-dimensional chemical reactor at various equilibrium temperatures $T_{eq}$. (b) The TCE dissociation rates with parameters from table 3 and comparison with DMS data.

3. Simulation details: geometry and flow conditions

Axisymmetric simulations were conducted on a $25^{\circ }/55^{\circ }$ double cone geometry. The same configuration has been used in several previous studies (Nompelis et al. 2003, 2010; Moss & Bird 2005; Holloway et al. 2020). The conditions in the free stream were selected based on Nompelis et al. (2010) and correspond to those in Run 88, as detailed in table 4. We point out that the nominal conditions used in the simulations of Nompelis et al. (2010) exhibited a slight degree of vibrational non-equilibrium in the free stream ($T_{\infty } = 570$ K and $T_{v,\infty } = 698$ K). Here, however, we assumed $T_{\infty } = T_{v,\infty } = 606$ K, i.e. a weighted average of the rotranslational temperature and the vibrational temperature assuming equipartition of energy and classical energy modes. Full thermal and momentum accommodation were modelled by imposing diffuse reflection at the double cone surface, whose temperature was held constant at 300 K. No wall-driven chemistry was modelled either, i.e. the wall is assumed to be chemically inert. Finally, a cell-independent particle weight is not optimal for axisymmetric simulations, as it generally leads to very small (vanishing) particle counts for cells near the axis of rotational symmetry and very large particle counts for cells away from it. To circumvent this difficulty, a per-cell weighting scheme was adopted, where each cell in the domain is assigned a weight proportional to the distance of its midpoint from the axis of symmetry. Particles are then cloned or removed at each time step to maintain the correct local mass density.

Table 4. Free stream conditions for Run 88.

Figure 4 shows a magnified view of the simulated flow around the triple point with labels for the main flow features. The magnitude of the density gradient, plotted on a logarithmic scale, reveals the complex structure of the various intersecting shocks, including shocklets that form downstream of the triple point. The separation and reattachment points define the extent of the recirculation zone whose streamlines are plotted in figure 4.

Figure 4. Magnified view of the flow domain around the triple point with labels for the main flow features. Magnitude of the density gradient is shown on a logarithmic scale with streamlines of the flow inside the recirculation zone.

The entire simulated flow field is illustrated in figure 5, where Mach number contours are shown. A characteristic vibrational temperature of 2256 K (McQuarrie & Simon 1997, p. 739) was used in the harmonic oscillator approximation. The dotted line indicates the sonic condition in the flow. Subsonic flow is seen near the wall, where the gas velocity approaches zero, and in the separation zone. Furthermore, the subsonic flow after the quasinormal portion of the detached shock that is created near the triple point is rapidly accelerated to supersonic conditions, thus creating a supersonic jet very near the wall of the second cone. A three-dimensional view of the distribution of pressure is shown in figure 6(a). The pressure rise due to the laminar separation is seen near the intersection between the two cones, whereas the maximum pressure is reached downstream of the intersection. The heat flux, shown in figure 6(b), is highest at the nose of the first cone. The laminar separation causes a steep drop and the highest heat flux on the surface of the double cone is reached where the pressure peaks.

Figure 5. Mach number contours for the pure oxygen flow around a $25^{\circ }/55^{\circ }$ double cone geometry.

Figure 6. Three-dimensional view of the surface pressure and heat flux on the double cone.

3.1. Grid-independence study

The need for a rigorous grid-independence analysis has been recognized in various previous studies. For example, Nompelis et al. (2003) point out how insufficient resolution may lead to fortuitous good agreement with the experiments. On the other hand, insufficient resolution in early DSMC computations (Wang et al. 2001) led to significant discrepancies with the measurements of pressure and heat flux distributions on the surface of the double cone. In particular, the size of the recirculation zone was incorrectly predicted by the DSMC computations, particularly for the higher density cases (e.g. Run 35 or Run 42).

For CFD, a grid refinement study for these axisymmetric flows is relatively straightforward. For example, as shown by Nompelis et al. (2003) or Ninni et al. (2022) among others, the axial and normal grid spacings can be progressively refined until grid independence is observed. Such an approach is not computationally convenient for DSMC, which requires grids that resolve the local mean free path. For this reason, due to the steep density variations throughout these types of flows, a simple uniform refinement of the entire grid system would result in an unnecessarily fine grid throughout the domain. Consequently, the total number of particles would have to be extremely large in order to satisfy the usual requirement to have a minimum of ${\sim }20$ particles in each collision cell.

For these reasons, the SPARTA code implements adaptive mesh refinement capabilities (Plimpton et al. 2019). Cartesian cells are subdivided according to threshold criteria that are user-set. In this work, cells were progressively halved in each spatial direction $x$ and $y$ until the cell-based Knudsen number $Kn_c$ was greater than one nearly throughout the domain. Based on the definition of $Kn_c$,

(3.1)\begin{equation} Kn_c = \frac{\lambda_{VHS}}{\delta_c}, \end{equation}

where $\lambda _{VHS}$ is the mean free path based on the VHS model and $\delta _c$ the cell size, the condition $Kn_c > 1$ implies that $\delta _c<\lambda _{VHS}$.

An illustration of the complexity of the grid system is shown in figure 7, where the portion of the domain surrounding the triple point is magnified. The solid lines indicate the edges between levels of refinement. In general, the smallest cells are at the wall, where the mean free path is of the order of $10^{-6}$ m. However, the complex density gradients trigger grid adaptations around the triple point as well. As expected, the grid is also progressively refined around the shocks.

Figure 7. Magnified view of the flow domain around the triple point. Solid lines separate grid cells at different refinement levels.

In order to minimize the overall computational cost of the simulations, the maximum number of refinement levels (i.e. how many times a cell could be refined) was held fixed but progressively increased. In this manner, we could identify the two final grid systems that produced near-identical solutions, without unnecessarily triggering further refinement that would result in an increased total number of particles. For each grid system, the flow was advanced in time until steady-state, which was identified by both total number of particles and total number of collisions, and typically required ${\sim }50\,000$ time steps. A time step of 0.5 ns was used in the simulations, approximately half the smallest mean collision time found in the domain.

Using this procedure, we identified a coarse and a fine grid system. The coarse system was obtained by allowing seven levels of refinement from the initial uniform grid for a total of approximately $22\times 10^6$ collision cells. Around $4\times 10^9$ particles were used to guarantee a minimum number of particles/cell above ${\sim }20$, even though this is a conservative figure since solution independence was observed with $400\times 10^6$ simulated particles as well. The fine grid system was obtained by allowing two further refinement levels (nine total levels of refinement), resulting in approximately $30\times 10^6$ flow cells. The total number of particles was also increased to approximately $9\times 10^9$.

As shown in figure 8, both the pressure and the heat flux on the surface of the double cone obtained with the coarse and fine grids are nearly identical. We also plotted various isolines for temperatures (translational $T_t$, rotational $T_r$ and vibrational $T_v$), and pressure ($p$) in figure 9 that demonstrate that the two solutions overlap. We emphasize that the flow domain surrounding the triple point is the most critical flow feature, as it is affected by the strong nonlinearities associated with separation and detached shock locations. Therefore, the flow contour comparisons are quite stringent to demonstrate grid independence.

Figure 8. Comparison of (a) wall pressure and (b) heat flux obtained on the coarse and fine grid systems.

Figure 9. Comparison between coarse and fine grid isolines for (a) translational $T_t$, (b) rotational $T_r$, (c) vibrational $T_v$ temperature and (d) pressure in the flow domain surrounding the triple point.

As an additional heuristic, we looked at a statistical metric related to cell sizing for the coarse and fine grids. Specifically, we computed cumulative distribution functions (CDFs) of $Kn_c$ for cells within a prescribed distance ($d$) from the entire wall of the double cone. Similar to the region around the triple point discussed earlier, the near wall domain is also critical because of the relatively low temperature (300 K) that causes very strong flow gradients, particularly for mass density. It can be seen from figure 10(a) that only 20 % of cells within 0.1 mm of the wall do not satisfy the $\delta _c<\lambda _{VHS}$ condition in the fine grid system, whereas approximately 80 % of cells are larger than the mean free path in the coarse grid system. Farther from the wall, more and more cells satisfy the criterion, as shown in figure 10(b,c). While this analysis could be made more granular, for example, by conditioning the CDFs on axial location, the comparisons shown in figures 8 and 9 quite convincingly demonstrate grid independence of the solution. However, it is possible that a different wall boundary condition characterized by partial slip or by incomplete accommodation could require a more strict enforcement of the collision cell size requirement. Finally, we verified that the minimum number of particles in any cell of the domain was greater than 20. All results shown in § 4 were obtained on the fine grid.

Figure 10. The CDF of cell-based Knudsen number for cells within various wall distances $d$.

4. Results

Figure 11 shows the pressure field surrounding the double cone. The static pressure rapidly rises at $x\simeq 0.075$ m on the first conical surface, near its intersection with the second cone, due to the laminar separation. The highest pressures are reached on the surface of the second cone in the vicinity of the triple point, i.e. the intersection of the separation and detached shocks, where the transmitted shock impinges on the wall ($x\simeq 0.1$ m).

Figure 11. Pressure field around the double cone geometry. Flow field data are extracted along the dashed lines normal to the second cone and shown in figures 14–16.

As shown in figure 12, strong vibrational non-equilibrium, indicated by $T_t - T_v>0$, is seen in the shock layer around the first cone and is the result of the much slower relaxation of vibrational degrees of freedom compared with the translational modes. On the other hand, most of the shock layer on the second cone is in near-equilibrium ($T_t - T_v\simeq 0$), with the exception of the shock interfaces, the area near the triple point, and the strong expansion at the outflow. The rapid re-equilibration could be due to O + O$_2$ collisions, which are more effective at relaxing the vibrational modes, as shown in figure 2(b). This mechanism has been previously attributed to exchange processes, which more effectively redistribute the collision energy into internal energy (Valentini et al. 2016; Grover et al. 2019), thus reducing the characteristic vibrational relaxation time. We note that exchange reactions are not explicitly modelled in this work, but their effect on the relaxation process of vibrational energy is accounted for by a significantly reduced and nearly temperature-independent $\tau _v$ (figure 2b).

Figure 12. The DSMC prediction for the difference between translational $T_t$ and vibrational $T_v$ temperatures in the flow around the double cone. The dashed line identifies the sonic condition.

The free stream enthalpy (approximately $9\ {\rm MJ}\ {\rm kg}^{-1}$) is sufficient to cause significant molecular oxygen dissociation, as shown in figure 13. The atomic oxygen mass fraction $\chi _{\text {O}}= \rho _{\text {O}}/\rho$ (where $\rho _{\text {O}}$ and $\rho$ are the partial atomic oxygen density and the total density, respectively) reaches approximately 20 % and it is produced almost exclusively in the subsonic region emanating from the triple point in the shock layer surrounding the second cone. The weaker attached shock produced on the first cone is not strong enough to excite the internal modes of the gas sufficiently to cause any significant chemical dissociation.

Figure 13. Mass fraction of atomic oxygen around the double cone geometry predicted by DSMC. The dashed line identifies the sonic condition.

To better illustrate the degree of thermal non-equilibrium in the flow, several data lines normal to the second cone surface were extracted at various axial locations $x$, both upstream ($x=0.0932$ m and $x=0.0990$ m) and downstream ($x=0.1013$ m, $x=0.1047$ m, $x=0.1104$ m and $x=0.1162$ m) of the triple point, as indicated in figure 11. As figure 14 shows, the rotational and translational modes are in equilibrium throughout the shock layer, with the exception of the shock (shocklet) interfaces. On the other hand, the vibrational energy mode is significantly out of equilibrium across a large portion of the domain between the shock interface and the wall, as expected. Translational temperatures peak around 9000 K across the bow shock. At these energies, it is expected that transitions to electronically excited states for O and O$_2$ play a significant role. In this work, however, these transitions are not modelled.

Figure 14. Extraction lines normal to the surface of the $55^{\circ }$ cone at various axial locations $x$ for translational $T_t$, rotational $T_r$ and vibrational $T_v$ temperature.

The complex structure of the shock layer surrounding the triple point is further illustrated in figure 15, which shows pressure and Mach number of the gas along the same extraction lines identified in figure 11. The pressure at the wall rapidly rises downstream of the triple point ($x=0.0990$ m) because the transmitted shock impinges on the cone wall. The Mach number is shown to drop across the detached shock and rapidly rises downstream of the triple point ($x=0.1013$ m, $x=0.1047$ m, $x=0.1104$ m and $x=0.1162$ m) very near the wall ($d<0.003$ m). This is the supersonic jet illustrated in figure 5. Due to the diffuse wall boundary condition with full momentum and thermal accommodation, the gas then rapidly slows down at the cone surface.

Figure 15. Extraction lines normal to the surface of the $55^{\circ }$ cone at various axial locations $x$ for static pressure and Mach number.

The chemical reactivity in the shock layer is shown in figure 16, where the mass fraction of atomic oxygen is plotted along the selected extraction lines. Molecular oxygen dissociation occurs near the triple point ($x>0.0990$ m) and atomic oxygen is advected downstream. Similar to what we observed in the DMS solution for a chemically reactive flow around a blunt wedge (Grover et al. 2022a), mass fraction gradients are present at the wall, despite imposing no wall-driven chemistry. This phenomenon, known as the Soret effect (Eastman 1928) can be attributed to thermal gradients that cause a mass flux in gas mixtures. Although this closure term is usually neglected in the Navier–Stokes equations, it could instead be relevant for chemically reactive, high-speed flows.

Figure 16. Extraction lines normal to the surface of the $55^{\circ }$ cone at various axial locations $x$ for atomic oxygen mass fraction.

No previous DSMC predictions for this flow exist in the literature. Therefore, we can only compare our results with CFD data. Figure 17(a,b) shows a comparison with the CFD computations of Nompelis et al. (2010) and the more recent results of Holloway et al. (2020). Remarkable agreement with the results of Nompelis et al. (2010) is shown for both pressure and heat flux at the wall over the entire double cone geometry. The maximum pressure predicted by both methods differs by approximately 7 %. A similar level of agreement is found with the CFD heat flux predictions of Nompelis et al. (2010), with a maximum discrepancy of approximately 6 %–7 % with the DSMC results. However, as the magnified views of the regions around the peak pressure in figure 17(c) and the peak heat flux in figure 17(d) better highlight, the CFD data do not exhibit the pronounced trough shown in our DSMC results.

Figure 17. Comparison of DSMC and CFD predictions with measurements of wall pressure and heat flux.

The results of Holloway et al. (2020) for wall pressure are in very good accordance with the DSMC data, too. However, the first heat flux peak is approximately 14 % lower. Furthermore, the heat flux along the wall of the second cone is approximately 10 % lower compared with the DSMC result. Neither CFD data sets predict the significant trough shown in the DSMC heat flux around the transmitted shock impingement area.

In general, the level of agreement of the DSMC solution with the experimental data is similar to that obtained with CFD (figure 17). All computations underestimate the size of the laminar separation zone by approximately 9 mm. The peak heat flux measured in the CUBRC experiment is approximately 15 % higher than the DSMC result, which is clearly outside the ${\pm }5\,\%$ quoted uncertainty in the experiments (Holden 1999). Similarly, the heat flux along the surface of the $55^{\circ }$ cone is underestimated by the numerical predictions.

Figure 18 shows the density-based gradient-length local ($Kn_{GLL}$) Knudsen number investigated by Wang & Boyd (2003) to analyse the continuum breakdown. This parameter is defined as

(4.1)\begin{equation} Kn_{GLL} = |\boldsymbol{\nabla} \rho| \frac{\lambda}{\rho} \end{equation}

and it was found to be a good indicator of where the continuum assumption becomes inadequate to correctly describe high-speed flows. In particular, Wang & Boyd (2003) found that $Kn_{GLL} > 0.05$ reliably indicated such a failure for flows around similar geometries as that investigated here. From figure 18, it is clear that $Kn_{GLL}$ remains well below 0.05 along all the selected extraction lines, with the exception of the strong detached shock, which is identifiable by the sharp temperature jumps. However, because of the large gradients caused by the cold wall, $Kn_{GLL}$ is seen to slightly exceed 0.05 at the surface.

Figure 18. Continuum breakdown parameter along several extraction lines normal to the surface of $55^{\circ }$ cone at various axial locations $x$.

In figure 19(a,b), we plot the temperature and velocity jumps at the double cone wall. These were approximated by a difference between the wall imposed values and the cell-centred quantities for the cells adjacent to (i.e. intersected by) the wall. As expected, the largest temperature slip occurs on the surface of the first cone, particularly at the sharp nose. In the laminar separation zone ($0.075< x<0.09$ m), nearly no temperature slip is observed, due to the slow-moving, high-density gas. A more modest jump of ${\sim }50$ K is seen on the surface of the second cone. We note that the vibrational temperature jump is considerably less pronounced compared with rotational and translational temperature slip on the first cone. This is because the free stream vibrational temperature is set to 606 K and the attached shock on the first cone is not strong enough to excite the vibrational mode of the gas significantly, unlike the rotational and translational modes, thus producing a more gradual gradient at the wall. The slip velocity exhibits a similar behaviour, as shown in figure 19(b).

Figure 19. Velocity and temperature jumps along the surface of the double cone.

The trends found for this test case are in qualitative agreement with those shown by Moss & Bird (2005) and Tumuklu et al. (2018b) for Run 7. However, at variance with the conditions for Run 88, Run 7 is characterized by a much higher free stream vibrational temperature ($T_{v,\infty } = 1986$ K) that leads to a much more significant temperature jump for $T_v$ at the wall.

4.1. Velocity and energy distributions

In addition to macroscopic quantities, molecular-level distributions for velocities and internal energy can also be sampled from DSMC simulations at any location in the flow field. For example, Tumuklu et al. (2018b) report several velocity distribution functions across the bow shock for Run 7. The characteristic bimodal shape is observed at several locations within the shock front, as expected. Here, we extract distribution functions for molecular (O$_2$) velocities and internal energies at two locations in the flow, namely on the surface of the first cone at approximately 3.3 mm from the nose, magnified in figure 20(a), and near the triple point, magnified in figure 21(a). For the statistical samples, approximately 200 000 particles were collected at the first locations and approximately $1.4\times 10^6$ particles at the second location. In both cases, statistics were gathered over the course of 1000 time steps at intervals of 10 time steps.

Figure 20. (a) Extraction location for (b,c) velocity, (d) rotational energy and (e) vibrational energy distributions near the first cone's surface, approximately 3 mm from the tip.

Figure 21. (a) Extraction location for (b,c) velocity, (d) rotational energy and (e) vibrational energy distributions near the triple point.

As shown in figure 20(b,c), the $x$- and $y$-component centre-of-mass velocity distributions exhibit a significant departure from the equilibrium distribution obtained at the corresponding average translational energy. Furthermore, their skewness is indicative of the gas slip at the surface of the cone. Both rotational ($e_r$) and vibrational ($e_v$) energy distributions are shown to have non-equilibrium features, but more so $f(e_r)$, which has a slight bimodal shape. The location of the extraction, very near the cone tip, is significantly affected by rarefaction, which is revealed by the velocity distributions at the microscopic scale and by temperature jump and velocity slip, as illustrated in figure 19. The weak, thin shock still produces significant rotational excitation as demonstrated by the overpopulated tail of the distribution in figure 20(d) but very modest vibrational excitation (near-Boltzmann distribution), as shown in figure 20(e).

A much different microscopic picture is found in proximity of the triple point. As shown in figure 21(b,c), centre-of-mass velocities are shown to be in local thermodynamic equilibrium. Rotational energy also follows an equilibrium distribution, consistent with the profiles shown in figure 14. The vibrational energy distribution, however, is strongly non-Boltzmann, and it exhibits a significant over-population of the tail, indicative of the high-energy collisions taking place in the hot shock layer produced on the second cone.

5. Conclusions

In this work, a large-scale, fully resolved DSMC computation of a chemically reactive flow in thermochemical non-equilibrium in the near-continuum regime is investigated. As there is a wide range of length scales that must be resolved, the computational demands are found to be quite significant. The full resolution simulation required nearly 4000 CPUs and a total cost of around $5\times 10^6$ CPU-hrs. However, with robust grid adaption and efficient parallelization in SPARTA, we show that DSMC results for complex flows can be obtained in regimes where direct comparisons with continuum descriptions are possible.

Although selected for computational efficiency and simplicity, the thermochemical and transport models used in this work were parametrized based on previous ab initio data obtained from the highly accurate DMS method. Because SPARTA implements molecular collision models using collision-specific energies, the macroscopic internal energy relaxation rates and chemical dissociation rates were evaluated a posteriori via zero-dimensional heat bath simulations. Grid independence was demonstrated via convergence of surface properties, primitive flow variables comparisons and a heuristic based on cell-based Knudsen number.

The analysis of surface properties, namely heat flux and pressure, shows near-exact agreement with previous CFD results of Nompelis et al. (2010). Similar, but worse agreement with the CFD data of Holloway et al. (2020) is also demonstrated. Differences with the experimental data are in line with the CFD results. Specifically, both our DSMC and previous CFD predictions in the literature significantly underestimate the size of the laminar separation zone and the peak heat flux. In light of these comparisons, it appears that this flow is rather insensitive to the thermochemical and transport models, at least at the level of surface quantities. A clear demonstration of this was also offered by Gimelshein & Wysong (2019) for other double cone flows investigated with DSMC. The unresolved discrepancy with the measurements could be due to inconsistent free stream conditions with the actual experimental data or missing physical phenomena altogether, for example electronically excited states for both atomic and molecular oxygen, more complex gas–surface interactions, or three-dimensional effects. However, to the best of our knowledge, no bulk flow quantities have been previously reported, such as those shown in figures 14–16. The impact of thermochemical models could be much more significant for portions of the flow away from the cold wall.

The DSMC method naturally accounts for physical phenomena that are often ignored in CFD, namely non-equilibrium between rotational and translational modes, mass diffusion due to thermal gradients in multicomponent mixtures (Soret effect), and non-continuum effects near the wall. For this flow, however, the close agreement with CFD results suggests that these physical phenomena do not significantly influence surface properties. A similar conclusion cannot be reached for the gas flow away from the wall, where more marked differences with standard continuum solvers may be observed.

The advantages of obtaining a DSMC particle solution for these flows resides in the method's ability to be easily informed from first principles at the level of molecular cross-sections and to seamlessly describe internal energy non-equilibrium among all modes. With the advent of exascale computing and beyond, particle methods will be an increasingly important tool to verify the validity of physical assumptions in reduced-order models via experimental-scale simulations, down to the level of molecular-level distributions.