Kinetic modeling of unsteady hypersonic flows over a tick geometry Free

In the spirit of this memorial issue to honor the contributions of the late Prof. Graeme A. Bird, we present our investigations into hypersonic flow separation and shock-wave boundary layer interactions (SWBLIs) over a “tick” geometry, a unique application of the direct simulation Monte Carlo (DSMC) method to unsteady, continuum flows. Nearly 15 years ago, Moss and Bird¹ applied the DSMC method to statistically model SWBLIs of the flow around the hollow cylinder flare and double cone. By comparing their numerical results with the experiment of Ref. 2, they showed that the DSMC method is capable of predicting separated flow properties in the near-continuum regime accurately, given that the numerical parameter requirements, such as time step, cell size, and the total number of particles, are satisfied. At that time, the computational power was insufficient to model higher Reynolds numbers or time-accurate simulations, which DSMC is uniquely well suited for. Since then, Tumuklu et al.³ performed time-accurate DSMC simulations on the same double-cone geometry combined with linear stability analysis to show that the required time to reach the steady solution increases with the Reynolds number. Further analyses on the origin of the unsteady mechanics were revealed in Ref. 4.

The DSMC method,⁵ the well-known stochastic approach to solve the Boltzmann equation of transport, is applicable to a wide Knudsen number spectrum with the increase in computer speed and massively parallel computing facilities. Previously, Bird⁶ performed hypersonic three-dimensional (3D) reentry calculations of the space shuttle orbiter for altitudes from 170 to 90 km by taking translational, rotational, vibrational, and chemical nonequilibrium effects into account. Since then, computationally much more demanding unsteady DSMC calculations were conducted to simulate micro-Rayleigh-Bénard convection⁷ with the differences between two-dimensional (2D) and 3D configurations in the flow structures reported in Ref. 8. Similarly, Gallis et al.⁹ simulated the Taylor-Green vortex flow using DSMC and the turbulent kinetic energy dissipation rate was found to be in perfect agreement with the direct numerical simulation of Navier-Stokes (NS). Again, this important finding suggests that the DSMC method can also be used to study turbulent flows (i.e., in the continuum regime) since the DSMC simulations inherently account for the thermal fluctuations which can lead to turbulence.

In this paper, we consider hypersonic flow separation and laminar SWBLIs which have received considerable attention.^10,11 These interactions can lead to the laminar and turbulent transition, unsteadiness, and localized high pressure and heating regions that can have detrimental effects on the performance of hypersonic vehicles. Therefore, accurate predictions of these phenomena play a crucial role for design purposes. In this regard, many experimental,^12,13 theoretical, and numerical works^14,15 have been conducted to investigate the origin of SWBLIs. Chapman et al.¹⁰ conducted experimental and theoretical studies to investigate flow separation on bases, compression corners for Mach numbers changing from 0.4 to 3.6 with the Reynolds number between 4000 and 5 000 000. The pressure distribution was found to be strongly dependent on the location of transition relative to the reattachment and separation points. Experimental studies on separated hypersonic flows on double wedge and cone, hollow cylinder/flare, and biconic configurations^16,17 were conducted. However, the time characterization of such unsteady flows can be challenging in ground-based experiments primarily because of limited measurement times as well as the lack of a well-established metric to relate the establishment time of short-duration facilities with the time convergence of computations for separated flows.¹⁸ 

While many experimental^16,17 and numerical studies^19–21 on double cone, hollow cylinder, and wedge geometries have been conducted over decades, the leading edge separation which eliminates complicated boundary-layer development¹⁰ is the focus of recent studies. Experiments are being performed in the T-ADFA free-piston shock tunnel at the University of New South Wales (UNSW) Canberra.^22–24 This configuration, introduced first by Chapman,¹⁰ produces a large separation zone for high Mach and Reynolds numbers, commonly observed in hypersonic flight conditions. It provides an opportunity to test both continuum, NS and, once again, kinetic, DSMC computational approaches, especially for high stagnation enthalpies. Moss et al.²² showed that the impact of the slip velocity has a significant impact on the size of the separation region by comparing NS and DSMC solutions for a “tick” configuration in the near-continuum regime. In particular, the separation size predicted by NS was approximately 30 times larger than DSMC due to the zero slip assumption for a case with conditions similar to that of A2 that will be discussed in this work. Recently, Khraibut et al.²⁵ performed Computational Fluid Dynamics (CFD) calculations on the “tick” configuration at a steady state and found that they correlated with hypersonic triple deck theory. It was also observed that the wall temperature had a strong effect on the size of the separation region and for simulations run with wall temperatures between the free stream temperature of 165 K and the adiabatic wall temperature of 2770 K, the thickness of the boundary layer increased causing a larger separation region, consistent with the experimental observation of Ref. 26. Similarly, Prakash et al.²⁷ showed that the temporal evolution of the corner vortex of the tick geometry depends on the wall temperature by performing DSMC calculations using the SPARTA code.²⁸ 

