DSMC (Direct Simulation Monte Carlo) is a statistical simulation method based on kinetic gas theory. The method, introduced by Bird (Bird, 1970; Bird, 1994), models gas flows by collision and statistical evaluation of test particles to solve the Boltzmann equation. In a simulation, each test particle carries the physical properties of a large number of real molecules advancing through the flow. A key element of the method is the consideration of only binary collisions and the decoupling of collision and motion of the test particles. This can be assumed in the case of a rarefied gas and when the simulation time step Δ t is smaller than the physical collision time t _c . Conservation of momentum and energy is achieved by additional interactions with the boundaries of the numerical domain.

After initialization at the beginning of a simulation, the test particles are moved according to their individual Maxwellian velocities and collisions are performed based on binary and wall collision models. When a cell boundary or open domain boundary is reached, the test particles are transferred or deleted. Sampling of macroscopic flow properties is then performed on the statistical basis of the test particles present in a cell over a number of time steps. These steps are repeated until a steady state is reached.

The mass flow driven inflow model presented and validated in this work was developed within the free software package OpenFOAM (OpenFOAM Ltd, 2016) using the DSMC solver dsmcfoam implemented by Scanlon et al. (Scanlon et al., 2010). The NTC (No Time Counter) scheme (Bird, 2007) is used to evaluate the number of collision pairs in a time step, as well as subcell generation to achieve nearest neighbor collisions within a cell (Bird, 2007). To model binary collisions, the variable hard sphere collision model (Bird, 1994) is used with the addition of Larsen-Borgnakke energy.

Following Bird’s original formulation, the inflow and outflow of the computational domain is modeled based on an equilibrium state formulation (Bird, 1994). In a simulation, the approach describes molecular flux quantities of an equilibrium gas across a surface between the computational domain and an associated virtual volume (Bird, 1994). The approach can be considered velocity driven, since a flow velocity must be set to produce a mass flow normal to the surface in a given direction. It is well suited for geometries in free flow environments, such as a re-entering space capsule (Shang and Chen, 2013).

However, in nozzle flow simulations, this treatment of mass flow at the entry and exit boundaries poses some challenges, since the exact flow velocity must be known in addition to temperature and number density. In contrast, the mass flow entering the engine is often a specified and controlled quantity in experimental setups. Thus, an inflow model based on this quantity would take advantage of the precise knowledge of the inflow mass flow rate as a boundary condition. In Figure 1 the physical inflow and outflow characteristics required by the present mass flow driven inflow model are shown as a flow chart. At the inlet, the mass flow rate must be known from the experiment, as well as the surface area of the inflow and the stagnation temperature of the case. The outflow is modeled by absolute pressure, temperature and velocity at some distance from the nozzle exit.In the model, the inflow particle flux N ̇ in into the numerical domain is modeled as N ̇ in = A in n in , (1) in which A _in is the nozzle’s inlet cross section. The quantity n _in is directly derived from the mass flow rate m ̇ in known from the experiment, the molar mass M of the used gas and the Avogadro constant N _A with n in = m ̇ in A in M N A . (2) to create a steady inflow of molecules across the inlet boundary, all molecules hitting an open domain boundary are additionally reflected. The molecular velocity of the test particles generated at the inlet is based on the Maxwellian velocity distribution for the given temperature T _in. Assuming a subsonic, slow inflow condition at the inflow boundary, the flow velocity is based on the pressure drop across the nozzle and the Maxwellian velocity distribution of the test particles after a certain distance in the downstream direction.

The outflow of the numerical domain is modeled on the basis of Bird’s equilibrium formulation (Bird, 1994) with N ̇ out = A out n out c m p exp − q 2 + π 1 + erf { q } 2 π . (3) Here, c _mp is the most probable thermal velocity of a particle and n _out is the number density at the boundary. The most probable thermal velocity is given by the Boltzmann constant k _B , the macroscopic temperature T (= T _out) and the molecular mass m with c m p = 2 k B T m . (4) the term q in Eq. 3 is given by q = u out c m p cos θ , (5) with the local stream velocity u _out and the angle θ to the surface normal of the outflow boundary.

Experimental data on flow properties in de Laval nozzles for low Reynolds numbers and high Knudsen numbers are hard to find. A well-known case is Rothe’s nozzle experiments (Rothe, 1971), where density and temperature data were measured inside a de Laval nozzle by an electron beam technique, analyzing the emission of the excited nitrogen molecules of the gas flow (Rothe, 1971; Rothe, 1965). The data have been used in a variety of numerical studies (Kim, 1994; Zelesnik et al., 1994; Chung et al., 1995; Ivanov et al., 1997; Arlemark et al., 2012), ranging from the 1970s to the present day, and are chosen as a test case to benchmark and validate the present inflow model.

