The argon outflow from a conical nozzle into the background gas is considered (Figure 1). The computational domain consists of a chamber, in which the stagnation gas parameters (p ₀ and T ₀) are fixed (domain A), a diverging conical nozzle (domain B), and the region of jet expansion outside the nozzle (domain C). The nozzle is a divergent cone with length L _n = 6 mm, throat diameter d _∗ = 1 mm, outlet diameter d _a = 3.6 mm, and nozzle half-angle θ _a ≈ 12.2°. The throat corresponds to the minimum cross-section of the nozzle’s conical part. In the considered geometry, it coincides with the coordinates’ origin, as shown in Figure 1.

The chamber (stagnation) pressure, p ₀, and the background gas pressure, p _∞ , are varied at close values of gas temperatures in the chamber and background T ₀/T _∞ ∼ 1 (Table 1).

The flow in the nozzle is determined by the Reynolds number Re_∗, calculated from the parameters of the gas in the throat or by the corresponding Knudsen number Kn_∗. The interaction of the outflowing gas with ambient gas is characterized by the Reynolds number Re that is determined by environmental parameters, the characteristic jet size, and the characteristic outflow velocity. Re_∗ and Re are related to each other. When the gas in the chamber and background gas are the same and their temperatures are equal, the additional dimensionless parameters that determine the flow are the expansion ratio n, Mach number at the nozzle outlet, M _a , and Prandtl number, Pr (Dulov and Lukyanov, 1984). The considered flow parameters are summarized in Table 1.

The Reynolds number in the nozzle throat is determined as R e ∗ = ρ ∗ v ∗ d ∗ μ T ∗ , (1) where ρ _∗ and v _∗ are the gas mass density and velocity in the throat (averaged over the section), respectively, and μ(T _∗) is the gas viscosity at the throat temperature T _∗. The Knudsen number (Bird, 1994) can be estimated as K n ∗ = 2 5 − 2 ω 7 − 2 ω 15 π γ 2 M ∗ R e ∗ , (2) where γ = 1.67 is the specific heat ratio, ω is the exponent in the viscosity–temperature dependence, and M_∗ is the Mach number in the throat. For estimation of the main parameters, it is assumed that M_∗ = 1.

The Reynolds number characterizing the flow regime in the expansion region is defined as Re = ρ _∞ v _m L/μ _∞ = αRe _L (Dulov and Lukyanov, 1984). Herein, ρ _∞ and μ _∞ are the density and viscosity of the ambient gas, respectively; v _m is the maximum outflow velocity; L is the characteristic size of the underexpanded jet, estimated as the distance to the Mach disk, x _c , or to the corresponding intersection point of an oblique shock and axis (see Section 3.2); Re _L is defined as R e L = R e ∗ / N , (3) where N = p ₀/p _∞ is the nozzle pressure ratio. The parameter α is a weak function of the Mach number at the nozzle outlet, specific heat ratio, and half-angle and is of the order of unity. Therefore, instead of Re, parameter Re _L is frequently used (Rebrov et al., 1971; Dulov and Lukyanov, 1984). The corresponding Knudsen number (Dulov and Lukyanov, 1984) can be estimated as K n L = λ ∞ / L ≈ 1 / R e L . (4) where λ _∞ is the mean free path of molecules in the ambient gas. Equation 4 assumes an equality of the stagnation and ambient temperatures.

As shown in Table 1, the considering regimes of the nozzle flow are hydrodynamic, while the regime of the jet with a background interaction is transitional in terms of the Knudsen number.

The Bird’s classical algorithm using the no time counter (NTC) scheme (Bird, 1994; Bird, 2013) is employed in DSMC simulations. Elastic collisions are described by the variable hard sphere (VHS) model, for which the collision diameter was determined by the particle relative velocity c as d = d ref c ref / c ω − 1 / 2 , (5) where, for argon, d _ref = 4.17 ⋅ 10^–10 m, c _ref = 558.3 m/s, (T _ref = 273 K), and ω = 0.81 (Bird, 1994). This equation implies that μ ∼ T ^ω (Eq. 6).