In this work, we use the same time-accurate properties of the DSMC method to model hypersonic, expanding flows over the tick geometry for two conditions investigated at UNSW. The two conditions, known as “A2” and “E,” span a Knudsen number range from transitional to continuum, a Reynolds number range that increases only by a factor of four, and a free stream enthalpy for the first condition that is sufficiently high to cause chemical reactions in the working gas, air, but is too low to be significant for the continuum case, E. Again, without the efficient, massively parallel DSMC algorithm and present improved computational power, a time accurate study of the latter case would not have been possible. By investigating these two flow cases with different simulation setup conditions with a single gas dynamic formalism, we are able to evaluate the effect of tick geometry in terms of sharp vs blunted leading edges, the presence of multiple air species, chemical reactions, and the unsteady characteristics of the flows and surface fluxes.

Using the time-accurate DSMC solutions with the residual algorithm (RA) of Theofilis,²⁹ we present the general time characteristics of the tick, leading edge separation flows and are able to assess the experimental efforts in terms of measurement times vs time required for these range of flows to reach the steady state for the first time. Earlier modeling of Moss et al.²² using DSMC for steady state simulations only revealed large differences in terms of the size of the separation region due to rarefaction effects for a case with conditions similar to A2 and geometries with different bevel radii. In this work, we investigate for the first time the effect of leading edge bluntness, wall temperature, and surface accommodation coefficient on the time characteristics of the flowfield parameters and the evolution of the size and structure of the separation region using the DSMC/RA approach to model case A2 and the higher density and Reynolds number case E. The application of our coupled DSMC/RA technique to the higher Reynolds number case E also allows us to present a unique comparison of our velocity-slip, temperature-jump, and heat fluxes with those obtained from the CFD simulations of Khraibut et al.²⁵ Furthermore, we show that despite differences between the DSMC and CFD models of velocity-slip and temperature jump surface values, both gas dynamic approaches essentially predict similar heat flux profiles for condition E when the flow is fully developed. However, the time-accurate DSMC at a time earlier than the steady state (i.e., the CFD result) gives a bit better agreement with the experiment. The use of the RA analysis demonstrates that the measurement time for condition E is too short for the tick flow to have reached a steady state; however, it is sufficient for case A2.

The outline of the remainder paper is as follows. The collisional models employed in the DSMC simulations, the selection of numerical parameters, particularly for case E which is a continuum simulation, and the different tick geometries are described in Sec. II. To understand the time-dependent and steady state flow and surface parameters, the general flow features such as the expansion fan, leading edge shock, shear layer, triple point, recirculation region, and reattachment shock for cases A2 and E are identified and their Mach number contours are compared. The variation in local rarefaction is examined in the spatial distribution of the local Knudsen number for two geometries for case A2 (a sharp and beveled leading edge) and for the sharp leading edge for case E. In Sec. III, the effects of leading-edge bluntness and chemical reactions are presented for case A2. It can be seen that the dissociation of molecular oxygen and the diffusion of atomic oxygen affect the spatial distribution of species mole fractions and translational temperatures. These effects are less important for case E because the free stream energy is lower and the flow has a higher collision frequency. The effect of geometry, however, is important for case E, as discussed in this section. Section IV discusses the sensitivity of surface parameters to wall temperature, incomplete vibrational accommodation coefficient, and chemistry for case A2. For case E, a comparison is presented with the CFD simulations, as mentioned earlier. Finally, in Sec. V, the time evolution of the separation zone and its effect on surface parameters are discussed for both cases A2 and E. Using the RA, decay rates are compared for A2 and E and the regions of flow time activity are identified. It is found for the first time in a hypersonic expanding flow that a centrifugal instability exists in the region of the laminar separation bubble for case E.

In this work, particle-based numerical analyses were performed using the DSMC method,⁵ which models the Boltzmann equation of transport, provided that the simulation numerical parameters are properly selected. The Boltzmann equation of transport is given as

where f(c_i) and n are the distribution function of class c_i and number density, respectively.³⁰ Note that the external force field acting on gas molecules, F_i, is assumed to be zero in the current work. DSMC is well suited for resolving the transient characteristics of unsteady flows toward transition due to its time-accurate nature. It also provides the highest fidelity of molecular thermal nonequilibrium phenomena, an accurate model for viscosity, thermal conductivity, and diffusivity.³¹ The DSMC method inherently captures velocity and temperature slip with the Maxwell gas-surface interaction model without any further dependence on an a priori slip model. The numerical analyses were performed using the Statistical Modeling In Low-density Environment (SMILE)³² solver employing the majorant frequency scheme for modeling the molecular collision frequency³³ and the variable hard sphere (VHS) model for modeling the interaction between particles.⁵ The VHS cross section can be modeled as

where k is the Boltzmann constant, g is the relative speed, m_r is the reduced mass, and Γ is the gamma function. The running gas is air, and the corresponding VHS parameters of collision diameters, d_ref, and viscosity indices, ω, of each species at a reference temperature, T_ref, of 1000 °K are given in Table I.