Further comparison of numerical results of Rothe’s nozzle geometry, obtained by Ivanov et al. (Ivanov et al., 1999), are also considered, since their study includes results for both a Navier-Stokes solver GASP and a DSMC solver SMILE.

Rothe’s study covers several cases of nitrogen de Laval nozzle flows for two throat diameters, different mass flow rates, and varying nozzle pressure ratios between the stagnation chamber inside the nozzle and the background pressure of the vacuum chamber. Each case is specified by a so-called throat Reynold number notation B = 2 h in ρ in μ in r * , (6) where the stagnation chamber gas density and viscosity are ρ _in and μ _in, while the throat radius is represented by r*. The specific enthalpy given by h _in forms the adiabatic escape speed of the case with u ̂ = 2 h in .

The investigated nozzle geometry, shown in the lower part of Figure 2, consists of a converging cone with a half-angle of 30° and an inlet radius of 8.3 mm, connected to a diverging cone with a half-angle of 20°, which forms the supersonic part of the nozzle. The throat of the nozzle is rounded by half the throat radius, which in the chosen case is r* = 2.55 mm. The diverging cone has a length of about 50.7 mm, resulting in a nozzle exit radius of about 21 mm. From this geometry, the numerical mesh shown in the upper part of Figure 2 is derived. As described by Arlemark et al. (Arlemark et al., 2012), inlet and outlet domains are added to the nozzle geometry to account for upstream and downstream effects of the nozzle flow. The inlet domain has a length of 10 mm with no change in radial size. As stated by Ivanov and Markelov (Ivanov et al., 1999) and later by Arlemark et al. (Arlemark et al., 2012), a length of 25 mm and a radius of about 29 mm for the applied outlet domain is sufficient to obtain physical results. The numerical mesh consists of about 79,000 cells. Due to the axial symmetry of the case, a wedge mesh with an opening angle of 2°, 30 cells in radial direction and one cell layer in rotational direction is used. The inflow domain has been refined in the radial direction to improve numerical performance. Note that for a better representation of the mesh, only every 10th line in axial direction and every 5 ^th line in radial direction is shown in Figure 2.

The B = 590 case of Rothe is chosen for simulation because it provides a complete set of experimental results and has also been studied by Ivanov and Markelov (Ivanov et al., 1999). With the given throat radius r* = 2.55 mm, an experimental stagnation chamber pressure of p _in = 473.54 Pa is derived from (Rothe, 1971), while the given pressure ratio p _in/p _∞ = 310 leads to a background pressure of about p _∞ = 1.53 Pa. The experimental stagnation gas temperature is set to T _in = 300 K (Rothe, 1971; Ivanov et al., 1999; Arlemark et al., 2012).

For all present DSMC simulations, the Maxwellian velocity distributions at the entrance and exit of initialized test particles are based on the stagnation temperature T _in. At the walls a temperature of T _w = T _in is applied. Wall collisions are treated by full diffuse scattering of the test particles according to their Maxwellian velocity distribution, while the side planes of the mesh reflect the simulation particles in a specular manner. Binary collisions are modeled by the Larsen-Borgnakke variable hard sphere model (Borgnakke and Larsen, 1975; Bird, 1994). The number of equivalent simulation particles is set to at least seven test particles per cell throughout the nozzle.

In the numerical study, the present mass flow driven inflow model is validated and compared with a velocity driven inflow model based on Bird’s equilibrium state formulation. Thus, two sets of analytically and numerically derived boundary conditions are presented.

