Low-frequency molecular fluctuations in the translational nonequilibrium zone of one-dimensional strong shock waves are characterized for the first time in a kinetic collisional framework in the Mach number range 2M10. Our analysis draws upon the well-known bimodal nature of the probability density function (PDF) of gas particles in the shock, as opposed to their Maxwellian distribution in the freestream, the latter exhibiting two orders of magnitude higher dominant frequencies than the former. Inside the (finite-thickness) shock region, the strong correlation between perturbations in the bimodal PDF and fluctuations in the normal stress suggests introducing a novel two-bin model to describe the reduced-order dynamics of a large number of collision interactions of gas particles. Our model correctly predicts two orders of magnitude differences in fluctuation frequencies in the shock vs those in the freestream and is consistent with the small-amplitude fluctuations obtained from the highly resolved direct simulation Monte Carlo computations of the same configuration. The variation of low-frequency fluctuations with changes in the conditions upstream of the shock revealed that these fluctuations can be described by a Strouhal number, based on the bulk velocity upstream of the shock and the shock-thickness based on the maximum density gradient inside the shock, that remains practically independent of Mach number in the range examined.

It is well known that numerical solutions of the deterministic Euler and Navier–Stokes equations (NSE) require special treatment if one is to obtain reliable flow features containing shocks. Solutions of differential equations describing fluid flow motion, in either a shock-capturing1,2 or a shock-fitting context,3,4 can predict both location and unsteadiness of shocks accurately, but the description of the internal shock structure and predictions of the gas properties inside a shock must be modeled in the context of such solutions. Information about the interior of shocks may be obtained from the numerical solution of systems of equations that are more general than the Navier–Stokes equations. The first documentation of differences in shock structure predicted by the Navier–Stokes equations and solution of the more general Bhatnagar–Gross–Krook model for the Boltzmann equation of transport is due to Liepmann et al.,5 who showed that the Navier–Stokes equations are inadequate in describing the internal structure of a shock wave at high Mach numbers (M3) and that the differences between the two solutions increase in the low-pressure region of the shock layer with an increase in Mach number. Soon afterward, Bird6 solved the exact Boltzmann equation using the stochastic DSMC method7 and quantified the level of strong translational nonequilibrium in the interior of shocks.

The significance of the correct description of not only the location and motion of shocks but also of their internal dynamical structure cannot be overstated in fluid mechanics. Among other classes of flows, shocks are responsible for the production of sound in supersonic jets as shown by the extensive studies of their interaction with turbulence and isolated vortices.8–20 The shock-wave/boundary-layer interactions (SBLIs) on compression ramps, cones, and flat plates are widely investigated for their role in generating separation bubbles, unsteadiness, and large surface heat and pressure fluxes.21–23 It is also recognized that the study of receptivity of shocks to freestream or induced disturbances is of importance in the investigation of transition in hypersonic boundary layers.24–28 

However, the vast majority of research in laminar to turbulent transition and SBLIs focuses on the physics of transition in the boundary layer and not the shock.29–33 In fact, research on the effect of kinetic fluctuations on triggering transition from laminar to turbulent flows has been explored using Landau–Lifshitz's theory of fluctuating hydrodynamics34 in several works on the receptivity of boundary layers in the incompressible35,36 and high-speed compressible37,38 regimes. In the phenomenological approach of fluctuating hydrodynamics, stochastic fluxes are added to the stress tensor and heat flux vector in the NSE formulation to account for the effect of molecular fluctuations on the flow field. Using this approach, Luchini36 showed that the presence of kinetic fluctuations is a non-negligible source of disturbances in a boundary layer and sets an upper bound on the transition Reynolds number. Fedorov and Tumin37 showed that in a compressible flat-plate boundary layer, the kinetic fluctuations could trigger random wave packets of Tollmien–Schlichting waves or a Mack second-mode instability. However, this approach cannot be easily extended to account for the kinetic fluctuations in shocks, where the fluctuation amplitude is known to be larger than the predictions of the equilibrium theory.39 

Therefore, to understand the role of molecular fluctuations in the shock layer, kinetic methods that provide insight into the molecular nature of the shock in terms of particle velocities and energy distribution function are important. However, understanding detailed molecular collisions in two- and three-dimensional (2D and 3D) flows is challenging because of the tremendous length and time scales related to boundary layers, instabilities, and unsteadiness of SBLIs. Our study aims at closing this theoretical gap by examining the origin of molecular fluctuations in the well-known one-dimensional (1D) shock of a monatomic gas, argon, in a kinetic framework. We show the surprising result that even simple one-dimensional shocks exhibit low-frequency fluctuations that have been, up until now, ignored, yet their presence may well contribute to laminar to turbulent transition in supersonic and hypersonic flows.

Historically, solutions of the internal structure of strong shocks were sought using kinetic models5 as they allow for anisotropy of stresses and heat fluxes, which are not accounted for by the traditional Navier–Stokes–Fourier constitutive relations. Bird6 successfully modeled 1D argon shock structures using the DSMC method and obtained shock-thicknesses that were in good agreement with the experiments of Alsmeyer.40 Since then, the method has been widely used for modeling the steady-state behavior of strong shocks41–45 as well as unsteady flows without46–50 and with shocks.51 Stefanov et al.39 simulated Mach 26 bow shock over a 2D cylinder using DSMC and attributed an increase in velocity fluctuations inside the shock layer to “thermal nonequilibrium effects”; however, they did not delve deeper into the reason for these effects, the changes in timescales of fluctuations in comparison with an equilibrium state, and their dependence on the strength of the shock wave. Our work attempts to answer these questions by investigating the fluctuations in molecular velocity and energy distribution functions obtained from the solution of the Boltzmann equation. This work will show that the major role of bimodality in the distribution function inside the shock is to change the dominant frequencies of molecular fluctuations compared to an equilibrium freestream.

Our past work has exploited the fidelity of DSMC to understand the critical role of shocks in hypersonic SBLIs. Tumuklu et al.52,53 simulated laminar SBLIs in a Mach 16 axisymmetric flow over a double-cone and observed strong coupling between the complex shock structure and the separation bubble. For the freestream unit Reynolds number of Re1=3.74×105 m−1, they found oscillations of the detached (bow) and separation shocks characterized by a Strouhal number of 0.078. In this work, we will show that the Strouhal number associated with the low-frequency fluctuations in an isolated 1D shock falls within a similar range of St= 0.001–0.02. Sawant et al.54 extended the capability of DSMC in modeling the computationally challenging 3D laminar SBLIs of a Mach 7 flow over a double-wedge at near-continuum input conditions through the use of adaptively refined octree grids. In addition, we have investigated the stability of the spanwise-periodic laminar separation bubble in the double-wedge flow to self-excited, small-amplitude, spanwise homogeneous perturbations.55,56 We demonstrate the existence of the strong coupling of the shock structure with the separation bubble and show that linear amplification of the most unstable 3D flow perturbations leads to synchronization of low-frequency unsteadiness of the separation bubble with the unsteadiness of the manifold of shock structures corresponding to a Strouhal number of 0.028.

The remainder of the paper is organized as follows: Sec. II describes the DSMC numerical implementation and the nonequilibrium aspects of a 1D shock structure of argon. It then describes the fluctuations in the overall stress inside the shock with that in the freestream and contrasts the differences in their frequencies. From additional analyses of cross-correlations, it is hypothesized that the differences in fluctuations are caused by the long-time collision interaction of particles in two modes of the bimodal energy distribution of particles. To prove this hypothesis, we construct a simplified two-energy-bin ordinary differential equation (ODE) model in Sec. III similar to the predator–prey model of Lotka–Volterra57–59 but with modified terms accounting for intermolecular collisions between the two bins. The evaluation of rate coefficients used in the model is also described in this section, and the simplification of the collision processes is discussed in detail in  Appendix B. In Sec. IV, we discuss the dynamics of the two-energy-bin model and compare with the broadband frequencies obtained from the power spectral density (PSD) analyses of fluctuations in the DSMC residuals. Section V establishes a range of Strouhal numbers for Mach numbers 2–10 as well as for variations in upstream temperature at a given Mach number. Using this range, an estimate of low frequencies is made and qualitatively compared with experiments. Finally, Sec. VI summarizes the present findings.

The DSMC solution of the shock is obtained using a 1D version of the Scalable Unstructured Gas-dynamic Adaptive mesh-Refinement (SUGAR) DSMC solver,54 which makes use of a binary adaptive mesh refinement structure, where each computational Cartesian “root” cell is recursively refined into smaller cells until its size is smaller than the local mean-free-path. The smallest cells, referred to as “leaf” or “collision” cells, are used to select neighboring collision partners to perform binary elastic collisions between particles of a monatomic gas using the majorant frequency scheme.60 The macroscopic flow and transport parameters are computed based on the statistical equations of kinetic theory7 and shown on larger root cells. In terms of traditional DSMC solvers, each reference to a computational cell in this paper means a Cartesian root cell. The gas is assumed to follow a variable hard sphere (VHS) molecular model7 with viscosity index ω=0.81, mass m=6.637×1026 kg, and reference diameter dr=4.17×1010 m at reference temperature Tr = 273 K. We compared our DSMC profiles of normalized viscous stress and heat flux along the direction of shock propagation with the results of Ohwada61 obtained by using a finite-difference Boltzmann solver for a 1D, Mach 3 flow of argon and obtained excellent agreement.

The DSMC simulation was initialized using the Rankine–Hugoniot jump conditions at X/λ1=0, where X is the streamwise direction normal to the shock and λ1 is the upstream mean-free-path. At the upstream boundary, X/λ1< 0, computational particles are introduced using a directed local Maxwellian flow at a number density, bulk velocity, and temperature of n1=1×1022 m−3, ux,1=3572.24 m s−1, and Ttr,1=710 K, respectively, to generate hypersonic flow. At the subsonic downstream boundary, X/λ1> 0, a specularly reflecting surface moving at a bulk velocity of ux,2, equal to the downstream Rankine–Hugoniot velocity is used (see Ref. 7, Chap.12). The domain length is 0.02 m, that is, X/λ1=57.45 to 57.45, which is sufficiently large so that the boundaries are placed sufficiently far from the shock layer but small enough so that the shock does not move from its initialized central location during the early transient period. In this domain, the flow takes only ∼6 μs to transition from the jump condition to form a typical shock structure, a time much shorter than the low-frequency molecular fluctuations that we will discuss. Bird's STABIL boundary condition7 is used to prevent random walk effects from causing the shock to move; however, because of the high number density and the use of a large number of computational particles (475 000), the drift was found to be very small. Nevertheless, the STABIL algorithm has not been previously used in a 1D shock simulation where we wish to characterize molecular fluctuations that generate low-frequency unsteadiness. For this reason, additional clarification is provided in  Appendix A to demonstrate that the frequencies that we will observe from DSMC in Secs. II B and V are not influenced by the STABIL algorithm.