In terms of collision models, the Larsen-Borgnakke (LB)³⁵ model for rotational-translational (R-T) energy transfer, with the rate proposed by Parker,³⁶ and vibrational-translational (V-T) energy transfer, with the rate proposed by Millikan and White,³⁷ were used (see Ref. 38 for details of the rotational and vibrational relaxation numbers). The correction factors of Lumpkin et al.³⁹ and Gimelshein et al.⁴⁰ were employed to convert the continuous (R-T) and (V-T) relaxation rates to relaxation probabilities. Furthermore, the gas-surface interactions were modeled using the Maxwell model with full and partially complete energy accommodation, as will be discussed in Sec. III. Finally, the total collision energy (TCE) approach⁵ was employed to model chemical reactions in the current work. Recently, Bird⁴¹ developed the quantum-kinetic (Q-K) theory to provide a more accurate chemistry model for highly nonequilibrium flows. However, since we expect small deviations from equilibrium distribution functions due to the relatively high-densities and low stagnation enthalpies in the tick-flow, the TCE approach was used. Furthermore, due to relatively low freestream stagnation enthalpies, as presented in Table III, the number of dissociation events of nitrogen will be negligibly small and therefore this reaction was also not modeled in the present work. However, the lower dissociation energy of molecular oxygen suggests that this chemical reaction will affect the flowfield properties and shock structure significantly. Therefore, in the current study, the exchange reactions of O atoms with N₂ molecules and O₂ dissociations through collisions with nitrogen and oxygen molecules, represented with M in the following chemical reactions

are modeled in order to test the effects of these chemical reactions using the Arrhenius rate parameters given in Table II. It should be noted that the recombination reaction rates^5,42 are significantly lower than the dissociation and exchange reaction rates under the current flow conditions. Therefore, recombination reactions involving atomic oxygen and nitrogen species are neglected in the current work.

Figure 1 shows the sharp “tick” geometry and three bevel radii used in the experiments and simulations. Due to the manufacturing imperfections and possible damage in the experiments, a bevel radius of 0.015 mm is considered first to determine the effect of bevel radius on the flow and surface parameters. Also, two larger radii of 0.1 and 0.2 mm bevels are studied to quantify the effect of the bevel size on the flow. It should be noted that the bevel radius, |R_i|, and the three different bevel configurations (i.e., i = 1, 2, and 3) are shown in the inset of Fig. 1. Also, the length of the expansion surface for different bevel sizes is shortened in reference to a sharp leading edge of 19.53 mm, as can be seen in the inset of Fig. 1.

Table III presents the two different freestream conditions studied in this work. The lower density case, A2, has a larger stagnation enthalpy, providing an opportunity to test the influence of chemical reactions on massively separated flows. The higher density case, condition E, has about a four times larger Reynolds number compared to A2 case, enabling us to study the time characteristics of massively separated flows. It should be noted that both these cases have regions of the continuum, transitional, and slip regimes which necessitates the use of a numerical approach applicable in such a wide Knudsen number spectrum. Additionally, the flow has been accelerated to high freestream velocities through the T-ADFA free-piston shock tunnel, resulting in a vibrationally frozen flow. Therefore, the freestream vibrational temperature of each species is different, as presented in Table III. In particular, at the exit of the nozzle in the reflected shock tunnel, the flow has a high level of thermal nonequilibrium and NO produced could be potentially used to measure the planar laser-induced fluorescence (PLIF) measurements of temperature and velocity field.

The general nature of the flow structures for cases A2 and E is shown in Fig. 2. The resulting flowfield for condition A2 shown in Fig. 2(a) has a diverse set of physical phenomena such as a Prandtl-Meyer expansion fan which enhances the thermal nonequilibrium, a leading edge shock emanating from the leading edge terminating at the triple point, and a shear layer emanating after the expansion and reattachment on the compression surface. The density for condition A2 drops by a factor of 64 near point 4 and increases by a factor 5 near location 16 with respect to the freestream value. The curved profile of the shear layer can be observed for condition A2 due to larger viscous interactions. Moreover, a small region of further expansion below the leading edge shock is presented for condition A2. Condition E, which is a higher density compared to condition A2, results in a more continuumlike, sharper shock structure with a larger separation region, as shown in Fig. 2(b). The higher Mach number causes a larger adverse pressure gradient, resulting in a larger separation region and a subsonic recirculation region, as shown in Fig. 2(c).

The sharp gradients in these flow conditions result in a range of Knudsen numbers at different locations. The local Knudsen number is a dimensionless parameter that determines the validity of the continuum formulation of fluid dynamics. In particular, if the Knudsen number exceeds the value of 0.05, the Navier-Stokes equations may not provide a correct representation of the flow.⁴³ The gradient based local Knudsen defined as^5,22

where λ and ρ are the local mean free path and density, respectively, is shown for the two conditions and a geometry with a bevel in Fig. 3. It can be seen from Fig. 3(a) that the Knudsen number increases up to a value of 0.8 in the vicinity of the leading edge which coincides with the region where the density gradients are the highest and shows that an appropriate numerical approach is needed to accurately simulate a flow at such high Knudsen numbers. Once the bevel is added to the geometry, the maximum Knudsen number observed in the flow decreases to a value of about 0.5 in the expansion surface and along the compression surface. On the other hand, when condition E is used, the Knudsen number decreases significantly and the flowfield can be considered to be in continuum except for the leading edge and the shear layer.

To resolve these sharp gradients that result from a sudden expansion of the flow through the leading edge and shock boundary layer interactions especially at the compression surface, the DSMC numerical parameters are carefully selected to ensure that the presented results have no dependency on the numerical parameters. Details of the correct numerical criteria for DSMC simulations were summarized in Ref. 3, and the numerical parameters in this work were selected accordingly. Table IV presents them for conditions A2 and E for the sharp leading edge geometry. Similar values were used for the bevel geometries.