The dimensions of domain (A) are L ₀ = 2d _∗ and H ₀ = 1.35d _∗, and the dimensions of domain (C) are L _e = 330d _∗ and H _e = 45d _∗. At the boundaries 0 and 1, the gas is assumed to be at rest. Diffuse reflection with total energy accommodation and temperature T ₀ is assumed at boundaries 2, 3, and 4. The conditions at boundaries 5 and 6 are similar to the “gas at rest” condition: incident particles are removed from the simulation, and the particles flowing into this region correspond to a semi-Maxwellian distribution function with parameters p _∞ and T _∞ .

A zone mesh with a different number of cells within each zone is used in the simulations. The zones are of a rectangular or triangular shape and consist of rectangular and triangular cells. The DSMC mesh for the whole domain and nozzle region is shown in Supplementary Figure S1. The number of zones varies from 272 to 400, the total number of cells varies from 10⁶ to 3.5 ⋅ 10⁶, and the number of computing particles varies from 5 ⋅ 10⁸ to 10⁹, depending on the simulation case. Additional verification simulations, in which the number of cells in the nozzle and the area adjacent to it varies, showed the practical absence of dependence of the solution with a further decrease in cell size. The in-house DSMC code (Bykov and Gorbachev, 2019) is employed, and the algorithm is parallelized by the spatial decomposition method and performed on the Polytechnic Tornado cluster of the supercomputer center of the Peter the Great St. Petersburg Polytechnic University.

The code verification was performed in several previous works, such as in Bykov and Gorbachev (2019) and Bykov and Zakharov (2022). Herein, we present a comparison of flow through nozzle DSMC results against solution results of the NS equations.

Supersonic jet modeling was performed using the OpenFOAM package, version 2112, distributed under the GNU GPL license. A set of equations describing the transient flow of viscous compressible gas was used—the mass conservation equation (continuity equation), momentum conservation equation (motion equation), and energy conservation equation—which were closed by the ideal gas law and Newtonian law for viscous stress tensor (Greenshields et al., 2009).

The power law of dependency of viscosity on temperature for argon was used as follows: μ = A T ω , (6) where ω is as given in (5), A = 2.25 ⋅ 10^–7, and the viscosity units are Pa⋅s.

An unsteady viscous gas flow solver rhoCentralFOAM was used. It is based on the Kurganov–Tadmor central-upstream difference scheme. The rhoCentralFOAM solves momentum and mass conservation equations using the time separation approach. The inviscid equations are solved explicitly to obtain the predicted variable value. Then, the diffusion terms are added in the form of implicit corrections to the original inviscid equations (Greenshields et al., 2009).

The computational domain for NS simulations was taken to be the same as that for DSMC (see Figure 1). At the inlet and outlet sections (borders 0 and 6), the soft boundary conditions are used as ∂ V n / ∂ n = 0 , ∂ p / ∂ n = 0 , ∂ T / ∂ n = 0 . (7) where n ⃗ is the unit normal vector of the surface, V ⃗ is the gas velocity, and V _n is the normal component of the surface velocity vector. On the walls (1, 2, 3, and 4), isothermal with T ₀ temperature and no-slip boundary conditions are applied as T = T 0 , V ⃗ = 0 . (8) Although the conditions at boundary 1 inside the chamber between the DSMC and NS approaches differs, they lead to the same mass flux through the throat.

The outer boundary (5) represents a far field with fixed pressure. In the case of an inflow, the temperature boundary condition is switched from a zero normal gradient to a fixed value at T = T _∞ : O u t fl o w : ∂ V n / ∂ n = 0 , p = p ∞ , ∂ T / ∂ n = 0 , (9) I n fl o w : ∂ V n / ∂ n = 0 , p = p ∞ , T = T ∞ . (10)

The problem is solved using the transient approach, where at the initial moment, discontinuity in the pressure field is assumed: a high pressure p ₀ at chamber (A) and a low pressure p _∞ in nozzle (B) and expansion region (C). During computation, steady fields are obtained.

A 2D structured mesh with approximately 200,000 cells with near-surface layer resolution is used. The computational mesh is shown in Supplementary Figure S2. To see the cell size effect, additional simulations are performed, which show that the solution obtained on the 200,000 cells mesh is not dependent on the cell size.

The simulation results, of both the DSMC and NS, show a non-monotonic behavior of the gas dynamic parameters inside the nozzle (Figures 2, 3). In the paraxial region of the nozzle, the density and temperature maxima and the minimum Mach number are observed. The behavior of the parameters in the flow regions near the nozzle wall and near the axis significantly differs. These features of the flow are due to the following factors.