For the velocity driven inflow model, temperature, number density and inflow velocity are derived from analytical theory. When analyzing the change of physical quantities over the axial course of the described Laval nozzle, a relation between local velocity and local sound velocity is established, taking into account the conservation of total enthalpy and the definition of sound velocity. With equivalence of the total enthalpy at the nozzle inlet (in), the nozzle throat (*) and the nozzle outlet (out) h in + u in 2 2 = h * + u * 2 2 = h out + u out 2 2 (7) and substitution of the enthalpy depending on the relation of the adiabatic exponent κ, the enthalpy h and the speed of sound, a = κ p ρ = κ h − e = κ c p − c v T = κ − 1 c p T = κ − 1 h − h ref (8) the described relation of velocity and sound velocity at the reference positions results corresponding to the speed of sound far away from the nozzle with u _∞ = 0: a ∞ 2 κ − 1 ⏟ h ∞ − h ref = a 2 κ − 1 ⏟ h − h ref + u 2 2 ⇒ a ∞ 2 κ − 1 = a in 2 κ − 1 ⏟ h in − h ref + u in 2 2 = a * 2 κ − 1 ⏟ h * − h ref + u * 2 2 ⏟ a * 2 2 = a out 2 κ − 1 ⏟ h out − h ref + u out 2 2 (9) accordingly, the relationship between the speed of sound at long range and the local speed of sound as a function of Mach number can be described as follows: a ∞ a = 1 + κ − 1 2 u 2 a 2 = 1 + κ − 1 2 Ma 2 . (10) by means of this relation, the local Mach number at the nozzle inlet and outlet can be predicted as a function of the local sound velocity. a = κ c p − c v T ⇒ T ∞ T = a ∞ a 2 = 1 + κ − 1 2 Ma 2 (11) for thermodynamically ideal gases, isobaric and isochoric heat capacity c _p and c _v results in a direct relationship between the relation of the sound velocities and the corresponding temperature relation.

In order to make more concrete statements about the local change of thermophysical state variables, pressure, density and temperature changes for thermodynamically ideal gases are related to each other as follows: T ∼ p ρ ⇒ ln T T ∞ = ln p p ∞ − ln ρ ρ ∞ + C ⇒ d T T = d p p − d ρ ρ . (12) accordingly, it is known that for isentropic changes of state, as in the Laval nozzle, pressure and density changes are described for ideal gases with 1/T = (c _p − c _v )ρ/p in the power relation to each other described as follows: 0 = ln p p ∞ − c p c v ln ρ ρ ∞ = ln p ρ κ ⇒ p ∼ ρ κ . (13) for thermodynamically ideal gases, the corresponding relation to the temperature change follows: p ∼ ρ T ⇒ T ∼ p ρ ∼ ρ κ − 1 ∼ p κ − 1 κ . (14) hereby, over the entire course of the investigated nozzle, the temperature and density changes result from the axial pressure change depending on the corresponding temperature relation (Eq. 11). With the general equations ρ ∞ ρ = T ∞ T 1 κ − 1 = 1 + κ − 1 2 Ma 2 1 κ − 1 (15) p ∞ p = T ∞ T κ κ − 1 = 1 + κ − 1 2 Ma 2 κ κ − 1 (16) the relative density and pressure change is determined.

In order to be able to make reliable statements regarding the velocity change within the nozzle with varying cross-sectional areas in the axial direction, velocity and density relations must be substituted on the basis of mass conservation as follows: ρ u A = ρ * a * A * ⇒ A A * = ρ * a * ρ u = 1 Ma ρ * a * ρ a = 1 Ma ρ * ρ ∞ ρ ∞ ρ a * a ∞ a ∞ a . (17) from the homogeneity of the total enthalpy, the desired relations result for the local changes of sound velocity a ∞ 2 κ − 1 = a * 2 κ − 1 + a * 2 2 = κ + 1 2 κ − 1 a * 2 ⇒ a * a ∞ = 2 κ + 1 (18) and density following the developped relations: ρ * ρ ∞ = T * T ∞ 1 κ − 1 = a * a ∞ 2 κ − 1 = 2 κ + 1 1 κ − 1 . (19) After appropriate substitution of the nozzle-throat specific relations (Eqs 18, 19) and the general transsonic relations (Eqs 10, 15) to the mass conservation (Eq. 17), the final equation for the local Mach number, A A * = 1 Ma ρ ∗ ρ ∞ ⏟ 2 κ + 1 1 κ − 1 ρ ∞ ρ ⏞ 1 + κ − 1 2 Ma 2 1 κ − 1 a ∗ a ∞ ⏟ 2 κ + 1 a ∞ a ⏞ 1 + κ − 1 2 Ma 2 = 1 Ma 2 κ + 1 1 + κ − 1 2 Ma 2 κ + 1 2 κ − 1 (20) which can only be solved implicitly, is obtained as a function of the ratio of the local cross-sectional area to that at the nozzle throat.