Figure 4 shows local gas kinetic properties based on the numerical parameters presented in Table IV. The SMILE code uses a two-level fixed size Cartesian grid, namely, a coarse grid, or background grid, in which the sampling of the flowfield parameters takes place, and a refined grid constructed by subdividing background cells based on the local density gradient. In present work, the number of background cells is 680 by 350 in the X and Y directions, respectively. As can be seen in the inset of Fig. 4(b), these number of cells are sufficiently large to resolve macroparameter gradients especially at the leading edge. Note that the black rectangle at location A shows the same zoomed region in the inset. On the other hand, the size of collision cells is adapted based on the local mean free path and the maximum number of cell subdivisions of 60 × 60 is allowed for condition E. In particular, the number of subdivisions at locations A, B, and C are 121, 225, and 3600, respectively. Figures 4(a) and 4(b) show the ratio of the local mean free path to the collision cell size for conditions A2 and E. As shown, the cell size is sufficiently small such that the aforementioned ratio is larger than unity throughout the domain. In particular, for condition E, this criterion is met except for a small location of C with a value of about 0.8 where the shock interactions result in large density gradients. Nonetheless, our earlier simulations showed that this does not have an influence on results. It will be shown in Sec. III that the surface parameters are essentially the same for a degraded case for condition A2 with a cell size larger than the case given in Table IV. In fact, for the degraded case, the ratio of the local mean path to the size of collision cell becomes comparable with condition E shown in Fig. 4(b). Therefore, the size of collision cells is small enough to calculate the collision frequency accurately.

Since the number of collision cells required to obtain a resolution on the order of a mean free path is high, the corresponding number of simulated particles becomes very large especially for condition E where more than 1.8 × 10⁹ simulated particles were required. As shown in Figs. 4(c) and 4(d), the total number of particles per collision cell in the shock region was found to be more than 15. More specifically, the number of particles per species is about 11, 4, and 1 for N₂, O₂, and NO, respectively, for condition E, consistent with the mole fraction of each species given in Table III. Similar to condition E, the number of particles for each species especially per collision cell at location B is about 7, 15, 5, and 3 for O, N₂, O₂, and NO, respectively, for condition A2. Since the effects of chemical reactions are negligible for both cases, as will be discussed in Sec. III, and NO is a trace species, the low number of NO particles has no impact on the overall results. The number of particles for other species is sufficient since the majorant collision frequency algorithm is used and gives correct transport properties of each species in the selection of collision pairs in this work.³³ Finally, time steps of 2.5 × 10⁻⁹ and 2.0 × 10⁻⁹ s were chosen such that the ratio of the mean collision time to the time step is greater than unity, as can be seen in Figs. 4(e) and 4(f).

The effects of the leading edge bluntness on the size of the separation and the nonequilibrium are found to be significant for the lower density case. Figure 5 shows the spatial distribution of the structure of the separation zone and the translational temperatures at condition A2. It can clearly be seen that the size of the separation zone strongly depends on the bevel radius. In particular, the separation zone for the sharp leading edge case [Fig. 5(a)] is confined to a small region in the vicinity of the vertex and the separation point moves in the upstream direction and the reattachment point shifts toward the downstream direction with increasing bevel size. Due to the presence of the leading edge shock and the presence of the compression region due to the second surface, the translational temperatures in the downstream region of the leading edge are increased [i.e., between locations A and B of Fig. 5(b)]. In particular, when a larger bevel configuration is used, a localized bow shock occurs and changes the leading edge structure, resulting in an elevated temperature field in the vicinity of the leading edge, as can be seen in Figs. 5(c) and 5(d) (locations D and E). As shown in Fig. 2(a), the region between the leading edge shock and the shear layer experiences an expansion region where the temperatures slightly decrease at location C of Fig. 5(b). For the higher bevel size, the temperature values are found to decrease significantly in the expansion region due to interactions with the cold wall and the increase in the interaction time. The translational and vibrational temperatures for each species can be seen in Figs. 6 and 7, respectively. The translational temperature field shows only a slight difference between each species due to the high relaxation rate of the translation mode. It should be noted that the translational temperatures of atomic oxygen in the vicinity of the reattachment point were found to be lower in comparison to the other species. This can be attributed to the fact that due to the smaller atomic oxygen diameter, few collisions occur. Therefore, the high-temperature region (between 5000 and 6000 K) does not reach in the upstream direction for atomic oxygen. On the other hand, the vibrational temperature field shows considerable variations due to the difference in freestream initialization given in Table III and the low relaxation rate of the vibrational mode. These results emphasize the nonequilibrium nature of the flow.

Due to high vibrational temperatures, the vibrational energy may not fully accommodate especially for cold metallic surfaces during the short interaction time of molecules with the surface.⁴⁴ Therefore, to model a more realistic gas-surface interaction, an incomplete vibrational accommodation of 0.1 is used. In Fig. 8, the vibrational temperature field shows slight variations in the vicinity of vertex where the flow residence time is the highest. Nonetheless, the overall flow structure is essentially unchanged when Figs. 5(a) and 9(a) are compared. In other words, the effect of the surface vibrational accommodation coefficient is found to be negligible on the size of the separation region (and the surface parameters as discussed below), consistent with the previous work of Ref. 45.