First, an oblique shock is formed inside the conical nozzle (Ivanov et al., 1999). A sharp expansion of the flow directly behind the throat and further interaction of the flow with the wall leads to the appearance of a positive pressure gradient along the wall near the throat. A local increase in pressure near the surface leads to the formation of a shock wave. With the increase in chamber pressure, the intersection zone of shock with the flow axis shifts toward the nozzle outlet (Figure 3). The intensity of shock also increases. In cases 2 and 3 (inside the nozzle, cases 1 and 2 are identical), for which the stagnation pressure in the chamber differs by the factor of 3, the ratio of the Mach numbers on the axis behind and before shock is approximately 1.22 (case 2) and 1.78 (case 3), respectively. After the shock, the axial parameters inside the nozzle vary monotonically; the density and temperature decrease, while the velocity and Mach number increase. In the general case, shock reflected from the axis can reach the nozzle wall, be reflected from it, and “return” to the downstream axis.

Second, the flow in the near-wall region is affected by the boundary conditions on the wall. In the axial region, the temperature decreases with increasing distance from the throat, except at the vicinity of the shock. There is practically no temperature gradient along the wall surface. The gas temperature is close to the wall temperature. At the wall, the gas decelerates and a subsonic flow region is formed. The sound line starts at the top of the throat and closes at the nozzle edge.

Despite the rather high values of the Reynolds number, there is no region of one-dimensional isentropic expansion, and viscous effects manifest themselves in the entire region of the flow through the nozzle.

For the Re_∗ numbers given in Table 1, the flow regime through the nozzle is close to hydrodynamic. This makes it possible to carry out a calculation using the NS equations and compare the modeling results obtained by different approaches. The axial distributions of the parameters obtained by solving the NS equation for case 2 are shown in Figure 3, and an almost complete agreement of the results is seen.

The flow in a supersonic underexpanded jet for high Reynolds numbers Re _L has been widely studied (Dulov and Lukyanov, 1984; Avduevski et al., 1989; Rebrov, 2001; Franquet et al., 2015; Sutton and Biblarz, 2017). Depending on the value of Re _L , the flow in the jet expansion region can be laminar or turbulent. In the initial section behind the nozzle, a barrel shock is formed, which confines the region of free expansion of the jet. At small expansion ratio values, the jet consists of some periodically repeating structures, such as the barrels (Mach diamonds). This mode is characterized by a regular reflection of a barrel shock from the axis (no Mach disk). The pressure ratio increases above a certain value n _∗, which leads to the appearance of a central shock wave, the Mach disk. An increase in n leads to the gradual vanishing of all barrels, except the first one. The shock wave structure includes the barrel shock, Mach disk, and reflected oblique shock (Supplementary Figure S3A). The distance to the Mach disk at high Reynolds numbers Re _L can be estimated using the Lewis–Carlson expression (Lewis and Carlson, 1964) as follows: x c / d a = 0.69 M a γ n , (11) where M _a is the Mach number at the nozzle outlet. Equation 11 is the result of experimental data processing for nitrogen, helium, and carbon dioxide jets. For sound outflow through an orifice in the hydrodynamic regime, the following expression (Ashkenas and Sherman, 1965) can be used: x c / d ∗ = 0.67 N . (12) The last equation is convenient because it does not require a knowledge of the parameters in a nozzle outlet section.

With a significant decrease in the Re _L number down to approximately 100, two trends are observed. The first trend is in the thickening of the barrel shocks and Mach disk, as well as the jet gas and ambient gas mixing layer at the jet boundary. The flow is restructured from the laminar regime to the scattering regime (Muntz et al., 1970), in which the observed flow in the far field of the jet acquires a diffusion character.

The second trend is in the evolution of a shock wave structure. In the experiments with a decrease in Re _L below a certain value (e.g., Re _L ≈ 350 for CO₂), a transition to a quasi-regular reflection of an incident oblique shock, the so-called X-shaped configuration, is observed (Kisliakov et al., 1975; Rebrov, 2001) (Supplementary Figure S3B). In the region of the X-shaped configuration, there is a sharp increase in density that is relative to the maximum, which is possible for a direct shock in an inviscid gas.

Herein, we present the simulation results of an underexpanded jet outflow for extremely low Re _L (Re _L ≤ 30, Table 1).