Using the aforementioned computational approach, Fig. 1 shows time-averaged profiles of macroscopic flow and transport parameters obtained from the DSMC simulation of a 1D, Mach 7.2 shock in argon. Unless otherwise indicated, subsequent figures and results are based on these 1D shock conditions. The time-averaged or mean quantities discussed in this work were calculated after the 6-μs transient period for additional time steps until ∼0.2 ms. Figure 1(a) shows deviation of the X and average lateral directional temperatures, Ttr,x and based on axial symmetry, Ttr,n=(Ttr,y+Ttr,z)/2, respectively, from the overall translational temperature, Ttr, within a region of 16<X/λ1<6, indicating the presence of translational nonequilibrium in the shock. In the equilibrium regions (X/λ1<16 and X/λ1>6), it can be seen that all directional temperatures are equal to the overall temperature. Within the region 2<X/λ1<6, the gradient of overall translational temperature asymptotes to zero, yet there is a strong deviation of directional temperatures from each other, similar to that in the DSMC simulation of a Mach 8 argon shock by Bird.6Figure 1 also shows the locations of the numerical probes F (X/λ1=29), P (X/λ1=1), and D (X/λ1=29), based on which later discussions will follow.

FIG. 1.

(a) Ttr,x,Ttr,n, and Ttr normalized by the upstream translational temperature Ttr,1. (b) Normalized p, σxx, τxx, and τnn (see the text). (c) and (d) show normalized streamwise bulk velocity, (uxux,1)/(ux,2ux,1), and number density, (nn1)/(n2n1), respectively, as well as the comparison of DSMC-derived standard deviation in their fluctuations, σux and σn, in percent of respective unnormalized macroscopic quantities with estimates from equilibrium statistical mechanics. Subscripts “1” and “2” denote the upstream and downstream macroscopic quantities, respectively.

FIG. 1.

(a) Ttr,x,Ttr,n, and Ttr normalized by the upstream translational temperature Ttr,1. (b) Normalized p, σxx, τxx, and τnn (see the text). (c) and (d) show normalized streamwise bulk velocity, (uxux,1)/(ux,2ux,1), and number density, (nn1)/(n2n1), respectively, as well as the comparison of DSMC-derived standard deviation in their fluctuations, σux and σn, in percent of respective unnormalized macroscopic quantities with estimates from equilibrium statistical mechanics. Subscripts “1” and “2” denote the upstream and downstream macroscopic quantities, respectively.

Close modal

In addition, Fig. 1(b) shows the deviation of time-averaged pressure, p, from the X-directional overall stress, σxx, inside the shock as well as the normal viscous stress components, τxx and τnn=(τyy+τzz)/2, calculated based on the relationship

τij=(σijpδij),
(1)

where δij is the Kronecker delta function. The pressure and stress components are normalized by the upstream parameter ρ1β12, where β=m/2κbTtr is the inverse of the most probable speed of molecules, m is mass, κb is the Boltzmann constant, ρ=nm is the mass density, and n is the number density. The existence of these non-zero stresses leads to a finite-thickness of the shock wave, which, when modeled by accounting for molecular thermal fluctuations using a stochastic method, such as DSMC, reveals nonequilibrium bimodal velocity and energy distributions of particles. This paper will show that these bimodal distributions exhibit a dominant frequency that is two orders of magnitude lower than those found in the freestream.

With respect to the equilibrium regions of the shock, Figs. 1(c) and 1(d) show excellent agreement between the DSMC-computed velocity and number density fluctuations and theory, as was observed by Stefanov et al.39 For a dilute gas, theory predicts that at equilibrium, statistical fluctuations in macroscopic quantities of bulk velocity, ux, and number of particles, N, have a standard deviation of κbTtr/mN and N, respectively, where the brackets denote time or ensemble average of macroscopic quantities.34,62 The DSMC standard deviations are computed for each computational cell of width Δx=1×104 m for instantaneous data collected at every time step of Δt=3 ns from 6 μs to 0.2 ms. Note that since the velocity fluctuations shown in Fig. 1(c) are expressed in terms of DSMC computational particles, the amplitude of actual thermal fluctuations is obtained by division by a factor of Fn, where Fn=107 is the number of dilute gas molecules represented by each DSMC computational particle. Therefore, the 0.344% standard deviation of velocity fluctuations in the freestream in Fig. 1(c) corresponds to an actual value of 1.09×104 %, a small yet significantly large number on the scale of small amplitude perturbations considered in shock-dominated flows.

We are particularly interested in the nonequilibrium zone of the shock layer, where a significant deviation from equilibrium is seen in DSMC-computed standard deviations, as shown in Figs. 1(c) and 1(d). Note that the standard deviations peak at the location of absolute maximum gradients of respective flow parameters, that is, at location X/λ1=0. In contrast to the findings of Stefanov et al.,39 the density fluctuations inside a shock also deviate from the equilibrium Poisson law, although their magnitude is much smaller than the velocity or energy fluctuations.

To understand the frequencies of fluctuations in the shock layer, we focus on probe P located in the vicinity of maximum absolute bulk velocity gradient in the shock to obtain the time history of fluctuations of the normalized X-directional overall stress, σxxβ12ρ11, about the time-averaged mean. These data are important because fluctuations in σxx correspond to fluctuations in Ttr,x given that σxx=ρRTtr,x and density fluctuations are negligible. The original instantaneous signal is window-averaged with a moving time-window of 0.3 μs (100 timesteps) such that frequencies greater than 3.33 MHz are removed. The power spectral density (PSD) of this mean-subtracted, window-averaged signal reveals frequency content in fluctuations of σxx, as shown in Fig. 2. For spectral estimation here and elsewhere in the paper, Welch's method63,64 is used in SciPy65 software with two Hann-window weighted segments of data sampled with a frequency of 333 MHz prior to the Fast Fourier Transform (FFT) such that the frequency resolution is 0.9 kHz. In  Appendix A, the sensitivity of the PSD obtained using DSMC particle data to the value of Fn is discussed and it is shown that when Fn is set to unity (i.e., no possible numerical amplification of noise), the main features of the non-dimensionalized low-frequency spectrum do not change.

FIG. 2.

PSD of mean-subtracted, window-averaged σxx is shown inside the shock layer at probe P in (a), in the entire X/λ1 space in (b), and in the freestream at probe F in (c). (a) and (c) also show NCE (see the text).

FIG. 2.

PSD of mean-subtracted, window-averaged σxx is shown inside the shock layer at probe P in (a), in the entire X/λ1 space in (b), and in the freestream at probe F in (c). (a) and (c) also show NCE (see the text).

Close modal

At probe P, the PSD shows a broadband of low frequencies that goes up to 90 kHz frequency, as shown in Fig. 2(a). This upper bound of the broadband is defined as the frequency at which the normalized cumulative energy (NCE), obtained from normalizing and cumulatively summing the PSD spectrum, exhibits an inflection point. This low-frequency broadband is defined by its weighted average of 37.5 kHz and a standard deviation of 21.4 kHz. Figure 2(b) shows the contours of PSD plotted on the axes X/λ1 vs frequency created by interpolating the PSD data at each X/λ1 location spatially separated by X/λ1=1 in the Tecplot-36066 software using the inverse-distance algorithm with default parameters (exponent = 3.5, point selection=Octant, Number of points = 8). The figure reveals that such low-frequency broadband is observed within 6<X/λ1<3, the region of strong translational nonequilibrium seen earlier in Fig. 1.

In contrast to the nonequilibrium region of the shock, Fig. 2(c) shows that the PSD at probe F in the freestream contains widely distributed energy across the whole spectral limit of 166 MHz and does not exhibit a noticeable inflection point within this limit. It also shows that most of the spectral energy is concentrated near the frequency of 3.5 MHz, which is the inverse of the mean-collision-time of 0.284 μs. This is consistent with the analysis of Kogan (see Ref. 67, pp. 82 and 124) for a gas in absolute equilibrium (zero bulk velocity, constant temperature, and number of particles), which showed that small perturbations to the velocity distribution function and its moments decay exponentially with a characteristic timescale of mean collision time. In addition, from the comparison of PSDs in Figs. 2(a) and 2(c), it can be noted that the low-frequency fluctuations in the shock have about two orders of magnitude larger amplitude than the freestream, and therefore, low-frequency fluctuations of shock may be even more important in laminar to turbulent transition of high-speed boundary layers than, as commonly believed, freestream kinetic fluctuations. The order of magnitude differences in these timescales of fluctuations will be explained by the new dynamics model discussed in Sec. IV.

The key to understanding the fluctuations in the overall stress, σxx, is to correlate the fluctuations in the PDF fξx of the normalized X-directional energy of particles, ξx, and its mean, μ. Note that ξx=(vx2ux2)β2, where vx is the X-directional instantaneous molecular velocity of particles, which is the sum of their bulk and thermal components, ux and cx, respectively. The mean, μ, is related to overall stress, σxx, as

μ=ξxfξxdξx=σxxβ2ρ.

Figure 3 shows that the behavior of the PDF changes at different locations in the flow, that is, from the upstream to downstream regions relative to the shock, fξx changes from a nearly symmetric equilibrium distribution [Fig. 3(a)] to a one-sided, asymmetric equilibrium distribution [Fig. 3(b)]. At the location of maximum density gradient in the shock [Fig. 3(c)], it can be seen that the density function is bimodal with an inflection point at ξx=1.33, where the portions 0.86<ξx<1.33 and 1.33<ξx<3.52 have shapes similar to PDFs in the subsonic and hypersonic regions, shown in Figs. 3(b) and 3(a), respectively. Note that the bimodal PDF can be expressed as a linear combination of upstream and downstream contributions of equilibrium PDFs,68 similar to the Mott–Smith model of the bimodal velocity distribution.69 

FIG. 3.

The PDF fξx at probes (a) F in the freestream, (b) D downstream, and (c) P inside the shock layer.

FIG. 3.

The PDF fξx at probes (a) F in the freestream, (b) D downstream, and (c) P inside the shock layer.

Close modal

Toward that end, the cross-correlation coefficient of μ and the number of particles as a function of ξx are defined as

c(ξx)=w=0w=Wξx=minmax[N(ξx,w)N(ξx)w][μ(w)μw]WΣN(ξx)Σμξx,
(2)

where

ΣN(ξx)=[N(ξx,w)N(ξx)w]2W,Σμ=[μ(w)μw]2W.

Note that ξx is discretized into 200 energy bins from its minimum to maximum value. N(ξx,w) is the total number of DSMC particles within the normalized X-directional energy space ξx and ξx+Δξx from time window w to w +1, and N(ξx)w denotes the number of particles in the same energy space but averaged over time-windows, W. Similarly, μ(w) is the instantaneous mean of the PDF fξx, from time-window w to w +1, whereas μw is the mean computed by averaging over all time-windows. ΣN(ξx) and Σμ are the standard deviations in the fluctuations of N(ξx,w) and μ(w) about their respective means. Figures 4(a) and 4(b) show the fluctuations in number of particles [N(ξx,w)N(ξx)w] in the energy space ξx at probes P and F, respectively.70 