By understanding one of the two positions in the equation as an iteration step of the Mach number, two iteration rules result for all permutations of iteration steps k and k + 1. The first one serves for the development of the Mach number Ma_in < 1 at the nozzle inlet with A = A _in in Eq. 20 as follows: A in A * = 1 Ma in k + 1 2 κ + 1 1 + κ − 1 2 Ma in k 2 κ + 1 2 κ − 1 ⇒ Ma in k + 1 = A * A in 2 κ + 1 1 + κ − 1 2 Ma in k 2 2 κ − 1 κ + 1 . (21) The second iteration rule results from the alternative form of resolution according to the nozzle outlet with A = A _out as follows: A out A * = 1 Ma out k 2 κ + 1 1 + κ − 1 2 Ma out k + 1 2 κ + 1 2 κ − 1 ⇒ Ma out k + 1 = 2 κ − 1 κ + 1 2 Ma out k A out A * 2 κ − 1 κ + 1 − 1 . (22)

The result corresponds to the Mach number at the outlet of the nozzle Ma_out > 1. With a nozzle inlet diameter of 16.6 mm, a nozzle throat diameter of 5.1 mm, and a nozzle outlet diameter of 42.0 mm, the following Mach numbers result with the adiabatic exponent κ = 1.4 at the inlet and outlet: Ma in = 0.05472 (23) Ma out = 6.48801 . (24) both iteration cycles start with Ma⁽⁰⁾ = 1. The inlet velocity u _in = 19.31 m/s results from the sound velocity a _in = 353.07 m/s at the inlet temperature T _in = 300 K.

Ivanov and Markedov (Ivanov et al., 1999) reported that no difference in flow properties is observed in the respective nozzle geometry whether simulation particles are inserted with a radially constant velocity profile or with a physically more correct parabolic velocity profile when using an inlet domain of the given size as in Figure 2. Therefore, in the present study, test particles with a radially uniform velocity are inserted on the left side of the inlet domain, as described by Arlemark et al. (Arlemark et al., 2012).

In contrast to the velocity driven inflow, the boundary conditions for the present inflow model must be set differently. Using Eq. 2, the number of initialized test particles per time step is derived from the mass flow rate, the molecule mass and the area of the inlet. Since Rothe does not give mass flow rates in his studies (Rothe, 1971), the inlet number flux N ̇ in is derived indirectly from the mass flow rate results of two simulations by Ivanov and Markedov (Ivanov et al., 1999) and the received mass flow rate of our own simulation using the velocity driven inflow model. The number density n _in at the inlet is calculated from the steady-state pressure equation p = n k B T , (25) with the given values for stagnation pressure p _in and temperature T _in and the Boltzmann constant k _B . The mass flow rate m ̇ at the inlet may be calculated using ρ in = n in m , (26) for the stagnation density ρ _in and m ̇ in = ρ in u in A in . (27) The vacuum condition at the velocity driven outflow is modeled using the inlet number density value, but with a high flow velocity out of the outlet region. Exceptionally the value has to be larger than the ratio of normal cell length and time step. Assuming that the exit temperature in the diffuser is not greater than the entrance temperature, the exit velocity should not be greater than the product of the expected Mach number at the exit Ma_out (Eq. 24) and the speed of sound at the entrance a _in = 353.07 m/s, so as not to disturb the display in postprocessing. The product of the two actual values would be approx. 2,290 m/s. For the velocity driven simulations, the outlet velocity is set to u _out = 2000 m/s for each of the three outlet surfaces in Figure 2.

Thus, the results for the present inflow model are indicated by the notation A, B, or C. Case A uses the numerical mass flow rate from the GASP solver (AeroSoft, Inc., 1996), while case B uses the mass flow rate from the SMILE solver (Ivanov et al., 1999), both values presented by Ivanov and Markedov in (Ivanov et al., 1999) for the B = 590 test case. For case C, on the other hand, the mass flow rate is calculated from own DSMC simulations, performed with basically the same solver, but with the velocity driven inflow model, using Eq. 27. For a better understanding, the validation process for the present study is shown in Figure 3. The outlet boundary number density is set using Eq. 25 according to the background pressure p _∞ of Rothe’s experiments. Additionally, the macroscopic flow velocity at the outlet is set as described for the velocity driven inflow model. An overview of the used boundary conditions is given in Table 1.

Experimental and numerical results of the rotational temperature for the B = 590 case of the Rothe nozzle are shown in Figure 4. The temperatures of all data are scaled by the experimental stagnation chamber temperature T _in = 300 K. On the left side of Figure 4 are centerline temperatures of Rothe’s experiment, numerical results for GASP and SMILE solvers, and results of four corresponding simulations using the present and the velocity driven inflow models. The radial position is scaled by the throat radius r*, where zero indicates the position of the nozzle throat.