For case 3, Re _L = 30 (Kn _L = 0.033), the density field is shown in Figure 4. The formation of a shock structure is clearly seen (Supplementary Figure S3C). The jet core is confined by the oblique shock. Shock is smeared (Figures 4, 5, 6) and located close to the jet axis at a distance that is less than Y < 10d _∗ = 2.8d _a . The influence of the reflected shock on the flow is limited by the radial coordinate Y ∼ 20d _∗. In the expansion region, the density minimum is close to the environment density. A minimum on temperature distribution and maximum on Mach number distribution along the axis are observed in the region before the shock–axis intersection region. The temperature and density minima and Mach number maximum positions do not coincide. According to the data presented in Figure 6, the Mach number reaches its maximum value on axis M ∼ 22 at the position X/d _∗ ∼ 50, then decreases to value M ∼ 3 behind the shock. The flow behind the shock on the jet centerline is supersonic. The maximal value of M does not exceed 5. The oscillations are typical for the axial pressure distribution. The pressure change in the shock region exceeds two orders of magnitude. Starting from X/d _∗ ∼ 150, the pressure changes slightly. The density, velocity, and temperature changes in the expansion region are also non-monotonic (Supplementary Figure S4). The density in the shock–axis intersection region increases by an order of magnitude, and the velocity undergoes a smaller change, decreasing by 10%.

Due to rarefaction, the shock waves are smearing and no more the infinitely thin layers. For the considered flow case, the presence of a small radius Mach disk (Supplementary Figure S3A) is questionable. The shock wave configuration obtained in the simulations is very close to the X-shaped configuration, as previously described by Kisliakov et al. (1975) and Rebrov (2001) (Supplementary Figure S3B). The latter assumes an absence of the Mach disk, a quasi-regular intersection of the incident oblique shock with the flow axis and supersonic flow behind the X-point. A reflection of the incident shock from an axis should lead to the presence of the two density maxima of the longitudinal distribution near the axis that correspond to the incident and reflected shock, respectively. In Figure 5, at the coordinates Y/d _∗ = 3.5 and 5, three density maxima are observed. The first (number 1 in the figure) corresponds to the density change when the line with Y/d _∗ = 3.5 or 5 intersects the beginning of the oblique shock layer near the nozzle. The other two maxima (2 and 3) are the aforementioned line intersections at the end of the incident shock and reflected shock (near the interaction region), respectively. As Y/d _∗ decreases, starting from Y/d _∗ ∼ 1.5, maxima 2 and 3 degenerate into one maximum 4. For Y/d _∗ = 0.5 and 1, the density distributions in the shock region on the jet centerline practically coincide, demonstrating only one density maximum.

One of the goals of this work was to analyze the possibility of using the empirical equations for estimating the location of the intersection of shock and the axis. The location of shock may be predicted by using equation (11) or (12). The corresponding lines (marked as L1, L2, and A) are depicted in Figure 6 (also see Supplementary Figure S4). Equation 11 requires nozzle outlet parameters M _a and p _a . For the considered cases, such parameters vary by the radial coordinates at the nozzle exit. The parameters can be estimated as axial parameters (line L1 in Figure 6) or averaged parameters (line L2). The shock location can be defined as the position of the maximal density gradient within the shock. Equation 12, which is valid for sound outflow in the hydrodynamic regime, underestimates a shock location coordinate. The same situation exists with the assessment by (11) using the average parameters at the nozzle outlet (line A). On the contrary, the assessment by (11), with the use of axial parameters, overestimates a shock location coordinate (Supplementary Figure S4). It can be seen that all the shock location predictions are not satisfactory for the considered parameters’ range.

With a threefold decrease in the Reynolds number to Re _L = 9.9 (Kn _L = 0.1, case 2), the flow pattern changes significantly. The shock structure is no longer observed (Figure 7). However, the isolines of the parameters are significantly deformed in relation to the case Re _L → 0. The deformation of the isolines is the most pronounced for the density field (Figure 7). Thus, we can conclude that a smeared and degenerate shock wave structure exists for this situation. This statement somewhat corrects the conclusion (Dulov and Lukyanov, 1984) about the complete degeneration of such a structure at Re _L = 10. A plateau is observed on the axial density distribution corresponding to the coordinates 40 < X/d _∗ < 60 (Supplementary Figure S5), which can be identified with the interaction of smeared oblique shock with the flow axis. The Mach number reaches its maximum value M ∼ 13 at the location X/d _∗ ∼ 30 (Figure 8). Furthermore, the Mach number, without experiencing oscillations, decreases monotonically to zero at the boundary of the computational domain. The pressure in the jet becomes equal to the pressure of the background starting from X/d _∗ ∼ 100, and the density and temperature reach their asymptotic values at X/d _∗ ∼ 200 (Supplementary Figure S5). For this case, we can only conventionally discuss the position of shock on the axis and identify it with a plateau on the axial density distribution.