FIG. 4.

The fluctuations [N(ξx,w)N(ξx)w] as a function of ξx at probes (a) P inside the shock layer and (b) F in the freestream. Multimedia views: https://doi.org/10.1063/5.0065971.1; https://doi.org/10.1063/5.0065971.2

FIG. 4.

The fluctuations [N(ξx,w)N(ξx)w] as a function of ξx at probes (a) P inside the shock layer and (b) F in the freestream. Multimedia views: https://doi.org/10.1063/5.0065971.1; https://doi.org/10.1063/5.0065971.2

Close modal

The analysis of the correlation functions combined with particle distributions suggests an approach for grouping of particles into a finite number of energy bins whose dynamics can then be further analyzed. Starting with Fig. 5, the distribution of c(ξx) at probes P, inside a shock, and F, in the freestream, is compared. c(ξx) is calculated using a total number of time windows, W, of 646, each representing a time interval of 0.3 μs starting from t =6 μs. At probe P, Fig. 5(a) shows that c(ξx) has a strong negative correlation between 0.86<ξx<1.33 and positive correlation for 1.33<ξx<3.52. Opposite signs of c(ξx) means that particles in these parts lose and gain particles, respectively, when μ and σxx increase and vice versa, and understanding the underlying mechanisms in these regions, such as intermolecular collisions and flux exchange with adjacent cells, is necessary to understand fluctuations in σxx. These two zones in the energy space are denoted as energy bins “A” and “B,” which are coarser than the size of the histogram bins, Δξx=0.0115, used to obtain smooth variations of fξx. We ignore the energy space beyond ξx>3.52, because c(ξx) remains negative but small in comparison with other two bins and the total average number of particles in this region is only 3.5% of total particles at probe P, which can be calculated from the cumulative sum of the distribution of average number of particles shown in Fig. 6(a). The demarcation between the coarse energy bins A and B is at c(ξx)=0, which is also the location of inflection point. Examination of c(ξx) for probe F in the freestream [see Fig. 5(b)] shows that, in contrast, the maximum magnitude of c(ξx) is not more than 0.2, indicating that the fluctuations are random in nature and even if a strong correlation exists, it must be on timescales smaller than the window size of 0.3 μs used to obtain c(ξx). Comparison of the correlation functions at probes P and F also shows that the distribution of c(ξx) at the latter location is almost symmetric about the inflection point at ξx=0, which is consistent with the X-directional energy distribution having a nearly symmetric shape, as shown in Fig. 3(a).

FIG. 5.

(a) and (b) show c(ξx) calculated using Eq. (2), at probes P and F, respectively. Here and in subsequent figures, the bracket notation, Qt, denotes the time-average of a macroscopic quantity Q, in contrast to its instantaneous value. (c) and (d) show the ratios of instantaneous to time-averaged number of particles in bins A and B, that is, NA/NAt and NB/NBt, at probes P and F, respectively. To reduce statistical scatter, the instantaneous number of particles is averaged over a very small time interval of 0.3 μs, whereas the time-averaged number of particles is averaged from 6 μs to 0.2 ms.

FIG. 5.

(a) and (b) show c(ξx) calculated using Eq. (2), at probes P and F, respectively. Here and in subsequent figures, the bracket notation, Qt, denotes the time-average of a macroscopic quantity Q, in contrast to its instantaneous value. (c) and (d) show the ratios of instantaneous to time-averaged number of particles in bins A and B, that is, NA/NAt and NB/NBt, at probes P and F, respectively. To reduce statistical scatter, the instantaneous number of particles is averaged over a very small time interval of 0.3 μs, whereas the time-averaged number of particles is averaged from 6 μs to 0.2 ms.

Close modal
FIG. 6.

(a) and (b) show a comparison of the average number particles as a function of X-directional energy, N(ξx)t, at probes P and F, respectively, with those in an adjacent cell to their left. The time averages are taken from t = 6 μs to 0.2 ms.

FIG. 6.

(a) and (b) show a comparison of the average number particles as a function of X-directional energy, N(ξx)t, at probes P and F, respectively, with those in an adjacent cell to their left. The time averages are taken from t = 6 μs to 0.2 ms.

Close modal

If we now group particles into two X-directional energy bins, A and B, which are denoted in Figs. 5(a) and 5(b) for probes P and F, respectively, we can study the ratios of the instantaneous to time averaged particles, NA/NAt and NB/NBt as a function of time at these two probes, as shown in Figs. 5(c) and 5(d). The average number of DSMC computational particles per computational cell in energy bins A and B is NAt=1088 and NBt=611 at probe P, and NAt=510 and NBt=493 at probe F, which can be converted to number density (m3) by multiplication of an Fn=1×107 and division by the volume of computational cell, Δv=1×1012m3. From the comparison of these two figures, it can be seen that fluctuations in these ratios are larger in magnitude and exhibit longer timescales (0.01–0.05 ms) at probe P inside the shock in comparison to probe F in the freestream, where fluctuations appear to be random as their timescales are on the order of mean collision time of 0.284 μs, smaller than the time interval used to time-average the signal (0.3 μs). The standard deviations in the fluctuations of the ratios NA/NAt and NB/NBt are 4.5 and 2.1%, respectively, at probe P, vs 1.8 and 1.6% at probe F. In addition, a noticeable negative correlation between the fluctuations (out-of-phase fluctuations) of two energy bins at probe P can also be seen from Fig. 5(c), while such dynamics are absent at probe F.

Finally, to demonstrate the generality of the difference in properties between probes P and F, Fig. 6 shows a comparison of the time-averaged number of particles as a function of ξx at probes P and F, and at their respective left-adjacent computational cells. For probe P [Fig. 6(a)], two inflection points can be identified at ξx=1.33 and 3.52, where c(ξx) changes sign. We can also see that the average number of particles in energy bin A is larger at probe P than in the cell to its left, whereas in bin B it is lower. This is consistent with the fact that from the upstream to downstream region inside a shock, the contribution of the upstream symmetric distribution decreases, while the downstream asymmetric distribution increases. At probe F [Fig. 6(b)], however, the PDF of energies does not change between adjacent cells and there is no inflection point.

Although the distribution of X-directional particle energies is continuous, we propose a two-bin model to understand the role of particle collisions between the two bins, A and B, and how they induce low-frequency time dynamics in fluctuations of macroscopic flow parameters such as normalized overall stress, σxx, inside a shock. We first develop a simple two energy-bin model similar to the Lotka–Volterra's predator–prey model.57–59 We then evaluate the collisions rate coefficients for different energy transfer processes using the DSMC particle distribution data and show that instead of sustained particle oscillations, the solution to our ODE contains dampening terms that cause the periodic fluctuations to die out on timescales an order of magnitude longer than the period of oscillation.

The dynamics of the number of particles in energy bins A and B may be written as

dNAdt=[dNAdt]coll+FANA,dNBdt=[dNBdt]coll+FBNB,
(3)

where the first and second terms on the right are the change in number of particles due to collisions and fluxes, respectively. In addition, since the average number of particles in energy bins A and B constitutes 96.5% of the total particles, we can assume them to be mutually exclusive and write

[dNAdt]coll=[dNBdt]coll.
(4)

The flux coefficient Fj for bin “j” is defined as the difference between the net average influx from the left boundary and outflux from the right boundary of the “j”-type particles per second divided by the average number of “j”-type particles,

Fj=(Njl,inNjl,out)(Njr,outNjr,in)ΔtNjt,j{A,B},
(5)

where Δt is the DSMC time step, superscripts “l” and “r” refer to the left and right boundaries of the computational cell, respectively, and “in” and “out” refer to the incoming and outgoing particles, respectively, as shown in Fig. 7 for location P and energy bins A and B.

FIG. 7.

A sketch denoting number of particles coming in and going out through the left and right boundary in a two-energy-bin model.

FIG. 7.

A sketch denoting number of particles coming in and going out through the left and right boundary in a two-energy-bin model.

Close modal

Flux coefficients are required in the evaluation of Eq. (3). Using the values of A- and B-type particles given in Table I, flux values of FA and FB = −1 277 925.8 and 2 500 681.6 s–1, respectively, at probe P and zero in the freestream (as expected) are obtained. To make sure that the predictions of the model are not affected by instantaneous fluctuations or their numerical amplification, the DSMC-derived fluxes provided in Table I are temporally averaged from 6 μs to 0.21 ms, a time period ∼8 times longer than the timescale corresponding to the average low-frequency observed in the PSD. Note that at probe P, FA<0 because on average more A-type particles travel to a cell downstream than came in from the upstream, as also seen from Fig. 6(a). The opposite is true for type-B particles, resulting in FB>0. This is synonymous with the fact that if there were no collisions, that is, [dNAdt]coll=0, the number of A particles, those that mainly represent the subsonic part of the bimodal distribution, would decrease, while the number of B particles, those that mainly represent the upstream hypersonic flow, would increase.

TABLE I.

Average fluxes of DSMC computational particles per time step and computational cell area of 1×108 m2.

ProbesNAl,inNAl,outNAr,outNAr,inNBl,inNBl,outNBr,outNBr,in
P 41.46 7.19 47.40 8.96 69.84 0.68 65.37 0.80 
F 49.28 49.28 0.0 58.63 0.0 58.63 0.0 
ProbesNAl,inNAl,outNAr,outNAr,inNBl,inNBl,outNBr,outNBr,in
P 41.46 7.19 47.40 8.96 69.84 0.68 65.37 0.80 
F 49.28 49.28 0.0 58.63 0.0 58.63 0.0 

To evaluate the collision terms in Eq. (3), we need to identify all possible collision processes that can cause a loss or gain of type-A and type-B particles. Note that different collision processes can occur at different rates, which could change the dynamics of particle exchange between the two bins and consequently the timescale of fluctuations in σxx. Such collision processes are listed in column two of Table II, where the type of A particle is labeled based on the collision process that it undergoes, denoted by a subscript “i,” the relevant pre-and post-collisional states are specified where the latter is denoted by a primed superscript, and the rate coefficient for each fundamental collision process, ki, unless it is a “compound” rate, which is denoted by a tilde. The table shows that there are six fundamental processes (second column) and eight types of A particles, Aa,Ab,Ac,Ac,Ad,Ad,Ae, and Af, that can cause a net change in the number of particles in energy bins A and B. For example, Ab particles are the postcollisional particles located in energy bin A that are formed due to the collision process b, while Ac are the precollisional particles in energy bin A that take part in collision process c. Using the DSMC-derived PDFs for all sub-types of type-A particles of the normalized X-directional energy, ξx=(vx2ux2)β2, and total transverse energy, ξy+ξz=(vy2+vz2)β2, denoted by fξxA and fξy+ξzA, respectively (Fig. 8), we can reduce the number of types of A particles further, which allows us to write processes pi (column two) in the form of processes Pi (column three). In what follows, we only describe the grouping process for processes a and b, while details for the other processes are given in  Appendix B. Then, the Pi processes are further simplified to processes Qi (column 4), as described in  Appendix B.