Since the translational temperature was expected to reach a high value of about 5000 K as verified in Fig. 5, the aforementioned chemical reactions were modeled to determine their effect on flowfield and surface parameters. However, when the nonreacting and reacting simulation results are compared, see Figs. 5(a) and 9(b), for the sharp leading edge configuration, the temperature field is essentially the same in the separation and reattachment regions. It should be noted that the translational temperatures for the reacting air case are somewhat lower downstream of the oblique shock emanating from the leading edge because the endothermic chemical reactions remove energy from the flow.

Finally, for case A2, we consider the effect of tick geometry and nonreacting vs reacting air species on mole fraction distributions. Figures 10 and 11 compare mole fractions of each species for sharp and rounded leading edge configurations for condition A2. As shown in Fig. 3, the sharp leading edge configuration has a region of relatively high Knudsen number. This results in increased mean collision times making the diffusion process more predominant. The atomic oxygen species, which has the smallest molecular weight, increases in the expansion region compared to its free stream value, whereas the other species’ mole fractions show the opposite trend. As shown in Fig. 11, these effects are diminished when the geometry with a bevel size of R = 0.2 mm is used due to a decrease in the Knudsen number, as shown in Fig. 3. When chemical reactions are modeled, the dissociation of oxygen molecules results in an increase in the atomic oxygen mole fraction and the exchange region given in Eq. (4) slightly increases the mole fraction of NO while decreasing the concentration of nitrogen molecules in the expansion region, as can be seen by a comparison of Fig. 12(b) with Fig. 10(b). It should be noted that because the effect of the molecular oxygen dissociation reaction dominates over the exchange reaction, the overall mole fraction of atomic oxygen decreases, as can be seen by comparing Fig. 12(a) with Fig. 10(a). The chemical production of NO molecules in the A2 flow is sufficiently small such that reduction of the translation energy of the flow due to chemistry is not appreciable. However, the predicted distribution of NO mole fraction is important for interpreting PLIF measurements.

In our previous work,³ the size and structure of the separation region was found to be strongly dependent on Reynolds numbers. In particular, increasing the Reynolds number results in a larger separation region and a secondary vortex in the vicinity of the double-cone junction. Similarly, for the tick geometry, a higher density case, condition E, provides an opportunity to test the effects of the Reynolds number on the size of the separation region. Figure 13(a) shows that condition E yields a larger separation region in comparison to the lower density case of A2 due to larger adverse pressures on the compression surface. Therefore, the cold temperature region for the higher Reynolds number case of condition E is larger. In contrast to condition A2 shown in Fig. 5(a) where only a primary vortex exists in the separation zone, there is a secondary vortex at the vertex at 1.2 ms for condition E. For the case with a bevel R = 0.1 mm, Fig. 13(b) shows that the size of the secondary vortex is larger in comparison to the sharp leading edge case which has a significant effect on surface parameters and time characteristics of the flow, as will be discussed in Sec. IV. A comparison of the insets of Figs. 13(a) and 13(b) reveals that the presence of the bevel with a radius of 0.1 mm shifts the separation point in the downstream direction. Finally, since the effect of the chemical reactions was found to be negligible on the size of the separation region for the higher stagnation enthalpy of A2, it was not modeled for the lower stagnation enthalpy of condition E. Similarly, a smaller Knudsen number for condition E results in a lower degree of thermochemical nonequilibrium, and therefore, the effect of incomplete vibrational accommodation and diffusion would be negligible at condition E.

After examining the effects of the leading bluntness, an incomplete vibrational accommodation coefficient, and chemical reactions on temperatures and the size and structure of the separation zone, surface fluxes and slip at conditions A2 and E are investigated next. As shown in Fig. 14(a), as the bevel size increases, the surface pressure values also increase in the region from x = 0 to x = 0.03 m since the size of the separation becomes larger as shown in Fig. 5. On the other hand, there is a slight decrease around at x = 0.04 m caused by the fact that the recompression shock moves in the downstream direction, as can be seen in Fig. 14(a). Moreover, once the surface temperature is increased from 300 to 800 K for the bevel size of 0.2 mm, the pressure profile behaves approximately as if the bevel size was increased leading to a larger separation region.

Moving to the surface heating profile, a similar trend is observed as the bevel size is increased, as shown in Fig. 14(b). In the leading edge region, the heat fluxes increase with bevel size. The heat fluxes increase after x = 0.02 m much more rapidly compared to the pressure profile. Once again, the increase in the surface temperature results in an effective bevel size increase between x = 0.02 and x = 0.05 m. However, when the vibrational accommodation coefficient is changed from full accommodation (⁠$αvib$ = 1.0) to incomplete accommodation (⁠$αvib$ = 0.1), the surface pressure and heating profiles remain the same. As shown in Figs. 14(c) and 14(d), the surface pressure and the heating profiles are also the same regardless of whether chemical reactions are modeled. Finally, when the DSMC collision cell resolution and the number of particles are decreased by a factor of 2 (“degraded”) with the time step increased by a factor of 1.6, the surface pressure and heating values remain essentially the same, as can be seen in Figs. 14(c) and 14(d). This shows that the DSMC numerical parameters presented in Table IV generate converged solutions.