Further reduction of the Reynolds number to Re _L = 5 (Kn _L = 0.2, case 1) leads to a complete shock wave degeneration (Figure 9). Density field isolines’ deformation is no longer observed. In contrast to the outflow into vacuum, the Mach number and temperature fields have a maximum and minimum, respectively, but no oscillations of the parameters are observed. The gas dynamic parameters’ axial distributions are shown in Figure 10, Supplementary Figure S6. The plateau on the density distribution disappears, which confirms the complete degeneration of the shock waves. The distance to shock may be estimated by (11) or (12), but it does not make sense for this case.

The change in the shock wave structure in a strongly underexpanded jet at low Re _L numbers can be clearly seen in Supplementary Figure S7. In case 1, the density decreases monotonically; in case 2, a step appears on the density distribution; and in case 3, the structure shown in Supplementary Figure S3B becomes pronounced.

The calculations for all the cases were performed on the basis of the NS equations as well. The corresponding computational parameters and boundary conditions are described in Section 2.3.

Figures 11A, B show the axial parameters’ distributions in the jet, which include the flow deceleration region, for case 1 (Re _L = 5, Kn _L = 0.2). In the initial part of the jet, the solutions based on the NS and DSMC methods coincide with each other. However, the NS approach predicts a larger size for the initial jet region up to the point where the flow deceleration begins, which corresponds to the maximum velocity in the jet. Accordingly, the Mach number limiting values achievable in the jet, predicted on the basis of the NS equations, turn out to be higher, and the limiting density and temperature values are lower in relation to that predicted by the DSMC method. The thickness of the layer in which the flow deceleration occurs turns out to be smaller than that predicted by the NS equations. This conclusion is consistent with the one proposed by Greenshields and Reese (2007) and Scanlon et al. (2015), where the discrepancies between the experimental data and NS solution for the problem of the shock wave thickness determination in a rarefied gas are discussed. A larger density gradient in the shock wave obtained from solving the NS equations, with respect to the numerical solution of the Boltzmann equation or the DSMC results for Mach number M ≫ 1, was also observed in the study conducted by Shoev et al. (2020) and Poleshkin and Kudryavtsev (2023). The density distribution described by the NS equations clearly shows a low-intensity jump, which is not typical for the considered case (Re _L = 5 and Kn _L = 0.2), corresponding to the scattering regime (Dulov and Lukyanov, 1984; Avduevski et al., 1989). As discussed previously, the shock wave degenerates with Reynolds number reduction. Case 1 corresponds to the intermediate situation when it is completely smeared when simulated within the DSMC method, but not yet smeared when simulated by NS equations. Thus, in this case, there is a qualitative difference between the results of the methods used.

As Re _L increases, the situation changes (Supplementary Figure S8). At Re _L = 30, there is a qualitative agreement between the results of different methods. Both methods reproduce the shock structure that is typical for the considered parameters. However, as previously discussed, NS predicts higher Mach numbers.

The reason for the discrepancy between the NS equations and DSMC method solutions can be related with the high degree of non-equilibrium flows within the shock layers. Supplementary Figure S9 shows the distributions of the translational temperature components on the axis for cases 1 and 3. In case 1, the translational non-equilibrium turns out to be high and reaches the values R_n = (T _x − T _y )/T ∼ 50% in the jet deceleration region, where T _x and T _y are the x- and y-temperature components that are the mean values of the corresponding molecular kinetic energy components, determined by the molecular velocity components v _x and v _y , respectively. Within the shock layer for case 3, similar to the previous case, R_n increases up to ∼ 40 − 50 % . The distribution function in a rarefied shock layer can significantly differ from the Maxwellian one (Bykov and Lukyanov, 1998). The condition for the correct use of NS equations is some relatively small deviations from the equilibrium state (Cercignani, 1975).