TABLE II.

Detailed and simplified collision processes.

IdDetailed collision processCollision process after grouping alike A-particlesSimplified collision process
ipiPiQi
Aa + B ka B + B Aa + B ka B + B A + B k̃a B + B 
B + B kbAb + B B + B kb Aa + B B + B kb A + B 
Ac + B kcAc + Ac Ac + B kcAc + Ac A + B k̃cAc + Ac 
Ad + AdkdAd + B Ac + Ackd Ac + B ⋯ 
B + B keAe + Ae B + B keAe + Ae B + B keAe + Ae 
Af + Afkf B + B Ae + Aekf B + B ⋯ 
ga Ag + B kg B + Ag Ag + B kg B + Ag ⋯ 
ha Ah + AhkhAh + Ah Ah + Ahkh Ah + Ah ⋯ 
ia B + B ki B + B B + B ki B + B ⋯ 
IdDetailed collision processCollision process after grouping alike A-particlesSimplified collision process
ipiPiQi
Aa + B ka B + B Aa + B ka B + B A + B k̃a B + B 
B + B kbAb + B B + B kb Aa + B B + B kb A + B 
Ac + B kcAc + Ac Ac + B kcAc + Ac A + B k̃cAc + Ac 
Ad + AdkdAd + B Ac + Ackd Ac + B ⋯ 
B + B keAe + Ae B + B keAe + Ae B + B keAe + Ae 
Af + Afkf B + B Ae + Aekf B + B ⋯ 
ga Ag + B kg B + Ag Ag + B kg B + Ag ⋯ 
ha Ah + AhkhAh + Ah Ah + Ahkh Ah + Ah ⋯ 
ia B + B ki B + B B + B ki B + B ⋯ 
a

These processes cause no net change in the number of particles of energy bins A and B.

FIG. 8.

(a) and (b) show, at Probe P, the PDFs fξxA and fξy+ξzA of all sub-types of type-A particles distinguished based on detailed collision processes listed in Table II. We note that fξxA is obtained by further refining the X-directional energy space of bin A (0.86<ξx<1.33) with a resolution of Δξx=0.0115. Same resolution was used for the transverse energy space. (c) and (d) show the respective plots for probe F.

FIG. 8.

(a) and (b) show, at Probe P, the PDFs fξxA and fξy+ξzA of all sub-types of type-A particles distinguished based on detailed collision processes listed in Table II. We note that fξxA is obtained by further refining the X-directional energy space of bin A (0.86<ξx<1.33) with a resolution of Δξx=0.0115. Same resolution was used for the transverse energy space. (c) and (d) show the respective plots for probe F.

Close modal

Collision process pa describes a mechanism for type-Aa particles in energy bin A that collide with particles in energy bin B and now belong to bin B resulting in the loss of type-A particles. Simultaneously, process pb describes the collision process of two bin B particles that causes one of them to move to energy bin A leading to a gain of type-A particles, denoted by Ab. However, Figs. 8(a) and 8(b) show that at location P the energy distributions fξxA and fξy+ξzA for Aa and Ab particles are the same. Therefore, we can group Ab- and Aa-type particles and write the second reaction Pb as shown in column three of Table II. The same holds true for probe F as shown in Figs. 8(c) and 8(d); however, the distributions are different from those at probe P.

We can also group all collisions between type-A particles that do not lead to a net change in the energy distribution of particles within bin A as process ph and see from Fig. 8 that the energy distributions of all such type-A particles before and after collisions, Ah and Ah, is the same at probe P as well as F. Note that at probe F, the distribution fξy+ξzA for Aa-type particles matches with the Ah-type, that is, in the freestream the transverse energy of particles taking part in processes Pa and Pb is the same as those in Ph. The differences at probe P indicate the role of transverse translational modes in the collision processes Pa and Pb. Similar to process ph, process pg, is a bin-exchange process, which causes no net change in the number of particles in bins A and B. Particles simply swap energy bins by exchanging X-directional energies. Using similar logic, additional simplifications can be made to the detailed collision processes, pi, to construct the Pi column of Table II based on the energy distribution functions at locations P and F, as described in  Appendix B.

To complete the dynamical model for particles moving between energy bins A and B, an expression is required for the net change in the number of particles in energy bin A due to collisions from all possible processes related to A-type particles, that is,

[dNAdt]coll=[dNAadt]coll+[dNAcdt]coll+[dNAcdt]coll+[dNAedt]coll,
(6)

where NAa,NAc,NAc, and NAe are the number of particles of Aa, Ac, Ac, and Ae types, respectively, which are defined based on Pi processes. Using standard rate equation formulations for elementary reactions (see Ref. 71, Sec. 16.15), the terms on the right hand side can be expressed as

[dNAadt]coll=kaNAaNB+kbNB2,[dNAcdt]coll=kcNAcNB+kdNAc2,[dNAcdt]coll=2kcNAcNB2kdNAc2,[dNAedt]coll=2keNB22kfNAe2.
(7)

By substituting Eq. (7) in (6), we obtain

[dNAdt]coll=kaNAaNB+kbNB2+kcNAcNBkdNAc2+2keNB22kfNAe2.
(8)

The conditions in the 1D shock allow us to make additional simplifications to Eq. (8), as described in  Appendix B, so that using Eqs. (3), (4), and (B9) we obtain the final system of ordinary differential equations used to study perturbation dynamics,

dNAdt=(k̃ck̃a)NANB+(kb+2ke)NB2+FANA,dNBdt=(k̃ck̃a)NANB(kb+2ke)NB2+FBNB.
(9)

Finally, to solve the dynamical system of equations we need to evaluate the rate coefficients in Eq. (9). A summary of the expressions for the rate coefficients used in Eq. (9) and their values are presented in Table III for probes P and F based on the DSMC collision data collected during the simulation, summarized in Table IV. Note that Cit is the average number of DSMC collisions per time step per cell volume that take part in the ith detailed collision process of particle types R1 and R2 [see. Eq. (B1)]. Similar to the flux coefficients, Cit were temporally averaged from 6 μs to 0.21 ms. The rate coefficients k/(s1) shown in Table III can be converted to the traditional units of K/(m3s1) by multiplying them with the volume of computational cell, Δv=1×1012m3, and dividing by the parameter, Fn=1×107. By this conversion, it can be easily shown that the rates K̃a, Kb range from 1.85.0×1017 m3 s−1 at probes P and F.

TABLE III.

Rate coefficients per second for collision processes defined in Table II.

Rate coefficientFormulaProbe PProbe F
k̃a Cat/ΔtNAtNBt 188.9 497.0 
kb Cbt/ΔtNBt2 417.5 513.5 
k̃c (CctCdt)/ΔtNAtNBt 1843.5 
ke Cet/ΔtNBt2 75.0 
Rate coefficientFormulaProbe PProbe F
k̃a Cat/ΔtNAtNBt 188.9 497.0 
kb Cbt/ΔtNBt2 417.5 513.5 
k̃c (CctCdt)/ΔtNAtNBt 1843.5 
ke Cet/ΔtNBt2 75.0 
TABLE IV.

Average number of DSMC collisions per time step per computational cell volume, Cit.

ProbesCatCbtCctCdtCetCftCgtChtCit
P 0.377 0.468 5.074 1.396 0.084 0.0079 2.980 8.849 1.481 
F 0.375 0.374 0.382 0.383 0.0 0.0 1.028 0.893 0.887 
ProbesCatCbtCctCdtCetCftCgtChtCit
P 0.377 0.468 5.074 1.396 0.084 0.0079 2.980 8.849 1.481 
F 0.375 0.374 0.382 0.383 0.0 0.0 1.028 0.893 0.887 

Note that in the freestream, CatCbt, due to detailed balance in a local equilibrium condition, which dictates that each detailed collision must be balanced by its inverse collision process (see Ref. 72, pp. 38, 39, and 42). Furthermore, the net rate, k̃c, is nearly zero in the freestream because again by principle of detailed balance, CctCdt, as shown in Table IV. On the other hand, in the shock, a factor of 3.63 difference between the average values of these two types of collisions results in K̃c=1.83×1016 m3 s−1. Also, from Table IV it can be seen that both Cet and Cft are zero in the freestream, which means that these processes take place only in conditions of translational nonequilibrium. As a result, the rate Ke is zero in the freestream, whereas in the shock it is two orders of magnitude slower, 7.5×1018 m3 s−1, than K̃c.

The rate coefficients given in Table III can be compared with theoretical estimate of a rate coefficient, Kth, for a forward collision process of a quasi-equilibrium gas (see Ref. 72, pp. 213–216), which is defined as

R1+R2kthS1+S2.

The rate of change of R1-type particles per unit volume due to collisions is written as

[dnR1dt]coll=KthnR1nR2,

where nR1 and nR2 are number densities of R1 and R2 type particles, respectively. For a gas following the VHS molecular model, the rate coefficient Kth can be expressed as

Kth=1εσtgrtγ,
(10)

where

σt=πdr2(2κbTrmrgrt){Γ(2.5ω)}1,grt=(8κbTtrtπmr)12,γ=[1Γ(52ω)(Γ(52ω,EaκbTtrt)EaκbTtrtΓ(32ω,EaκbTtrt))],Ea=12mrgrt2,

where σt is the average equilibrium collision cross section for the VHS gas, grt is the time-averaged relative velocity, γ is the fraction of collisions for which the relative translational energy along the line of centers of colliding molecules exceeds the activation energy of Ea [see Ref. 7, Sec. 6.2 and Eq. (4.72))], and Γ is the gamma function. Also, mr is the reduced mass [see Ref. 7, Eq. (2.7)], d is the VHS molecular diameter (see Ref. 7, Eq. (4.63)], and ε is a symmetry factor equal to two because R1 = R2. Table V shows that the estimated parameters in Eq. (10) and average rate coefficients at probes P and F located inside the shock and the freestream, respectively, are consistent with the values presented in Table III.

TABLE V.

Parameters used in Eq. (10) for the theoretical estimate of average rate coefficient.