All solvers reproduce the experimental stagnation temperature of 300 K well, with only small numerical scatter in the results. Deviations begin to occur downstream of the nozzle throat with decreasing temperatures. While the first points of the GASP and SMILE results agree well with the present results, the deviations increase near the nozzle exit, which is at about x/r* = 20.

The DSMC solver SMILE reproduces the lowest rotational temperatures compared to the other solvers. While in a single comparison with the Navier-Stokes solver GASP, differences can be explained by the use of no-slip wall boundary conditions for the GASP code (Ivanov et al., 1999), which probably results in higher core flow velocities for the DSMC code and therefore lower rotational temperatures in the nozzle axis, the deviations compared to the results obtained in this paper cannot be explained by this. For the diffuser part of the nozzle, with x/r* > 0, it can be seen that the rotational temperatures predicted by the present simulations agree very well with the experimental data of Rothe. After the nozzle exit plane, the temperature curves of the present simulations drop slightly due to expansion effects of the applied vacuum boundary conditions.

The radial temperature profiles are shown on the right in Figure 4. The axial position of the profiles is near the nozzle exit at x/r* = 18.7. The radial position of the temperature is scaled by the local nozzle radius r _w . Deviations from experimental data are present in all simulations. The DSMC solver SMILE underpredicts temperature values near the centerline, but reproduces the experimental temperature profile quite well. All other solvers show large deviations near the nozzle wall. In particular, in the present simulations, the solvers do not reproduce the correct temperature values for almost 80% of the nozzle radius, and the temperature profile is lifted above the experimental results. This effect is not unexpected and is explained by the treatment of the temperature boundary conditions at the wall, which are calculated to be T _w = T _in. This behavior was also observed by Arlemark et al. in their study (Arlemark et al., 2012).

For both axial and radial temperature profiles, no major deviations are observed, which could be traced back to the three and slightly different mass flow rates of cases A, B and C used in the present simulations.

A comparison of the axial and radial mass density profiles for the discussed solvers is shown in Figure 5. The axial density profiles on the left side of the figure are scaled by the stagnation chamber density ρ _in from Rothe’s experiment, while the axial position is scaled by the throat radius r*. For better comparison the values are logarithmic.

All present simulations reproduce well the centerline density for the pressurized part of the nozzle. There is only a small overprediction of the centerline density. The deviations scale from about 1% to 4% depending on the mass flow rate used for cases A, B, and C. The smallest deviation of 1% is present for case C, which has the lowest mass flow rate. Deviations from the experimental results appear again in the second half of the diffuser. Here, the Navier-Stokes solver GASP overpredicts the experimental data, while all DSMC solvers tend to underpredict the centerline density. Deviations from experimental values can reach 10%, although this value can only be given for the present simulations. Behind the nozzle exit at about x/r* = 20, the mass density continues to decrease due to the suction effects of the applied boundary conditions and the expansion of the flow.

Radial mass density profiles for x/r* = 18.7 are shown in the right part of Figure 5. The mass density is scaled by the centerline value ρ _cl and the radial position is scaled by the nozzle radius at the given position. It can be seen that all solvers reproduce the density very well with only small deviations, while the DSMC results are slightly closer to the experimental ones.

Centerline pressure profiles for the B = 590 case of the Rothe nozzle are shown in Figure 6. All results are scaled by the stagnation pressure p _in from the experiment and plotted logarithmically. The axial direction is again scaled by the nozzle throat radius r*. While the pressure values for the present simulations are obtained directly, they have to be derived from the mass density and temperature for GASP, SMILE and experimental data.

It can be seen that the present simulations reproduce the centerline pressure very well for the pressurized part of the Rothe nozzle. Deviations from experimental data range from 0.5% to 2%. In the diffuser part of the nozzle there are deviations from the experimental data. The Navier-Stokes code GASP reproduces the centerline pressure profile best, with only a small overprediction for the last displayed value. The overprediction at the nozzle exit may be due to the fact that for the GASP simulations in the study by Ivanov et al. (Ivanov et al., 1999), the outlet patch is applied directly to the nozzle exit plane and therefore the numerical deviation is observed. Unfortunately, the good agreement with the experimental values is achieved by underpredicting the temperature and overpredicting the density, see the left sides of Figures 4, 5.

For all DSMC solvers, a slight underprediction of the experimental pressure can be seen in Figure 6. For the present simulations, it can be seen that the deviation is caused by an underprediction of the centerline density, since the centerline temperature values are in excellent agreement with the experimental data. For the DSMC code SMILE, the underprediction of the experimental centerline pressure can be caused by an underprediction of both temperature and density, as seen in Figures 4, 5.