Similarly, the effects of the chemistry modeling and vibrational incomplete accommodation coefficient do not appear to influence the velocity slip profile, as shown in Fig. 15(a), for condition A2. The velocity slip reaches considerably high values especially near the leading edge of the sharp tick geometry, associated with the higher Knudsen number as shown in Fig. 3. However, once a bevel of 0.015 mm is introduced, the velocity slip decreases in magnitude in the vicinity of the leading edge since the local Knudsen number decreases with the bevel, as can be seen in Fig. 3(b). As shown in Fig. 15(b), the vibrational temperature jump significantly increases when the vibrational accommodation coefficient is changed from $αvib$ = 1.0 to $αvib$ = 0.1. However, due to a low relaxation rate between the vibrational and translational modes, the translational temperature jump is not affected by the change in the vibrational accommodation coefficient.

For case E, the coefficient of skin friction at the steady state for the sharp leading edge case is found to be in very close agreement with the CFD computations of Khraibut et al.²⁵ for most of the sharp and R = 0.1 mm bevel cases, as shown in Fig. 16(a). Some discrepancies are observed at x = 0.05 m for the sharp leading edge geometry where the leading edge shock interacts with the recompression shock. The maximum percentage difference between the CFD and DSMC data is about 15% which can be attributed to difficulties related with resolving sharp gradients at this location for both approaches. The maxima for the blunt leading edge case have a tendency to decrease and to move downstream. It should be noted that the DSMC data at 2.5 ms and the flow do not reach a steady state, as will be shown in Sec. V, yet they are compared with the steady state CFD solution.⁴⁶ The discrepancy at x = 0.06 m for the bevel case can be attributed to the difference between the time dependent DSMC solution with the steady state CFD. However, overall, the DSMC and CFD solutions are in good agreement.

Finally, Ref. 46 analyzed the effects of velocity slip and temperature jump by performing CFD calculations with and without a slip model. It is interesting to compare the calculated slip values of CFD with DSMC since for DSMC, velocity slip and temperature jump are natural outcomes of gas surface interactions. As shown in Fig. 16(b), the calculated CFD slip values are found to be lower in magnitude as compared to the corresponding DSMC values especially at the leading edge and reattachment locations. The differences in the vicinity of the leading edge especially for the sharp leading edge case can be attributed to rarefied effects. As can be seen in Fig. 3(c) the local Knudsen number is sufficiently large so that a correct prediction of the slip and jump values could be a challenge for continuum methods. The discrepancies at the reattachment point could result from the aforementioned sharp gradients. On the other hand, both methods agree well with other locations and show a substantial decrease near the vertex due to the existence of the secondary vortex. Consistent with the decrease in the local Knudsen number for the bevel geometry, as shown in Fig. 3, the leading edge bluntness has a diminishing effect on the slip values especially at the leading edge. As can be seen in the inset of Fig. 16(b), the velocity slip becomes negative near the leading edge due to high reversed velocities and relatively large Knudsen numbers for the sharp leading edge case. By contrast, for the rounded leading edge, the slip velocities are almost zero at the separation point due to high viscous effects associated with a temperature increase after the leading bow shock as shown in Fig. 16(c). A comparison of the structure of streamlines for the sharp and rounded leading edge reveals that the flow immediately separates after the leading edge whereas it separates further downstream for the rounded leading edge case. Similar to the velocity slip, the estimated CFD temperature jump values underpredict those of DSMC for the sharp and rounded leading edge configuration, as shown in Fig. 16(d).

In our previous work,^3,4,47 the unsteady DSMC signal was postprocessed to obtain the damping rate of the least damped global eigenmode, the amplitude function of the corresponding eigenmode, and the steady-state solution. In Ref. 3, increasing the Reynolds number for hypersonic shock-dominated nitrogen flows over a double cone results in a decrease in the magnitude of the damping rate, which, in turn, results in an increase in the required time to reach the steady-state. A similar approach is used here to determine the time characteristics of massively separated flows with different bevel sizes.

Figure 17 shows the time variation of pressure for different leading edge configurations in the vicinity of the vertex for condition A2. Note that only the pressure field is shown since other flowfield macroparameters such as velocity and temperature fields have a similar time-dependency. It can be seen that the pressure signal for the different geometric leading edge configuration decays in an exponential manner with different rates. This implies that the required time to achieve a steady state solution is different for all cases. In particular, the sharp leading edge case immediately reaches the steady state in less than 0.1 ms whereas the bevel size of R = 0.015 mm case requires about 0.2 ms to obtain a time-independent solution. The larger bevel sizes, however, do not reach a steady state solution even at 1.0 ms which is approximately two times longer than the duration of the experiments for condition A2.