ParametersAt probe PAt probe F
Ttrt/(K) 10 440 710 
grt/(m s−13326.2 867.6 
σt/(m21.81×1019 4.16 × 10–19 
Ea/(J) 1.83 × 10–19 1.25 × 10–20 
γ 0.351 0.351 
Kth1/(m3 s−11.05 × 10–16 6.33 × 10–17 
ParametersAt probe PAt probe F
Ttrt/(K) 10 440 710 
grt/(m s−13326.2 867.6 
σt/(m21.81×1019 4.16 × 10–19 
Ea/(J) 1.83 × 10–19 1.25 × 10–20 
γ 0.351 0.351 
Kth1/(m3 s−11.05 × 10–16 6.33 × 10–17 

The numerical solution of the ODE system in Eq. (9) at probes P and F is shown in Fig. 9 as a function of time and in Fig. 10 as path lines in the phase space of NA vs NB. Starting with location P, Fig. 9(a) shows that the solution converges to the non-zero critical values of NA,crit and NB,crit, obtained by setting dNA/dt=dNB/dt=0, which are listed in Table VI. There are 18.2 and 7.5% differences between the critical values and the average number of DSMC particles per computational cell volume in bins A and B at probe P, which may be attributed to the simplifications of the two-energy-bin model discussed above. Yet, as will be shown, this simple model reveals two orders of magnitude disparity in frequencies at probes P and F. Figure 10(a) shows the streamlines of a vector field defined by the normalized growth rate vector G=[GNA,GNB]T with components,

GNi=1||G||dNidt,i=AorBand||G||=(dNAdt)2+(dNBdt)2.
(11)
FIG. 9.

(a) and (b) show the solution of the two-energy-bin ODE model at probes P and F, respectively. The initial values of NA and NB are taken to be +13.5 and −6.3% of the critical values, respectively, at probe P, and +5.4 and −4.8% of the critical values, respectively, at probe F. These values are three times the percent standard deviation observed from the DSMC simulation in the ratios of NA/NAt and NB/NBt.

FIG. 9.

(a) and (b) show the solution of the two-energy-bin ODE model at probes P and F, respectively. The initial values of NA and NB are taken to be +13.5 and −6.3% of the critical values, respectively, at probe P, and +5.4 and −4.8% of the critical values, respectively, at probe F. These values are three times the percent standard deviation observed from the DSMC simulation in the ratios of NA/NAt and NB/NBt.

Close modal
FIG. 10.

(a) and (b) show the streamlines of the normalized growth rate vector defined in Eq. (11) along with an overlaid path line of the instantaneous solution (NA, NB) shown in Fig. 9 at probes P and F, respectively. The starting and end points of the path lines are denoted by hollow and filled circles.

FIG. 10.

(a) and (b) show the streamlines of the normalized growth rate vector defined in Eq. (11) along with an overlaid path line of the instantaneous solution (NA, NB) shown in Fig. 9 at probes P and F, respectively. The starting and end points of the path lines are denoted by hollow and filled circles.

Close modal
TABLE VI.

Dynamics of the ODE at probes P and F.

SolutionProbe PProbe F
NA,crit 1286 511.5 
NB,crit 657.2 495 
L1 −281 762 + 1 765 303.3 j 0.0 
L2 −281 762 − 1 765 303.3 j −499 976 
SolutionProbe PProbe F
NA,crit 1286 511.5 
NB,crit 657.2 495 
L1 −281 762 + 1 765 303.3 j 0.0 
L2 −281 762 − 1 765 303.3 j −499 976 

We can see from the formation of a stable spiral that the critical point does not depend on the initial values. The figure shows that the streamlines are tangent to the direction of the maximum growth vector at a given location of (NA, NB) and show eventual decay to the stable critical point of 1245 and 681. The overlaid path line in Fig. 10(a) describes these dynamics for the instantaneous solution shown in Fig. 9(a).

The Jacobian matrix of Eq. (9) is given by

J(NA,NB)=[(k̃ck̃a)NB+FA(k̃ck̃a)NA+2(kb+2ke)NB(k̃ck̃a)NB(k̃ck̃a)NA2(kb+2ke)NB+FB].
(12)

The eigenvalues of matrix J(NA,crit,NB,crit) are listed in Table VI. At probe P, they are complex conjugates with a negative real part in the shock that gives the rate of decay from the initial to critical values. This is consistent with Fig. 9(a), which shows that at 0.01, 0.015, and 0.02 ms, the solution converges to 1%, 0.3%, and 0.07% of the critical values, respectively. These percentage values correspond to an approximate difference between the critical and instantaneous values of the number of A-type particles of 13, four, and one, respectively.

More importantly, these results suggest that inside the shock, the number of particles are subjected to recurrent microscopic thermal fluctuations, which always decay with timescales on the order of 0.01–0.02 ms with corresponding frequencies on the order of 100–50 kHz. These frequencies are close to the weighted average of the low-frequency broadband of 37.5 kHz observed in Fig. 2(a) from the PSD analysis. The instantaneous particles in the energy bins A and B, shown in Fig. 5(c), are a result of superposition of responses to recurrent low-frequency microscopic fluctuations. Therefore, even with the underlying simplifications of the two-energy-bin model, it is able to predict the existence of the low frequencies that originate due to the interactions in the average number of particles between the two modes of the bimodal PDF of particle energies in the presence of translational nonequilibrium in a 1D shock. In addition, the imaginary part of the pair of eigenvalues listed in Table VI at probe P also indicates oscillations with a time period of 2π/Im(L1) (frequency of 281 kHz), which may correspond to fluctuations observed in Fig. 2(a) in the mid-frequency range of ∼100–300 kHz outside of the low-frequency broadband.

Finally, we can relate the two timescales predicted by Eq. (9) in the shock to collision processes Qe and Qb that describe the relaxation of high energy particles in B to low energy bin A (see Table II). The frequency associated with process Qe is on the order fO(keNB), or equivalently fO(KenB=51kHz), where nB=mNBt=6.81×1021 m−3 is the average number density of particles in energy bin B. Similarly, the frequency associated with process Qb is fO(KbnB=284kHz). Therefore, it is seen that the two timescales present in the shock are a direct consequence of the bimodal energy distribution function. Also, note that although the analysis in this work is for a steady shock, these results hold for unsteady shock as well, because the PDF inside the shock remains bimodal.

Turning to probe F in the freestream, we want to demonstrate that the two-bin energy model is able to predict different dynamics than observed in the shock. Since the flux coefficients FA and FB and the rate coefficients k̃c and ke are zero, we obtain

NA,critNB,crit=kbk̃a,
(13)

which is satisfied by an infinite number of critical points corresponding to different freestream conditions. That is, the number densities are different but the temperature is the same since all solutions have the same shape of the PDF shown in Fig. 3(c), and therefore, the same cross-correlation coefficient is shown in Fig. 5(b). In this case, the solution depends on the initial values of NA and NB, and converges to a critical point that satisfies Eq. (13). Figures 9(b) and 10(b) show the solution of Eq. (9) as a function of time for a set of initial values of (NA, NB) and a streamplot along with an overlaid path line of the solution, respectively. In Fig. 10(b), the dependence on the initial values can be clearly seen. For the given conditions, the solution converges to critical values of NA,crit=511.5 and NB,crit=495, which are within 0.5% of the average number of DSMC particles per computational cell volume, of NAt=509.9 and NBt=492.8, and satisfy the ratio of 1.033 given by Eq. (13).

For the freestream conditions, there are two real eigenvalues, as shown in Table VI with the zero eigenvalue indicating again a non-unique solution of Eq. (13). The negative eigenvalue of -499 976 gives the linear rate of decay, the inverse of which gives the characteristic timescale of 2 μs. This fast decay rate is consistent with the PSD at probe F, discussed in Sec. II, which showed that ∼94% of the total spectral energy is distributed in frequencies greater than 1 MHz. More importantly, in the freestream the model predicts the absence of two orders of magnitude lower dominant frequencies having significantly large energy, which appear to be unique to the region of strong nonequilibrium inside a shock.

With the observation of dominant low-frequency perturbations of macroscopic flow parameters inside the shock, a natural question arises as to whether one can define a nondimensional Strouhal number,

St=fLsux,1,
(14)

that would be constant for different shock strengths. We propose to define the timescale as that required for the flow to traverse a distance equal to the shock-thickness with the upstream bulk velocity. To justify this timescale, cases ranging from Mach two to 10 were run by varying the upstream bulk velocity but keeping the upstream temperature constant (Ttr,1=710 K). Starting with the selection of frequency, f, the PSD of the instantaneous pressure spectrum at the location of the maximum gradient of bulk velocity in the shock for different Mach numbers is shown in Fig. 11(a). A broadband of low frequencies is seen, the boundary of which is defined by the inflection point in the NCE, shown in Fig. 11(b).

FIG. 11.

(a) The contours of PSD obtained by interpolating the pressure spectrum at the center of the shock at each Mach number. For interpolation, the inverse-distance algorithm in the Tecplot-36066 software is used with default parameters (exponent = 3.5, point selection = Octant, Number of points = 8). The overlaid white solid and dashed lines show the upper limits and the weighted-averages with ±1 standard deviation of the low-frequency broadbands for Mach 2–10, respectively. (b) NCE of the PSD for Mach 2–10 used to define the upper limit of the low-frequency broadband.

FIG. 11.

(a) The contours of PSD obtained by interpolating the pressure spectrum at the center of the shock at each Mach number. For interpolation, the inverse-distance algorithm in the Tecplot-36066 software is used with default parameters (exponent = 3.5, point selection = Octant, Number of points = 8). The overlaid white solid and dashed lines show the upper limits and the weighted-averages with ±1 standard deviation of the low-frequency broadbands for Mach 2–10, respectively. (b) NCE of the PSD for Mach 2–10 used to define the upper limit of the low-frequency broadband.

Close modal

The characteristic length, Ls, defined by the density-gradient shock-thickness is calculated as the overall density change divided by the maximum density gradient (e.g., see Ref. 72, Chap. X, Sec. 9) and is shown in Fig. 12(a). We note that λ1, the upstream mean-free-path, is obtained through the VHS model with a viscosity index of ω=0.81. The SUGAR-1D DSMC shock-thickness values match within 2% with the DSMC calculations of Macrossan and Lilley42 and Bird (see Ref. 7, p. 294) for the same viscosity index. We note that their results are scaled by 76% to account for the differences in λ1, which they defined based on the hard-sphere model.7 The noticeable discrepancy between the SUGAR-1D DSMC results for ω=0.81 with the experiments of Alsmeyer40 is due to the choice of viscosity index, while qualitatively the variation of shock-thickness with Mach number is consistent with the experimental results. Good agreement is observed between SUGAR-1D DSMC and experiment for a Mach 8 shock simulated with ω=0.75. Based on these values, the calculated ratio of shock-thickness to upstream bulk velocity, Lsux,11, decreases with Mach number, as shown in Fig. 12(b).

FIG. 12.

(a) Reciprocal of the density-gradient shock-thickness normalized by the upstream mean-free-path as a function of Mach number. The SUGAR 1D DSMC results use a viscosity index of ω=0.81 (open symbols) and ω=0.75 (filled circular symbol) at Mach 8. (b) Lsux,11 and Strouhal number, St, as a function of Mach number, including the St at Mach 8, ω=0.75 (filled circular symbol). The standard deviation in the Strouhal number is based on the standard deviation in weighted average frequency shown in Fig. 11(a).

FIG. 12.

(a) Reciprocal of the density-gradient shock-thickness normalized by the upstream mean-free-path as a function of Mach number. The SUGAR 1D DSMC results use a viscosity index of ω=0.81 (open symbols) and ω=0.75 (filled circular symbol) at Mach 8. (b) Lsux,11 and Strouhal number, St, as a function of Mach number, including the St at Mach 8, ω=0.75 (filled circular symbol). The standard deviation in the Strouhal number is based on the standard deviation in weighted average frequency shown in Fig. 11(a).

Close modal

The Strouhal number defined based on the above quantities and the weighted average frequency, f, shown in Fig. 11(a), is relatively constant within a range of St =0.007–0.011 and a standard deviation from 0.001 to 0.02. For Mach 8 case simulated with ω=0.75, St was found to be in the same range and equal to 0.0092. This range of Strouhal numbers is similar to those corresponding to low-frequency unsteadiness of shock and separation bubble, as observed in the literature for highly compressible flows over embedded bodies, indicating that the two phenomena could be related. For example, in the case of shock-dominated separated flows the Strouhal number associated with low-frequency shock motion ranges from 0.02 to 0.05 (e.g., Refs. 2, 23, 53, and 73–75), assuming that the characteristic length and velocity scales are given by the length of the separation bubble and the upstream bulk velocity, respectively. In the study of oblique SBLIs, Nichols et al.76 defined the Strouhal number based on the upstream boundary-layer thickness and freestream velocity and found it to be within a range of 0.0003–0.05. Therefore, we hypothesize that these fluctuations may play a major role in shock-dominated flows having shock-thickness comparable to other important length scales in the flow, such as the size of the boundary layer in SBLIs and turbulent length scale in shock-turbulence interaction.

To test whether the established range of Strouhal numbers holds for variations in upstream temperature, another set of DSMC cases were simulated in which the upstream temperature was changed from 710 K by fractions of 1/8, 1/4, 1/2, and 2, while the upstream bulk velocity was varied accordingly to maintain a constant Mach number of 7.2. Figure 13 shows a decrease in the ratio of Lsux,11 with increase in temperature, due primarily to the increase in ux,1 rather than the shock-thickness. The Strouhal numbers obtained using the aforementioned approach are within the range of St = 0.005–0.011, with a small reduction with decrease in temperature. Nonetheless, the ±1 standard deviation of broadband frequencies is contained within the limits of St =0.001 to 0.02.  Appendix A also shows that the Strouhal number remained in the same range for a case run by increasing the freestream number density, n1, by three orders of magnitude, where the mean-collision-time of molecules decreased, frequency increased, and shock-thickness decreased by the same order. In Ref. 68, we show that a semi-empirical relationship to predict low-frequency oscillations in shocks can be derived from bimodal probability distribution functions based on noncentral Chi-squared functions for a broad range of Mach numbers and freestream temperatures.

FIG. 13.

Strouhal numbers for a range of cases simulated at Mach 7.2, as a function of freestream temperature Ttr,1. A viscosity index of ω=0.81 was used in the DSMC simulations.

FIG. 13.

Strouhal numbers for a range of cases simulated at Mach 7.2, as a function of freestream temperature Ttr,1. A viscosity index of ω=0.81 was used in the DSMC simulations.

Close modal

Finally, our Strouhal number formulation was used to qualitatively explain the ∼300 kHz frequencies observed inside a shock in the experiments of Chynoweth et al. (see Ref. 77, Fig. 7) of a Mach 6 air flow and Re1=6.73×106 m−1 over a 7° half-angle cone attached with a 10° fin. At freestream conditions experienced by the fin downstream of the leading-edge shock, that is, ρ1=4.636×102 kg m−3, Mach 5, and freestream temperature, Ttr,1=71.67 K, the normal shock-thickness, Ls 6.6 μm, was estimated from Fig. 12(a) and λ1 calculated using the VHS molecular model. Using Eq. (14) and estimating St from Fig. 12(b), the low-frequency broadband was predicted to have a weighted average frequency of 1.2 MHz with one standard deviation limit from 0.4 to 2.1 MHz, which was also verified from a separate DSMC simulation. The broadband from the simulation also revealed high spectral energy near 50 kHz frequency, which is also present in the experimental PSD spectrum. This estimate is qualitatively reasonable, even though it ignores factors that could reduce the range of dominant frequencies, such as differences in gas properties, reduced strength of the oblique shock generated by the fin compared to the normal shock, and surface effects.

The investigation of macroscopic fluctuations in a DSMC-computed Mach 7.2 shock layer revealed two orders of magnitude lower dominant frequencies of fluctuations than found in the freestream. In the peak overshoot region in the shock, a low-frequency broadband with a weighted-average of ∼37 kHz was observed as opposed to the freestream region where peak frequencies are concentrated at ∼3.5 MHz corresponding to the inverse mean-collision-time. These disparities were attributed to the differences in particle distribution functions in the nonequilibrium zone of shock vs the equilibrium regime upstream. The fluctuations in the x-directional normalized overall stress component, which is the mean of the PDF fξx, were found to be correlated with perturbations in the entire X-directional energy space of PDF. Two distinct energy bins were identified depending on the sign of the cross-correlation coefficient, and a Lotka–Volterra type two-energy bin ODE model was constructed.

The model accounted for the interaction of two energy bins through intermolecular collisions and particle fluxes from neighboring computational cells. At a location inside the shock, the model predicted two disparate timescales: a longer timescale (50–100 kHz) associated with the time of decay of fluctuations in the number of particles in energy bins to critical values and an order of magnitude smaller timescale (281 kHz) associated with oscillations in the number of particles. In the freestream, the model predicted faster decay of fluctuations, O(MHz), and more importantly, an absence of low-frequency fluctuations consistent with the DSMC spectral analysis.

Finally, a Strouhal number was defined based on the density gradient shock-thickness and upstream bulk velocity to nondimensionalize the low-frequency broadband characterized by a weighted average frequency and ±1 standard deviation across a wide range of Mach numbers from 2 to 10. The Strouhal number was found to range from 0.007 to 0.011 with ±1 standard deviation between 0.001 and 0.02, consistent with values found in the literature corresponding to low-frequency unsteadiness of shock and separation bubbles in shock-dominated flows. The range of Strouhal numbers was also found to hold for a Mach 7.2 shock simulated with different upstream temperatures and number density, and explain qualitatively the frequencies inside a shock observed in the work of Chynoweth et al.77 

The presence of low-frequency fluctuations of shock layers suggests that these disturbances may play a key role in the receptivity process of transition in hypersonic flows, SBLIs, and shock–turbulence interactions, especially when the shock-thickness is comparable to other important length scales in the flow. Our work implies that the leading-edge shock interactions with freestream disturbances, such as tunnel-induced turbulence, can result in selective amplification of disturbances in the low-frequency range, and consequent excitation of the downstream boundary layer. Future numerical investigations focused on understanding the effect of freestream disturbances on transition should not ignore shock-induced perturbations in the boundary layer within the range of Strouhal numbers established in this work.

The research conducted in this paper is supported by the Office of Naval Research under Grant No. N000141202195 titled “Multi-scale modelling of unsteady shock-boundary layer hypersonic flow instabilities,” with Dr. Eric Marineau as the program officer.

This work used the STAMPEDE2 supercomputing resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE) at the Texas Advanced Computing Center (TACC) through Allocation No. TG-PHY160006. In addition, S.S.S. would like to thank colleague Nakul Nuwal for helpful discussions regarding the two-energy bin model.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

First, we discuss the frequencies associated with potential disturbances introduced due to the STABIL algorithm for the Mach 7.2 shock case described in Sec. II B. The STABIL algorithm works as follows: when the number of particles in the domain boundary increases by a small, pre-set number, then all particles are translated by a small distance, δx, in the direction opposite to the drift of the shock. After translation, molecules that move outside of the boundary are deleted and those that are within a δx distance from the opposite boundary are replicated. One-way spurious numerical frequencies could be introduced at the boundaries is due to the recurring frequency at which additional particles are added within a δx distance from the boundaries to stabilize the shock. Starting with the upstream side, it was found that particles were not added near the upstream boundary during stabilization. On the downstream side, Fig. 14 shows the number of particles added during stabilization as a function of time and time passed since the last instance of particle addition. The mean values of these quantities are 236 (or ∼6.2% of the downstream computational particles in a root cell) and 0.0018 ms, respectively. We note that mean time between adding particles corresponds to a frequency of ∼553 kHz, which is much higher than the frequency of low-frequency fluctuations that are discussed in Sec. II. It is possible that the disturbances generated due to these particle additions could propagate to the shock at the downstream speed of sound, but the frequency associated with such disturbances would be 205 kHz, which is again outside of the low-frequency broadband. Therefore, this analysis supports the conclusion that molecular frequencies at the peak over shoot of the shock are not numerical artifacts of the STABIL algorithm.

FIG. 14.

For stabilization of the Mach 7.2 shock, additional particles added in the vicinity of the downstream boundary and the time passed since the last instance of particle addition.

FIG. 14.

For stabilization of the Mach 7.2 shock, additional particles added in the vicinity of the downstream boundary and the time passed since the last instance of particle addition.

Close modal

Next, we analyze the possible effect of numerical amplification of fluctuations on our low-frequency predictions caused by the stochastic nature of DSMC for Fn>1. The best way to assess this is to study its effect on the PSD, as shown in Fig. 15. The DSMC simulation, discussed in Sec. II B, was re-run first with Fn reduced by four orders of magnitude, which now requires 33 million particles, 1.2 × 106 time steps with Δt = 1 ns, for the same number density of 1022 m−3. An additional DSMC simulation was also performed for the case of Fn = 1 (i.e., no numerical amplification of fluctuation amplitudes) with the freestream number density now raised by three orders of magnitude, n1=1025 m−3. We note that we have changed the abscissa from frequency to Strouhal number using Eq. (14) because when we change number density we do not expect the “line-by-line” PSD spectrum to remain the same (due to the molecular nature of the flow), but we do expect the mean Strouhal numbers to be similar if we are looking at a phenomena related to fluctuations induced by the motion of the shock and not computational noise. A comparison of the new PSD spectra in the shock with the original case for the 1022 m−3 upstream condition is shown in Fig. 15(a) where indeed we find only a change of 2.7% in the weighted average frequency of the low-frequency broadband in the shock such that the mean Strouhal number is St0.012. Finally, Fig. 15(b) shows the PSD obtained from the simulation run by increasing the freestream number density by three orders of magnitude, n1=1025 m−3 and Fn = 1. Comparison of this result with Fig. 15(a) shows that both conditions generate PSDs consistent with the Strouhal number scaling discussed in Sec. V and produce similar quantities such as mean Strouhal number and broadband limits.

FIG. 15.

Comparison of PSD of the X-directional normal stress, σxx, inside the shock near the location of absolute maximum velocity gradient, P, for Mach 7.2 argon shock. (a) n1= 1022 m−3 and Fn = 107 and 103 and (b) n1= 1025 m−3 and Fn = 1.

FIG. 15.

Comparison of PSD of the X-directional normal stress, σxx, inside the shock near the location of absolute maximum velocity gradient, P, for Mach 7.2 argon shock. (a) n1= 1022 m−3 and Fn = 107 and 103 and (b) n1= 1025 m−3 and Fn = 1.

Close modal

We first describe the grouping of detailed collision processes pi into like-Pi processes of Table II based on the energy density functions at locations P and F. Starting with the reduction of collision processes pi to like-Pi processes, Fig. 8 shows that the energy distributions of Ac and Ad are very similar at probes P and F, so that Ad-type particles can be considered to be the same as those of Ac. In addition, although the distributions fξxA and fξy+ξzA of Ac and Ad particles are the same at probe F, they are slightly different at probe P at high ξx and low ξy+ξz energies. This occurs because when Ac- and B-type particles collide, the B particle loses its X-directional energy, which results in the increase in not only the ξx of the Ac-type particle but also a noticeable increase in its ξy+ξz energy, as can be seen in the transverse energy distributions of Ac-type particles in Fig. 8(b) at probe P. Also, when two Ad particles collide, one of them is switched from energy bin A to bin B at the expense of higher transverse energy of the other Ad particle that remains in bin A. In the freestream, since these two Ad-type particles are the same Ac-type particles that were generated from process pc, their transverse energy distributions exactly match, as seen in Fig. 8(d). However, at probe P, it is also possible that only one of the two Ad particles is an Ac-type, and the other one has lower transverse and higher X-directional energy than the Ac-type particle, as can be seen in the slightly larger X-directional energy distribution of Ad-type particles than the Ac-type in Fig. 8(a) between 0.5<ξx<1.33. We demonstrate that small distinctions such as this do not affect the dynamics of the system and assume that AcAd to construct processes Pc and Pd given in Table II.

Similarly, we can describe the process pe, where the particles of type-Ae have a very different transverse energy distribution shown in Fig. 8(b) than the other type-A particles at probe P. This process does not occur at probe F. Ignoring the differences at low transverse energies, the distribution of Ae and Af particles is the same, as shown in Fig. 8(b), and therefore, they can be denoted by the same identifier Ae.

We move on to present the derivations for the simplifications to Eq. (8) that allow us to reduce the six collision processes (Pi) given in Table II to four processes (Qi). To construct compound collision processes, Qi we start with the rate coefficient ki for process Pi, which can be evaluated from the DSMC simulation as

ki=CitΔtNR1tNR2t,
(B1)

where the collision takes place with between particle types R1 and R2. Cit for collision processes defined in Table II are listed in Table IV. Using the definition of rates ka and kb, the first kinetic equation in Eq. (7) can be written as

[dNAadt]collCatΔtNAtNBtNANB+CbtΔtNBt2NB2,[dNAadt]collk̃aNANB+kbNB2,
(B2)

where we have used the fact that

NAaNAatNANAt,
(B3)

as shown in Fig. 16(a). The ratio NAa/NAat is estimated from Fig. 8(a) by obtaining the number of particles between 0.0<ξx<1.33 at probe P. At probe F, the ratio of NAa/NAt cannot be directly estimated from Fig. 8(c) because the distribution of Aa-type particles overlaps with other A-type particles. Yet, Eq. (B3) is assumed to hold true under the assumption that the entire X-directional energy zone 28<ξx<0 is coarsened into a single bin A.

FIG. 16.

(a) and (b) show, at probe P, the comparison of time variation of ratio NAa/NAat with NA/NAt and the time variation of ratio of collisions Cc/Cd, respectively, to justify assumptions in Eqs. (B3) and (B5), respectively. The data points in (a) and (b) are obtained by taking time averages over a time window of 0.3 and 3 μs, respectively, to reduce statistical scatter.

FIG. 16.

(a) and (b) show, at probe P, the comparison of time variation of ratio NAa/NAat with NA/NAt and the time variation of ratio of collisions Cc/Cd, respectively, to justify assumptions in Eqs. (B3) and (B5), respectively. The data points in (a) and (b) are obtained by taking time averages over a time window of 0.3 and 3 μs, respectively, to reduce statistical scatter.

Close modal

Using the definitions of rates kc and kd, the second kinetic equation in Eq. (7) can be written as

[dNAcdt]coll=CctΔtNActNBtNAcNB+CdtΔtNAct2NAc2.
(B4)

Since the average number of Ac particles is not known, it is difficult to estimate the second term on the right-hand side. However, this term can be simplified by observing from Fig. 16(b) that the ratio of rates Cc/Cd does not vary significantly from the average rates at probe P, that is,

CcCdCctCdt,
(B5)

and the standard deviation in the fluctuation of the ratio of Cc/Cd is only 3.7% of the average ratio of 3.63 for probe P and 4.2% of an average ratio of one for probe F. Using Eq. (B1) for Cc and Cd in B5, and the definition of kc based on the average rate, we obtain

kdNAc2CdtΔtNActNBtNAcNB.
(B6)

In addition, we can see from Figs. 8(a) and 8(b) that the Ac particles have the same energy distributions as the Ah at probe P, which implies that any particle that takes part in process Pc can take part in process Ph, and together, they constitute 72.37% of total collisions in which the particles in energy bin A take part. Therefore, we can expect the changes in the ratio NAc/NAct are correlated with the changes in ratio NA/NAt and write

[dNAcdt]collCctCdtΔtNAtNBtNANBk̃cNANB.
(B7)

The same assumption holds at probe F because the Ac-type particles constitute the majority of particles from 20<ξx<5, which is a significant portion of energy bin A (28<ξx<0), as seen from Fig. 8(c).

By using the aforementioned assumptions, the third kinetic equation in Eq. (7) becomes

[dNAcdt]coll=2kcNAcNB2kdNAc22k̃cNANB.
(B8)

We also drop the second term on the right-hand side of the fourth kinetic equation in Eq. (7) because CetCft. Note that both Cet=Cft=0 at probe F.

By substituting Eqs. (B2), (B7), and (B8), into Eq. (8), we obtain

[dNAdt]coll=k̃aNANB+kbNB2+k̃cNANB+2keNB2,
(B9)

which when substituted into Eqs. (3) and (4) gives the final system of dynamic equations, Eq. (9).

1.
S.
Lee
,
S. K.
Lele
, and
P.
Moin
, “
Interaction of isotropic turbulence with shock waves: Effect of shock strength
,”
J. Fluid Mech.
340
,
225
(
1997
).
2.
S.
Priebe
,
J. H.
Tu
,
C. W.
Rowley
, and
M. P.
Martín
, “
Low-frequency dynamics in a shock-induced separated flow
,”
J. Fluid Mech.
807
,
441
477
(
2016
).
3.
X.
Zhong
, “
High-order finite-difference schemes for numerical simulation of hypersonic boundary-layer transition
,”
J. Comput. Phys.
144
,
662
(
1998
).
4.
J.
Sesterhenn
, “
A characteristic-type formulation of the Navier–Stokes equations for high order upwind schemes
,”
Comput. Fluids
30
,
37
(
2000
).
5.
H. W.
Liepmann
,
R.
Narasimha
, and
M. T.
Chahine
, “
Structure of a plane shock layer
,”
Phys. Fluids
5
,
1313
(
1962
).
6.
G.
Bird
, “
Aspects of the structure of strong shock waves
,”
Phys. Fluids
13
,
1172
(
1970
).
7.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
, 2nd ed. (
Clarendon Press
,
1994
), p.
484
.
8.
H. S.
Ribner
, “
Convection of a pattern of vorticity through a shock wave
,”
Technical Report No. NACA-TR-1164
(
National Advisory Committee for Aeronautics
,
1954
).
9.
H. S.
Ribner
, “
Shock-turbulence interaction and the generation of noise
,”
Technical Report No. NACA-TR-1233
(
National Advisory Committee for Aeronautics
,
1954
).
10.
F. K.
Moore
, “
Unsteady oblique interaction of a shock wave with a plane disturbance
,”
Technical Report No. NACA TN 2879
(
National Advisory Committee for Aeronautics
,
1953
).
11.
L. S. G.
Kovasznay
, “
Turbulence in supersonic flow
,”
J. Aeronaut. Sci.
20
,
657
(
1953
).
12.
C.-T.
Chang
, “
Interaction of a plane shock and oblique plane disturbances with special reference to entropy waves
,”
J. Aeronaut. Sci.
24
,
675
(
1957
).
13.
M. V.
Morkovin
,
Mécanique de la Turbulence
(
Centre National de la Recherche Scientifique
,
Paris
,
1962
), p.
367
.
14.
K.
Mahesh
and
S.
Lee
, “
The interaction of an isotropic field of acoustic waves with a shock wave
,”
J. Fluid Mech.
300
,
383
(
1995
).
15.
K.
Mahesh
,
S. K.
Lele
, and
P.
Moin
, “
The influence of entropy fluctuations on the interaction of turbulence with a shock wave
,”
J. Fluid Mech.
334
,
353
(
1997
).
16.
Y.
Andreopoulos
,
J. H.
Agui
, and
G.
Briassulis
, “
Shock wave—turbulence interactions
,”
Annu. Rev. Fluid Mech.
32
,
309
(
2000
).
17.
J.
Larsson
and
S. K.
Lele
, “
Direct numerical simulation of canonical shock/turbulence interaction
,”
Phys. Fluids
21
,
126101
(
2009
).
18.
K.
Koffi
,
Y.
Andreopoulos
, and
C. B.
Watkins
, “
Dynamics of microscale shock/vortex interaction
,”
Phys. Fluids
20
,
126102
(
2008
).
19.
H.
Xiao
and
R.
Myong
, “
Computational simulations of microscale shock–vortex interaction using a mixed discontinuous Galerkin method
,”
Comput. Fluids
105
,
179
(
2014
).
20.
S.
Singh
,
A.
Karchani
, and
R.
Myong
, “
Non-equilibrium effects of diatomic and polyatomic gases on the shock-vortex interaction based on the second-order constitutive model of the Boltzmann-Curtiss equation
,”
Phys. Fluids
30
,
016109
(
2018
).
21.
D. S.
Dolling
, “
Fifty years of shock-wave/boundary-layer interaction research: What next?
,”
AIAA J.
39
,
1517
(
2001
).
22.
H.
Babinsky
and
J. K.
Harvey
,
Shock Wave–Boundary-Layer Interactions
(
Cambridge University Press
,
2011
).
23.
D. V.
Gaitonde
, “
Progress in shock wave/boundary layer interactions
,”
Prog. Aerosp. Sci.
72
,
80
(
2015
).
24.
A. V.
Fedorov
, “
Receptivity of a high-speed boundary layer to acoustic disturbances
,”
J. Fluid Mech.
491
,
101
129
(
2003
).
25.
Y.
Ma
and
X.
Zhong
, “
Receptivity of a supersonic boundary layer over a flat plate. Part 1. Wave structures and interactions
,”
J. Fluid Mech.
488
,
31
(
2003
).
26.
Y.
Ma
and
X.
Zhong
, “
Receptivity of a supersonic boundary layer over a flat plate. Part 2. Receptivity to free-stream sound
,”
J. Fluid Mech.
488
,
79
(
2003
).
27.
Y.
Ma
and
X.
Zhong
, “
Receptivity of a supersonic boundary layer over a flat plate. Part 3. Effects of different types of free-stream disturbances
,”
J. Fluid Mech.
532
,
63
(
2005
).
28.
C.
Hader
and
H. F.
Fasel
, “
Towards simulating natural transition in hypersonic boundary layers via random inflow disturbances
,”
J. Fluid Mech.
847
,
R3
(
2018
).
29.
L. M.
Mack
,
Boundary-Layer Linear Stability Theory
, AGARD Report No. 709 (
Jet Propulsion Laboratory, California Institute of Technology
,
1984
).
30.
N.
Zavol'skii
and
V.
Reutov
, “
Thermal excitation of waves in a boundary layer
,”
Fluid Dyn.
18
,
701
(
1984
).
31.
A. V.
Fedorov
and
A. P.
Khokhlov
, “
Prehistory of instability in a hypersonic boundary layer
,”
Theor. Comput. Fluid Dyn.
14
,
359
(
2001
).
32.
A.
Fedorov
,
A.
Shiplyuk
,
A.
Maslov
,
E.
Burov
, and
N.
Malmuth
, “
Stabilization of a hypersonic boundary layer using an ultrasonically absorptive coating
,”
J. Fluid Mech.
479
,
99
(
2003
).
33.
A. V.
Fedorov
, “
Receptivity of a supersonic boundary layer to solid particulates
,”
J. Fluid Mech.
737
,
105
(
2013
).
34.
L. D.
Landau
and
E. M.
Lifshitz
,
Statistical Physics: Part 1
, 3rd ed. (
Pergamon Press
,
1980
), Vol.
5
.
35.
P.
Luchini
, in
Seventh IUTAM Symposium on Laminar-Turbulent Transition
(
Springer
,
2010
), pp.
11
18
.
36.
P.
Luchini
, “
Receptivity to thermal noise of the boundary layer over a swept wing
,”
AIAA J.
55
,
121
(
2017
).
37.
A.
Fedorov
and
A.
Tumin
, “
Receptivity of high-speed boundary layers to kinetic fluctuations
,”
AIAA J.
55
,
2335
(
2017
).
38.
L. D.
Edwards
and
A.
Tumin
, “
Model of distributed receptivity to kinetic fluctuations in high-speed boundary layers
,”
AIAA J.
57
,
4750
(
2019
).
39.
S. K.
Stefanov
,
I. D.
Boyd
, and
C.-P.
Cai
, “
Monte Carlo analysis of macroscopic fluctuations in a rarefied hypersonic flow around a cylinder
,”
Phys. Fluids
12
,
1226
(
2000
).
40.
H.
Alsmeyer
, “
Density profiles in argon and nitrogen shock waves measured by the absorption of an electron beam
,”
J. Fluid Mech.
74
,
497
(
1976
).
41.
C.
Cercignani
,
A.
Frezzotti
, and
P.
Grosfils
, “
The structure of an infinitely strong shock wave
,”
Phys. Fluids
11
,
2757
(
1999
).
42.
M. N.
Macrossan
and
C. R.
Lilley
, “
Viscosity of argon at temperatures >2000 K from measured shock thickness
,”
Phys. Fluids
15
,
3452
(
2003
).
43.
T.
Ozawa
,
D. A.
Levin
,
I.
Nompelis
,
M.
Barnhardt
, and
G. V.
Candler
, “
Particle and continuum method comparison of a high-altitude, extreme-Mach-number reentry flow
,”
J. Thermophys. Heat Transfer
24
,
225
(
2010
).
44.
T. E.
Schwartzentruber
and
I. D.
Boyd
, “
A hybrid particle-continuum method applied to shock waves
,”
J. Comput. Phys.
215
,
402
(
2006
).
45.
T.
Zhu
,
Z.
Li
, and
D. A.
Levin
, “
Modeling of unsteady shock tube flows using direct simulation Monte Carlo
,”
J. Thermophys. Heat Transfer
28
,
623
(
2014
).
46.
G.
Bird
, “
Recent advances and current challenges for DSMC
,”
Comput. Math. Appl.
35
,
1–14
(
1998
).
47.
S.
Stefanov
,
V.
Roussinov
, and
C.
Cercignani
, “
Rayleigh–Bénard flow of a rarefied gas and its attractors. I. Convection regime
,”
Phys. Fluids
14
,
2255
(
2002
).
48.
S.
Stefanov
,
V.
Roussinov
, and
C.
Cercignani
, “
Rayleigh–Bénard flow of a rarefied gas and its attractors. II. Chaotic and periodic convective regimes
,”
Phys. Fluids
14
,
2270
(
2002
).
49.
S.
Stefanov
,
V.
Roussinov
, and
C.
Cercignani
, “
Rayleigh–Bénard flow of a rarefied gas and its attractors. III. Three-dimensional computer simulations
,”
Phys. Fluids
19
,
124101
(
2007
).
50.
M. A.
Gallis
,
T. P.
Koehler
,
J. R.
Torczynski
, and
S. J.
Plimpton
, “
Direct simulation Monte Carlo investigation of the Rayleigh-Taylor instability
,”
Phys. Rev. Fluids
1
,
043403
(
2016
).
51.
M. A.
Gallis
,
T. P.
Koehler
,
J. R.
Torczynski
, and
S. J.
Plimpton
, “
Direct simulation Monte Carlo investigation of the Richtmyer-Meshkov instability
,”
Phys. Fluids
27
,
084105
(
2015
).
52.
O.
Tumuklu
,
D. A.
Levin
, and
V.
Theofilis
, “
Investigation of unsteady, hypersonic, laminar separated flows over a double cone geometry using a kinetic approach
,”
Phys. Fluids
30
,
046103
(
2018
).
53.
O.
Tumuklu
,
V.
Theofilis
, and
D. A.
Levin
, “
On the unsteadiness of shock–laminar boundary layer interactions of hypersonic flows over a double cone
,”
Phys. Fluids
30
,
106111
(
2018
).
54.
S. S.
Sawant
,
O.
Tumuklu
,
R.
Jambunathan
, and
D. A.
Levin
, “
Application of adaptively refined unstructured grids in DSMC to shock wave simulations
,”
Comput, Fluids
170
,
197
(
2018
).
55.
S. S.
Sawant
,
O.
Tumuklu
,
V.
Theofilis
, and
D. A.
Levin
, in
The IUTAM Transition 2019 Proceedings
(
Springer
,
2019
).
56.
S. S.
Sawant
,
V.
Theofilis
, and
D. A.
Levin
, “
Kinetic modelling of three-dimensional shock/laminar separation bubble instabilities in hypersonic flows over a double wedge
,” arXiv:2101.08957 [physics.flu-dyn] (
2021
).
57.
A. J.
Lotka
, “
Contribution to the theory of periodic reactions
,”
J. Phys. Chem.
14
,
271
(
1910
).
58.
A. J.
Lotka
, “
Analytical note on certain rhythmic relations in organic systems
,”
Proc. Natl. Acad. Sci. U. S. A.
6
,
410
(
1920
).
59.
V.
Volterra
, “
Fluctuations in the abundance of a species considered mathematically
,”
Nature
118
,
558
(
1926
).
60.
M.
Ivanov
and
S.
Rogasinsky
, “
Analysis of numerical techniques of the direct simulation Monte Carlo method in the rarefied gas dynamics
,”
Russ. J. Numer. Anal. Math. Model.
3
,
453
(
1988
).
61.
T.
Ohwada
, “
Structure of normal shock waves: Direct numerical analysis of the Boltzmann equation for hard-sphere molecules
,”
Phys. Fluids A
5
,
217
(
1993
).
62.
N. G.
Hadjiconstantinou
,
A. L.
Garcia
,
M. Z.
Bazant
, and
G.
He
, “
Statistical error in particle simulations of hydrodynamic phenomena
,”
J. Comput. Phys.
187
,
274
(
2003
).
63.
P.
Welch
, “
The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms
,”
IEEE Trans. Audio Electroacoust.
15
,
70
(
1967
).
64.
O. M.
Solomon
, Jr.
, “
PSD computations using Welch's method
,”
Technical Report No. SAND-91-1533 ON: DE92007419
(
Sandia National Laboratories
,
Albuquerque, NM
,
1991
).
66.
Tecplot-360
, see https://www.tecplot.com/products/tecplot-360/ for the Tecplot software documentation (2020 R1).
67.
M. N.
Kogan
,
Rarefied Gas Dynamics
, 1st ed. (
Springer
US
,
1969
), p.
515
.
68.
S. S.
Sawant
,
D. A.
Levin
, and
V.
Theofilis
, “
Analytical prediction of low-frequency fluctuations inside a one-dimensional shock
,”
Theor. Comput. Fluid Dyn.
(published online
2021
).
69.
H. M.
Mott-Smith
, “
The solution of the Boltzmann equation for a shock wave
,”
Phys. Rev.
82
,
885
(
1951
).
70.
The movies loop over time windows from −20 to 646, where w=−20 to 0 correspond to the transient time, which is not used in the calculation of cross-correlation coefficient and time-averaged means.
71.
J. D.
Anderson
,
Modern Compressible Flow: With Historical Perspective
, 3rd ed. (
Tata McGraw-Hill
,
2003
).
72.
W. G.
Vincenti
and
C. H.
Kruger
,
Introduction to Physical Gas Dynamics
(
Wiley
,
New York
,
1965
).
73.
J.-P.
Dussauge
,
P.
Dupont
, and
J.-F.
Debiève
, “
Unsteadiness in shock wave boundary layer interactions with separation
,”
Aerosp. Sci. Technol.
10
,
85
(
2006
).
74.
S.
Piponniau
,
J.-P.
Dussauge
,
J.-F.
Debieve
, and
P.
Dupont
, “
A simple model for low-frequency unsteadiness in shock-induced separation
,”
J. Fluid Mech.
629
,
87
(
2009
).
75.
N. T.
Clemens
and
V.
Narayanaswamy
, “
Low-frequency unsteadiness of shock wave/turbulent boundary layer interactions
,”
Annu. Rev. Fluid Mech.
46
,
469
(
2014
).
76.
J. W.
Nichols
,
J.
Larsson
,
M.
Bernardini
, and
S.
Pirozzoli
, “
Stability and modal analysis of shock/boundary layer interactions
,”
Theor. Comput. Fluid Dyn.
31
,
33
(
2017
).
77.
B. C.
Chynoweth
,
C. A. C.
Ward
,
R. O.
Henderson
,
C. G.
Moraru
,
R. T.
Greenwood
,
A. D.
Abney
, and
S. P.
Schneider
, AIAA Paper No. 2014-0074,
2014
.