2.1. Viscosity cross-section model

The variable hard-sphere (VHS) cross-section model (Bird Reference Bird1994) 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. (Reference Valentini, Verhoff, Grover and Bisek2023) 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 Reference Bird1994), 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 Reference Lumpkin, Haas and Boyd1991; Haas et al. Reference Haas, Hash, Bird, III and Hassan1994; Gimelshein, Gimelshein & Levin Reference Gimelshein, Gimelshein and Levin2002a; Zhang & Schwartzentruber Reference Zhang and Schwartzentruber2013). 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 Reference Borgnakke and Larsen1975) 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 Reference Millikan and White1963), 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 Reference Bird1994) as it was shown to be relatively independent of temperature in previous MD studies (Valentini, Zhang & Schwartzentruber Reference Valentini, Zhang and Schwartzentruber2012; Kosyanchuk & Yakunchikov Reference Kosyanchuk and Yakunchikov2021). For vibration–translation energy exchange, the vibrational relaxation number $Z_v$ was modelled according to Bird (Reference Bird1994)

(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 Reference Bird1994) involved in the collision. As discussed by Eckert & Gallis (Reference Eckert and Gallis2022), (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. (Reference Valentini, Schwartzentruber, Bender, Nompelis and Candler2015) 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 (Reference Streicher, Krish and Hanson2020), the Millikan–White correlation (Millikan & White Reference Millikan and White1963), and the computational results of Grover, Torres & Schwartzentruber (Reference Grover, Torres and Schwartzentruber2019). For O + O$_2$, the model parameters in (2.3) were chosen to reproduce the recent ab initio computations of Grover et al. (Reference Grover, Torres and Schwartzentruber2019). 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 Reference Bird1994) is used in this work. The TCE reaction cross-section, specialized to the VHS model, is given by Bird (Reference Bird1994)

(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 (Reference Bird1994), $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 Reference Bird1994, 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 Reference Bird1994, 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. Reference Grover, Torres and Schwartzentruber2019) 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.1. Grid-independence study

The need for a rigorous grid-independence analysis has been recognized in various previous studies. For example, Nompelis et al. (Reference Nompelis, Candler and Holden2003) 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. Reference Wang, Boyd, Candler and Nompelis2001) 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. (Reference Nompelis, Candler and Holden2003) or Ninni et al. (Reference Ninni, Bonelli, Colonna and Pascazio2022) 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. Reference Plimpton, Moore, Borner, Stagg, Koehler, Torczynski and Gallis2019). 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$.