A similar DSMC signal analysis was performed for the higher Reynolds number case, condition E for the sharp and rounded bevel of R = 0.1 mm radius at the leading edge geometries. Figure 18 shows the time variation of the flow field quantities at location 6 near the leading edge separation point where the boundary layer becomes somewhat thicker for the sharp leading edge and the bevel of R = 0.1 mm. Overall, the flow field parameters keep changing in time and the flow does not reach the steady state even for a time period about four times the duration of the experiment for the sharp leading edge, as can be seen in Fig. 18(a). The rounded leading edge case shows a similar time-dependent behavior, and the flow does not reach the steady state at 2.5 ms, as shown in Fig. 18(b). The location where the y-component of the velocity changes its sign can be used to determine the time when the secondary vortex forms. For the sharp leading edge, the secondary eddy can be seen at about 2 ms at location 6, whereas in the rounded leading edge case, it occurs earlier at about 0.7 ms. The pressure values show a strong time dependency, whereas the other flowfield parameters seem to be changing slightly especially for the sharp leading edge case after 2.5 ms.

Based on the probed data at location 6, the formation of the secondary vortex, known as a Moffatt eddy, is captured by the DSMC method at the vertex where the number density is three times lower than the condition E freestream. Running the case further in time reveals that the small secondary vortex seen in Fig. 13(a) becomes larger and moves toward the leading edge of the “tick” geometry in time, as shown in Fig. 19(a). At the location of the secondary vortex, the translational and rotational temperature (not shown) values decrease compared to the earlier time due to additional collisions with the cold wall and with larger density values. Similarly, the vibrational temperature (not shown) is found to decrease in time due to further relaxation of the translational and rotational modes. A comparison of Figs. 19(a) and 19(b) reveals that for the rounded leading edge geometry, the cold region of the secondary vortex is larger. Moreover, the formation of a tertiary vortex at x = 0.025 m can be seen in the inset of Fig. 19(b) and will have a significant effect on the time convergence of the rounded leading edge geometry flow.

By contrast, the sharp leading edge solution almost reaches the steady state at about 5 ms. The size and location of the secondary vortex are found to be in very good agreement with recent CFD calculations^25,46 at the steady state as shown in Fig. 20(a). Similarly, as shown in Fig. 20(b), the translation temperatures are found to be in good agreement, and the small differences at locations A and B can be attributed to the relatively large Knudsen number regions shown in Fig. 3(c). It should be noted that the DSMC solution captures more spatial variation compared to the CFD solution which is somewhat smeared out in the shear layer.

Returning to the time characteristics of the unsteady probed data, in order to correlate the exponential decay rates with the linear stability theory, the residual algorithm²⁹ was used. According to this algorithm, any macroscopic flow parameters, $q=(u,v,Ttrn,Trot,Tvib,P)T$⁠, can be represented by the linear superposition of a steady solution, $q¯(x,y,z)$⁠, and a small 3D deviation from the steady state, denoted as residual, $q̃(x,y,z,t)$⁠, i.e.,

where ϵ is assumed to be small and the solution is near convergence. Assuming that perturbations are 2D and stationary, Eq. (6) can be written in the form of

where $q^r$ is the real part of the amplitude function and σ is the damping rate of the solution. Due to statistical fluctuations of the macroparameters, an exponential curve fitting tool was used to obtain the decay rate of the least damping rate.^3,4,47 The curve fitting results (not shown) indicate that the estimated dimensional damping rates based on the pressure values inside the separation zone (locations 6 and 22) are −24 170 and −22 290 s⁻¹ calculated between a starting sampling time of t₁ = 0.04 ms and the end of the sampling time t₂ = 0.125 ms at condition A2 for the sharp leading edge geometry. The magnitude of the decay rate decreases with leading edge bluntness which in turn increases the required time to reach the steady state. In particular, the decay rates for the leading edge geometry with the bevel of 0.2 mm for condition A2 are found to be −2124 and −2034 s⁻¹ (t₁ = 0.5 ms and t₂ = 1.0 ms) at locations 6 and 22, respectively.

A similar analysis to our previous work³⁸ revealed that the viscous dissipation rates are found to be larger especially in the vicinity of the leading edge and along the shear layer for the sharp leading edge case in comparison to the rounded leading edge case, as shown in Fig. 21. Larger viscous dissipation rates for the sharp leading edge case have a stabilizing effect on the flow. This, in turn, shortens the required time to reach the steady state. Additionally, the size difference of the separation zone can have an influence on the time characteristics of flows. In particular, the development of a larger separation zone takes more time since it is a relatively slow process. Different leading edge configurations of the tick geometry could possibly be used as an active control mechanism especially for design purposes of scramjet engines, similar to Ref. 49 where the impact of the subcavity on the pressure oscillations at different locations was investigated in a supersonic cavity flow.

Similarly, for condition E at the higher Re numbers and for the sharp leading edge geometry, a single damping rate of σ = −650 s⁻¹ is calculated based on the DSMC signal between t₁ = 2.5 ms and t₂ = 5.0 ms and is approximately constant for different locations in the domain and for different flow field parameters, $q=(u,v,Ttrn,Trot,Tvib,P)T$⁠, as can be seen in Fig. 22. Note that only the time variation of pressure, x-velocity, and translational temperature fields at specified locations are shown. The other flowfield parameters are also found to be in a good agreement with the single damping rate and not shown for the sake of brevity. It should be noted that at location 14, the differences between the flowfield parameters, especially for the x-field velocity and fitted curve, are somewhat larger. This could be attributed to the presence of other modes in the SWBLI region. Based on the trend seen for condition A2, it is expected that the magnitude of the decay rate of the least damped eigenmode would decrease further for condition E with a bevel geometry.

After obtaining the corresponding decay rate for conditions A2 and E, the converged steady-state solution can be estimated through the solution of the 2 × 2 system at t₁ and t₂ using Eq. (7). The corresponding amplitude functions for condition A2 are shown in Fig. 23. For the velocity components, a single structure is observed starting from the leading edge and traveling along the shear layer. On the other hand, the translational and rotational fields have a double structure confined to the same region. Due to the coupling between the separation zone and the recompression shock, the amplitude functions are generally found to be larger downstream of the reattachment shock. By contrast, for condition E, the presence of the secondary separation alters the spatial distribution of the amplitude functions considerably, especially inside the separation zone. As shown in Fig. 24, significantly wide regions occur at the primary laminar separation zone labeled as 1, 2, and 4 and in the vicinity of the reattachment point (i.e., location 3). The existence of local maxima in the perturbation velocity, temperatures, and pressure components of the eigenvector in the laminar separation bubble shows all features of a centrifugal instability, as commonly seen in a multitude of flows in which a vortex core is present.⁵⁰ This is the first time to our knowledge that this behavior has been identified and simulated for a hypersonic expansion flow.

Higher amplitude functions, especially in the vicinity of the reattachment point, result in longer time variations of the surface parameters. As shown in Fig. 25(a), the maxima of the surface heating profiles move downstream and increase in time for condition A2 with the bevel of 0.2 mm due to the movement of the triple point in the downstream direction. It should be noted that the time variation between 0.25 and 1.0 ms is discernible in comparison to that between 1.0 and 1.25 ms, showing that the flow still has not reached the steady state. Similarly, for condition E, Fig. 25(b) shows that the maxima of the heat flux values shift in the downstream direction in time. It was observed that the heat fluxes at 4.6 and 5 ms are essentially the same, indicating that the magnitude of the residuals becomes small and one can make a comparison between the steady state CFD and the time accurate DSMC solutions. As shown in Fig. 25(c), the surface pressure coefficient is found to be essentially the same for the DSMC solution at 5.0 ms with CFD whereas DSMC predicts a slightly higher surface coefficient of pressure. Taking the short duration of the experiment into account, which is approximately 1.2 ms, our time-accurate DSMC data at 1.2 ms agree somewhat better with the experiment. This finding suggests that a comparison of numerical simulations at the steady state with the experiment at a time when the flow is still developing can be misleading. Both calculations, however, show some discrepancy with the experiment especially for the prediction of the maxima. In particular, the underprediction of the numerical simulations can be attributed to the instrument limitation of pressure transducers.⁴⁶ 

Hypersonic leading edge separated flows were investigated by using the time-accurate DSMC method for low and high density tick-geometry cases. The effect of the leading edge bluntness was found to be significant on the flowfield and surface parameters for both cases. As the bevel size increases, the size of the separation region was also found to increase whereas the temperature and velocity slip values decrease due to the decrease in Knudsen numbers, especially at the leading edge. The effects of the chemistry modeling and vibrational incomplete accommodation coefficient were found to have little influence on the surface pressure and heat fluxes for case A2. On the other hand, for the same case, a higher surface temperature resulted in a larger separation region.

For the first time, the unsteady behavior of leading edge separated flows with different leading edge configurations was studied by applying the residual algorithm in order to obtain the average damping rate and amplitude function of the least damped eigenmode. The magnitude of the decay rate was found to decrease with the leading edge bluntness which, in turn, increased the required time to reach the steady state for case A2. For the sharp leading edge and lower density case A2, the flow reaches steady state in less than 0.1 ms whereas a value of 1.0 ms, comparable with the run time of the experiments, was insufficient to achieve a steady-state solution for the larger bevel sizes even at these lower densities. The magnitude of the corresponding amplitude functions was found to be larger in the shear layer and downstream of the reattachment point for both cases. The formation of the secondary vortex resulted in larger amplitude functions in the separation zone for the higher density case E, causing larger time variations of the surface parameters. Since the effects of the wall temperature and leading edge bluntness impact the time characteristics of flow, it is important to understand how experiments with short measurement times may affect our conclusions regarding leading edge roughness and variable surface temperatures. Nonetheless, the structure of the streamlines and surface pressure values obtained from the DSMC simulations were found to be in good agreement with CFD at the steady state for the higher density case. However, these flows are dominated by their expansion behavior and whether this conclusion holds for strong shock compression flows still needs to be studied.

The research of O.T. and D.A.L. was supported by the Air Force Office of Scientific Research through AFOSR Grant No. FA9550-11-1-0129 with a subcontract Award No. 2010-06171-01 to UIUC. O.T. and D.A.L. are also grateful for the computational resource provided on ERDC Topaz and Onyx, and AFRL Thunder and Centennial. The work of V.T. was sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under Grant No. FA9550-15-1-0387 Global transient growth mechanisms in high-speed flows with application to the elliptic cone and Grant No. FA9550-17-1-0115 Global Modal and Non-Modal Instability Analyses of Shock-Induced Separation Bubbles, with V.T. as the principal investigator and Dr. Ivett Leyva as the program officer. The authors are also grateful to Dr. Amna Khraibut, Dr. Sudhir L. Gai, and Dr. Sean O’Byrne for proving the experimental and CFD data and useful discussions.