Accuracy and Performance Evaluation of Low Density Internal and External Flows using CFD and DSMC

add1,add2,add3]Surya Kiran Peravali
add2]Vahid Jafari
add3,add5]Amit K. Samanta
add3,add4,add5]Jochen Küpper
add2]Philipp Neumann
add1]Michael Breuer cor1<EMAIL_ADDRESS>
[cor1]Corresponding author

[add1]Professur für Strömungsmechaik,
Helmut-Schmidt-Universität / Universität der Bundeswehr Hamburg,
22043 Hamburg, Germany

[add2]Professur für High Performance Computing,
Helmut-Schmidt-Universität / Universität der Bundeswehr Hamburg,
22043 Hamburg, Germany
[add3]Center for Free-Electron Laser Science CFEL, Deutsches Elektronen-Synchrotron DESY, 22603 Hamburg, Germany

[add4]Department of Physics, Universität Hamburg, 22761 Hamburg, Germany
[add5]Center of Ultrafast Imaging, Universität Hamburg, 22761 Hamburg, Germany

The Direct Simulation Monte Carlo (DSMC) method was widely used, among others, to simulate low density gas flows
with large Knudsen numbers. However,
DSMC has its drawbacks when the lower Knudsen number regime ($\Knud<0.1$) is considered. In this case, approaches from classical computational
fluid dynamics (CFD) relying on the continuum assumption are the methods of choice, providing
accurate solutions at acceptable computational costs. Here, we analyzed the advantages and
drawbacks of DSMC/CFD simulations of internal and external flow fields in terms of accuracy and
computational effort: We used DSMC and CFD
simulations to calculate the flow inside a convergent-divergent nozzle (internal flow) and the
flow around a conical body (external flow). CFD simulations utilized the software OpenFOAM and
the DSMC solution was obtained using the software SPARTA. The results of these simulation
techniques were evaluated in terms of (1) comparison with experimental results, (2)
time-to-solution, (3) energy consumption, and (4) feasibility of hybrid CFD-DSMC approaches.

DSMC; SPARTA; Continuum assumption; Transition regime; Rarefied flow;
high-performance computing

§ INTRODUCTION

The Boltzmann equation is valid at any Knudsen number $\Knud$, i.e., from very low $\Knud$ (hydrodynamic
flows) up to very high $\Knud$ (extremely rarefied gas flows). While flows at low $\Knud$ are
traditionally simulated by methods from computational fluid dynamics (CFD) solving the Navier-Stokes
equations, they fail to provide sufficiently accurate results for rarefied gases since the continuum
mechanics assumption is violated. However, in the rarefied gas flow regime, where the mean free path
of the molecules $\lambda$ becomes comparable to the characteristic length scale of the flow,
micro-scale effects displace the continuum hypothesis, such that particle-based methods often become
the method of choice. One of these methods is the direct simulation Monte Carlo (DSMC)
method [Boyd, 2003], which solves the
Boltzmann equation using Monte Carlo simulations. Unlike the molecular dynamics approach that
compute the trajectory of every particle in the fluid from Newton's equations of motion, the DSMC
method generates collisions stochastically with scattering rates and post-collision velocity
distributions determined from the kinetic theory of dilute gases [Alexander and Garcia, 1997]. DSMC provides
accurate results at high $\Knud$ [Bird, 1994, Shen, 2005].

However, none of these two methods – CFD and DSMC – is entirely
optimal for simulating gas flows in the intermediate-Knudsen-number regime $\Knud\approx0.01\ldots10$, especially when $\Knud$ changes drastically over the flow field. While CFD methods can be very inaccurate due to neglect of molecular effects, the computational efforts for proper
DSMC simulations increases with $\Knud^{-4}$ [La Torre et al., 2009]
such that it becomes computationally very expensive for $\Knud<0.05 $.
Consequently, there has been an increasing interest in using
hybrid methods, especially DSMC/CFD methods to simulate
rarefied gas flows with high Mach and intermediate Knudsen numbers [Wang and Boyd, 2003, Glass and Horvath, 2003, Ghazanfari et al., 2022, Virgile et al., 2022, Malaikannan et al., 2018, Espinoza et al., 2016, Farber et al., 2016, Abbate et al., 2008, Auweter-Kurtz et al., 2005].

Furthermore, to decrease the time-to-solution for the
arising pure or hybrid simulations, high-performance
computing techniques were often applied [Gallis et al., 2014, Plimpton et al., 2019, Li et al., 2022, Ivanov et al., 1997, Roohi and Darbandi, 2012].
While numerous works study the accuracy of the DSMC and CFD methods and how different parametrizations, e.g., the choice and parameters for the collision model or the number of simulation particles impact this accuracy [Fedosov et al., 2005, Alexander and Garcia, 1997, La Torre et al., 2009, Chung et al., 1995, Ghazanfari et al., 2022, Koura and Matsumoto, 1991, Liu et al., 2006, Moss et al., 1997], they do not take into account the cost and efforts of the methods and how much resources in terms of energy (memory and energy) they consume. On the other hand, in studies on
the computational performance of the DSMC method  [Roohi and Darbandi, 2012, Li et al., 2022, Braunstein and Cline, 2003, Gallis et al., 2014, Plimpton et al., 2019, Ivanov et al., 1997] the influence of the different parametrizations on
the simulation results was not investigated.

There is a long list of publications assessing the DSMC and CFD conundrum,
c.f. [Braunstein and Cline, 2003, Khanlarov and Lukianov, 1999, Roohi and Darbandi, 2012, Li et al., 2022, Gallis et al., 2014, Plimpton et al., 2019, Ivanov et al., 1997, Celoria, 2022, Chung et al., 1995, Ghazanfari et al., 2022, Koura and Matsumoto, 1991, Liu et al., 2006, Moss et al., 1997, Fedosov et al., 2005, Alexander and Garcia, 1997, La Torre et al., 2009]. However, a detailed explanation of the
impact of the great variety of parametrization possibilities of the solvers, in particular for DSMC,
is often incomplete or tailored to very specific setups that are not easily generalizable.
Furthermore, the precision of the results and the performance of the DSMC method are not considered
simultaneously to judging if the accuracy is worth the amount of time and energy that the simulation consumes.

Here, we briefly revise the CFD and DSMC methods and compare them in terms of accuracy and
computational effort in different flow regimes. Considering both an external flow around a
conical body and an internal nozzle flow (2D and 3D), parametrizations are detailed and
their impact on the solution is analyzed. The software packages OpenFOAM (CFD) [28] and SPARTA (DSMC) [Sandia National Laboratories, 2022] were used.
Besides accuracy, the computational costs (run-time and energy consumption) are analyzed.

After summarizing related work in <ref>, the
details of the test cases and the
parameters for both DSMC and CFD simulations are described in <ref>. Sections <ref>
and <ref> explain the CFD and DSMC solvers employed in this work, respectively, as well as
the fundamental physical modeling of these approaches. The simulation approaches described in
Sections <ref> to <ref> are underpinned by validation with experimental
results, c.f. Sec. <ref>. A conclusion and an
outlook to future work is provided in Sec. <ref>.

The main goals of this work are:

* Presenting a comprehensive study of rarefied gas flows for internal and external configurations

* Using a combination of 2D/3D configurations of DSMC, CFD and hybrid DSMC-CFD simulation

* Highlighting the effect of different parameters having major impact on the simulation such as
the mesh size, collision model, boundary conditions, and
computational speed

* Validating and evaluating the simulations against experiments

* Reporting the performances and computational efforts of the employed approaches in terms of
scalability and energy consumption

§ STATE OF THE ART

Numerous studies were performed to analyze gas flow fields, including rarefied and continuum flows,
most of them using DSMC or CFD methods. The majority focused on the accuracy and parameter studies
of DSMC simulations such as the choice of the mesh or the collision
model [Swaminathan-Gopalan and Stephani, 2016, Ze, 1705, Weaver and Alexeenko, 2015, Stefanov et al., 2005, Larsen and Borgnakke, 1974, Larsen and Borgnakke, 1975] and validated the results partially against
experimental results. Some studies were related to more practical cases for engineering
purposes [Ghazanfari et al., 2022, Xiao et al., 2014, Wang and Boyd, 2003, La Torre et al., 2009, Glass and Horvath, 2003, Liu et al., 2006, Moss et al., 1997]. Since the DSMC method can be very time
consuming, many studies targeted the performance of the approach: While the computational
performance of the DSMC method on different computer architectures in terms of runtime and energy
consumption is discussed in [Plimpton et al., 2019, Ivanov et al., 1997, Gallis et al., 2014, Khanlarov and Lukianov, 1999, Klothakis and Nikolos, 2015], algorithms and
optimizations to speed up the DSMC method can be found in [Roohi and Darbandi, 2012, Li et al., 2022, Braunstein and Cline, 2003, Celoria, 2022]. To the
authors' knowledge, energy consumption has not been addressed so far. Discussions on the impact of
the mesh dependency, collision model and boundary conditions is not included in the aforementioned

The flow inside micro-nozzles with slip and no-slip conditions are studied using DSMC and a compressible Navier-Stokes solver and the computational results are compared in [La Torre et al., 2009]. The most important outcome is that a better agreement between DSMC and CFD is
observed when a slip wall boundary condition is implemented. It is also shown that the CFD and DSMC
results differ in the divergent part of the nozzle, especially close to the outflow, where the
breakdown parameter $\Knud$ is relatively high or, in other words, the strong rarefied effects start to
appear. No evaluation of the results with regard to a comparison with experimental findings or
computational expenses are included. In [Liu et al., 2006] the effect of different parameters such as inlet
and wall boundary conditions and the Reynolds number on DSMC and Navier–Stokes approaches for a
micro-nozzle flow with a relatively small Knudsen numbers is studied. It has been further
investigated in which part of the micro-nozzle the DSMC and CFD provide the best results. It is
shown that the CFD results exhibit obvious deviations from the DSMC results as $\Knud $ exceeds
$ 0.045 $. The computational performance of large-scale parallel DSMC on homogeneous (CPU) and heterogeneous (CPU+GPU) systems is studied and different programming approaches (MPI, hybrid MPI+OpenMP and OpenACC) are discussed in [Li et al., 2022].

Extensive research on the development of hybrid DSMC-CFD methods
was carried out [Wang and Boyd, 2003, Glass and Horvath, 2003, Ghazanfari et al., 2022, Virgile et al., 2022, Malaikannan et al., 2018, Espinoza et al., 2016, Farber et al., 2016, Abbate et al., 2008, Auweter-Kurtz et al., 2005].
With regard to their application,
simulation cases were mostly very specific, and the comparison was focused typically either on
numerical accuracy [Wang and Boyd, 2003, Glass and Horvath, 2003, Ghazanfari et al., 2022, Farber et al., 2016, Abbate et al., 2008, Auweter-Kurtz et al., 2005, Wu et al., 2007] or on
performance [Virgile et al., 2022]. For example, the advantages of a hybrid DSMC-CFD approach over pure DSMC
is indicated in [Wang and Boyd, 2003]. For this purpose, the authors simulate a hypersonic flow over a
two-dimensional wedge. A comparison of the flow field predicted by a 3D CFD-DSMC simulation
considering a space capsule geometry is given in [Glass and Horvath, 2003], including a validation against
experimental results from wind tunnel tests. The study also investigates the effect of the mesh
dependency on CFD and DSMC methods. A comparison study [Ghazanfari et al., 2022] explores the performance difference between coupled CFD-DSMC and pure CFD methodologies in simulating a gas centrifuge handling $ ^{235}UF_6 $ gas.
Pure DSMC results around the Mars
pathfinder and Mars microprobe capsules are studied in [Moss et al., 1997], abstaining from validation
against experimental results. Computational efficiency of massively parallel (stand-alone) DSMC for
different cases is studied, amongst others, in [Klothakis and Nikolos, 2015] and [Plimpton et al., 2019]. Another comparison study between DSMC and CFD results of a low-density nozzle flow and the experimental evaluation is discussed in [Chung et al., 1995].

The accuracy of a DSMC simulation depends on a number of numerical parameters such as the time step,
the cell size, the number of samples etc., and the choice of the collision model. An analysis of (statistical and deterministic) numerical errors corresponding to numerical
parameters in the DSMC method based on a simple heat transfer problem between two parallel plates in
1D and 2D is provided in [Fedosov et al., 2005]. A further error analysis considering various numerical parameters (sampling cells,
sampling time step and sample sizes) is provided in [Plotnikov and Shkarupa, 2012].

<ref> summarizes which related
studies covered the different areas on DSMC/CFD simulations
in terms of accuracy and computational
performance and classifies the current paper accordingly.
It provides recommendations and suggestions on DSMC,
CFD and hybrid methods as well as a calibration of the DSMC
parameters having a major impact on the simulations to achieve sufficiently accurate results in an acceptable amount of time.

Paper   CM$^a$   $\Knud $    GSI$^b$   Run-   EC$^c$   MU$^d$   Validation

Time         against Exp.

Boyd [Boyd, 2003]

Wang and Boyd, 2003

Glass and Horvath, 2003

Ghazanfari et al., 2022

Gallis et al., 2014

Plimpton et al., 2019

Li et al., 2022

Roohi and Darbandi, 2012

La Torre et al., 2009

Chung et al., 1995

Hedahl and Wilmoth, 1996

Koura and Matsumoto, 1991

Liu et al., 2006

Moss et al., 1997

Khanlarov and Lukianov, 1999

Swaminathan-Gopalan and Stephani, 2016

Weaver and Alexeenko, 2015

Stefanov et al., 2005

Larsen and Borgnakke, 1974

Larsen and Borgnakke, 1975

Xiao et al., 2014

Klothakis and Nikolos, 2015

Plotnikov and Shkarupa, 2012

Pfeiffer et al., 2019

Falchi et al., 2017

Current study

Overview of related work on DSMC/CFD simulations. Here $\Knud$ refers to studies
related to the continuum breakdown using global and/or local Knudsen numbers ($^a$ Collision Model,
$^b$ Gas-Surface-Interaction,
$^c$ Energy Consumption).

§ TEST CASES

§.§ Internal flow – Low density nozzle flow

<ref> indicates the 2D simulation domain of the low density nozzle.
The test case is based on the experiment by Rothe [Rothe, 1970] where the low density flow properties are measured inside the nozzle using the electron beam fluorescence technique. The inflow boundary is located on the extreme left (ab), i.e., at the inlet of the pressure chamber. An outflow boundary condition is assumed along the extreme right  of the geometry (fghij). The segment (aj) in <ref> is the axis of symmetry and for 2D simulations, an axisymmetric boundary condition was imposed. 3D simulations were carried using the complete geometry.

Computational domain of the low density nozzle flow (cases i.I to i.III).

The test gas used in the simulation is nitrogen with a stagnation temperature of 300 K. Three different stagnation pressure configurations are tested according to the setup described in <ref>. That leads to the three internal flow cases i.I to i.III.
The nozzle Reynolds number $Re_N$ is calculated using $Re_N = {\rho_o \, \hat{u} \, r_*}/{\mu_o}$ where $\rho_o$ is the stagnation density, $\hat{u}$ is the adiabatic speed, $r_*$ is the radius of the throat and $\mu_o$ is the  viscosity based on the stagnation condition.
The Knudsen number is calculated based on the stagnation condition and the throat diameter.

Parameters   i.I   i.II    i.III

Test gas    $N_2$    $N_2$     $N_2$

Stagnation temperature, $T_o$ [K]     300      300       300

Stagnation pressure, $P_o$ [Pa]    474      209      141

Wall temperature, $T_w$ [K]    300      300      300

Nozzle Reynolds number, $Re_N$    590      260     175

Throat Knudsen number, $\Knud_t$    $2.3 \times 10^{-3}$     $6.17 \times 10^{-3}$      $9.15 \times 10^{-3}$

Back pressure, $P_a$ [Pa]    1.8      0.64       0.45

Flow-condition parameters
of the low-density-nozzle cases i.I to i.III.

§.§ External flow – Blunt cone & sharp cone

Figure <ref> demonstrates the 2D computational domains of the external flow simulations. For the external flow simulation, two geometries are considered. The first one is a cone-shaped half body with a blunted nose and the second one is a sharp nose. Both geometries are assumed to be infinitely long in positive $x$-direction. The symmetry line of the half body is aligned with the free-stream. For the first half body, only the first 0.05 m from the blunted nose with a radius of $6.35 \times 10^{-3}$ m (0.25 in) is considered and for the second half body, the first 0.09 m from the sharp nose is considered. The computational domain around these bodies is shown in Fig. <ref> (a).

The fluid considered in both cases is pure nitrogen and the free-stream conditions are listed in <ref>.
The ambient values at different altitudes are listed in <ref>.
Other parameters needed for CFD and DSMC simulations such as the pressure and number density were calculated correspondingly.
The flow direction was in  such that the
left and right surfaces of the domain indicate inlet and outlet, respectively, as shown in Fig. <ref>.
For the 2D external flow, half of the computational domain in $y$-direction was considered on the cut through the computational domain ($x$-$y$ plane) as visible in <ref> (b) and (d).
The axisymmetric boundary condition was imposed on the line a-d (see <ref>).

Parameters   e.I   e.II

Geometry    Blunt cone   Sharp Cone

Gas mixture    N2   N2

Ambient velocity, $ V_\infty~[\frac{\text{m}}{\text{s}}] $     2764.5    2072.6

Ambient temperature, $ T_\infty~[\text{K}] $     144.4    42.61

Ambient pressure, $ p_\infty~[\text{Pa}] $     21.91     2.23

Ambient density, $ \rho_\infty~[\frac{\text{kg}}{\text{m}^3}] $     $ 5.113\times10^{-4}$     $ 1.757\times10^{-4}$

Mach number, $M_\infty $    11.3     15.6

Wall temperature, $T_w~[\text{K}] $    297.2    297.2

Knudsen number, $\Knud_\infty$    $4.168 \times 10^{-3}$     $ 1.236 \times 10^{-2}$

Reynolds number, $Re_\infty$    4020     1862

Flow conditions for external test cases.

(a)   (b)

Computational domains for the external flow test cases. (a) Half body cone with blunted nose, 2D setup. (b) cone with sharp nose, 2D setup.

§ NUMERICAL APPROACH

§.§ Continuum Flow Solver:

For this study, CFD simulations were performed using OpenFOAM, version v2112.
OpenFOAM is a finite-volume based code which has a collection of
libraries dedicated to the solution of partial differential
equations (Navier-Stokes equations). Among the standard solvers available
in the OpenFoam library, the  solver is used for
investigating the present test cases.  is a
density-based transient solver used for transonic and supersonic
flow regimes of a compressible gas [Greenshields et al., 2010].

First/second-order schemes are used for the discretization of the governing equations. For gradient and divergence terms, a second-order Gauss linear scheme is used.
For the diffusive terms
in the governing equations, the Gauss scheme is the only choice of discretization and requires an interpolation scheme for the diffusion coefficient (linear). For the transient cases, the temporal derivatives
are discretized using either the first-order bounded implicit Euler scheme or the backward scheme. The Kurganov scheme [Kurganov and Tadmor, 2000, Kraposhin et al., 2015] is used to compute fluxes at the cell interfaces which prevents spurious oscillations around shocks.

In OpenFOAM, the computational domain is generally split up
into a set of patches and the boundary conditions are then
assigned as attributes to the patches and to the field variables
on a patch. Various kinds of boundary
conditions are available in the OpenFOAM library.
The boundary conditions assigned to the test cases
are listed in <ref> and <ref>
for the internal- and external-flow cases, respectively.

Boundary   U   P   T

Inlet    zero gradient     fixed value     fixed value

Outlet    zero gradient      wave transmissive     zero gradient

Wall    no-slip      zero gradient        fixed value

Boundary conditions for the
low-density-nozzle flow (internal flow).

Boundary   U   P   T

Inlet    fixed value     zero gradient     fixed value

Outlet    inlet outlet      wave transmissive     zero gradient

Boundaries     supersonic free-stream      zero gradient        fixed value

Solid walls    no-slip      zero gradient        fixed value

Boundary conditions for the flow around
a conical body (external flow).

§.§ Kinetic Approach: DSMC Method

The Direct Simulation Monte Carlo method is a discrete particle simulation technique that provides a numerical approximation to the solution of the Boltzmann equation. In this method each particle represents a large number of real molecules [Bird, 2013, Shen, 2005] while maintaining the overall phase space distribution. The inter-molecular and the molecule-surface collisions are calculated using probabilistic models.
The momentum term and  the collision term
in the Boltzmann  equation are solved in a decoupled
manner [Alexander and Garcia, 1997]. The DSMC method is a widely used methods
for simulating rarefied gas flows. Its accuracy depends on
the number of particles per grid cell, the size of the grid
cells, the choice of the time steps and the collision models,
i.e., for particle-particle and particle-surface collisions.
It is general practice to limit
the maximum cell size to one third
of the mean free path and the time step
to below one fourth of the
mean collision time [Shen, 2005]. In the subsequent
analysis process, the microscopic-flow properties are
sampled by averaging the particle properties per grid cell
to obtain the macroscopic-flow quantities. These sampled
values are stored for the geometric center of the DSMC grid

We used the SPARTA
(Stochastic PArallel Rarefied-gas Time-accurate
Analyzer [Corporation, 2022]) DSMC code to simulate
transitional and rarefied flows. SPARTA is a highly
benchmarked tool [spa, 2022] that can simulate
systems with a few to millions or billions of particles.
I exhibits a good scalability and memory usage [Corporation, 2022, Klothakis and Nikolos, 2015].

Presently, the no-time-counter (NTC) method is used as
collision-sampling technique [Bird, 2013] along with the
variable-hard-sphere (VHS) and variable-soft-spheres
(VSS) molecular models [Shen, 2005, Bird, 1994] as well
as the Larsen and Borgnakke (L-B) model [Larsen and Borgnakke, 1974, Larsen and Borgnakke, 1975] to handle
internal energy exchange. In the L-B model, only a fraction of the
collisions are assumed to be inelastic which is defined by an
average probability of the internal energy exchange $\phi$.
This parameter is used to determine the rate of the relaxation
process of the energy which can also be given as the inverse
of the relaxation collision number $Z$ ($\phi = 1/Z$) [Chung et al., 1995, Shen, 2005].
SPARTA uses the aforementioned Larsen and Borgnakke model with
constant and variable relaxation [Bird, 1994]. In the current study,
the vibrational mode and chemical reactions are assumed to be
frozen. The collision parameters used are given in

Gas   Molecular diameter   $\omega$   $T_\text{ref}$   $\alpha$   $\xi_{rot}$    $\phi_{rot}$

$N_2$    $4.07 \times 10^{-10}$ m     0.74     273.15 K     1.36     2     0.2

Collision model parameters of simulated gas molecules [Weaver and Alexeenko, 2015, Bird, 2013]. Here, $\omega$ is the power exponent of the temperature in the viscosity law, $T_\text{ref}$ is the reference temperature, $\alpha$ is the scattering coefficient determined from the inverse power law potential [Koura and Matsumoto, 1991] and $\xi_{rot}$ is the rotational degree of freedom.

For the variable relaxation model, the parameter $\phi_{rot}$ in SPARTA
is calculated using equation A.5 in ref. Bird, 1994.
The unknowns in the equation were inferred from experiments [Parker, 1959, Lordi and Mates, 1970].
Both models can be used for temperature ranges from 0 to 1200 K.

The boundary conditions in DSMC can be classified in two types:
Free-stream boundary conditions (inlet and outlet) and wall
boundary conditions, where the gas-surface interactions are
predominant. At the inlet the flux and the thermal state of
the molecules are defined according to the flow conditions
and the particles are introduced into the simulation
domain [Bird, 1994]. The outlet removes the simulation
particles leaving the simulation domain.

§.§ Gas-surface interactions

The gas-surface interaction model plays a dominant role
for the accuracy of the DSMC simulation and choosing an
accurate model is subjected to many factors. The most widely used
gas-surface interaction models in DSMC are
the Maxwell model and the Cercignani-Lampi-Lord
(CLL) model which include several unknown parameters.
We extensively discuss these models to establish some
guidelines in choosing the appropriate model along with
the corresponding parameters.

(a)   (b)

Schematic diagrams of the (a) Maxwell model and (b) CLL  models.

In the Maxwell model (<ref> (a)), it is assumed that
the molecules are either reflected diffusively with complete energy
accommodation ($\alpha_E =100~\%$) or reflected specularly with no
energy exchange. The Maxwell model takes a fraction $\epsilon$ of
incident molecules to be scattered diffusely and the remaining
fraction of molecules ($1-\epsilon$) are scattered
specularly. This model is useful for describing the thermodynamic
behavior of the gas, however, it does not describe the molecular
behavior of the gas which is frequently observed in fundamental
gas-surface scattering experiments [Hinchen and Foley, 1966].

The Cercignani-Lampis-Lord model (CLL) [Lord, 1991, Lord, 1995] is based on the assumption that there is no coupling of the normal and tangential velocity components
during gas-molecule reflections from a surface. Therefore, this model uses two coefficients $\alpha_n$ and $\alpha_t$, which represent energy-accommodation coefficients associated with normal and tangential components of the velocity, respectively. The scattering distribution of the molecules is centered around an average scattering angle $\overline{\theta}_r$ which is a function of the two accommodation coefficients (<ref> (b)). This scattering distribution has a lobular shape similar to the one observed in experiments [Hinchen and Foley, 1966]. This model also accounts the internal energy exchange by introducing accommodation coefficients for rotational and vibrational modes. Furthermore, CLL has the capability to produce diffuse scattering with incomplete energy accommodation ($\alpha_E = \alpha_n = \alpha_t < 100~\%$) by changing the scattering distribution. The details of the implementation of this model are described in [Shen, 2005].

§ HYBRID CFD/DSMC

The one-way coupling of CFD and DSMC was implemented as follows:
The fluid flow in the entire computational domain was simulated by the continuum solver (OpenFoam). An interface position (plane) is chosen in the computational domain and the fluid flow data are extracted at this position. The position is chosen according to a calculated continuum-breakdown parameter from the flow field. The DSMC computational domain was generated by splitting the former computational domain from the interface position toward the outflow exit. The extracted fluid flow data are introduced as an inflow boundary condition to the DSMC simulation. The DSMC simulation is then carried out to resolve the transitional rarefied region.

In order to quantify the continuum breakdown parameter,
different definitions of the Knudsen number are used such as
(a) a global Knudsen number $\Knud = \lambda / L$ and
(b) a local Knudsen number or Boyd's Gradient-Length-Local Knudsen number $\Knud_{GLL,Q} = \lambda |\nabla Q | / Q$ [Boyd, 2003]; $\lambda$ is the mean free path, $L$ is the characteristic length, $Q$ represents a macroscopic flow property such as the density $\rho$, the velocity $\textbf{v}$ or the temperature $T$. A breakdown Knudsen number $\Knud_B$ was calculated based on
the maximum of the local Knudsen numbers and the global Knudsen number
over in the computational domain [Pfeiffer et al., 2019, La Torre et al., 2009], ,

\begin{equation}\label{eq:Kn_break}
Kn_B = \max(Kn, Kn_{GLL,\rho}, Kn_{GLL, T}, Kn_{GLL,|\textbf{v}|}).
\end{equation}

When $\Knud_B>0.05$, the continuum breakdown was assumed [La Torre et al., 2009, Espinoza et al., 2016].

§ RESULTS AND DISCUSSION

§.§ Internal Flow: DSMC, CFD, Hybrid

Simulations were performed for the test cases presented in
<ref> using the continuum approach
(OpenFoam), the DSMC method and the hybrid CFD/DSMC.
The full DSMC simulations are inefficient in the low Knudsen
regime ($\Knud < 0.05$) and demands higher computational effort.
Therefore, the hybrid DSMC method was used in order to
accelerate the simulations. The results obtained by the continuum approach
were used to estimate the breakdown Knudsen number $\Knud_B$
(equation (<ref>)). By determining $\Knud_B$ throughout
the simulation domain, it was observed that the continuum breakdown
occurred in the simulation domain right after the throat.
Hence, the interface between the continuum and DSMC domain was
positioned at the throat using the
approach described in Section <ref>. The continuum simulations
were performed for a 3-D simulation domain
and the hybrid DSMC using both the 2D axisymmetric and 3D
configuration. The simulation results were validated against
experiments from Rothe [Rothe, 1970] and
differences are presented in the following figures.
The sensitivity studies of different simulation parameters
are described in Section <ref>. Due to the limited
availability of experimental data, the cases i.I and i.III
were studied in more detail than case i.II, where only the
centerline temperature data is available.

to <ref> show various flow
parameters obtained using the continuum approach and the
hybrid DSMC approach
and their validation with experimental data [Rothe, 1970].
<ref> (a) and (b) show
the centerline density profiles of the nozzle for cases i.I and
i.III, respectively. Here, the
densities $\rho$ were normalized by the stagnation density $\rho_o$
and this ratio was plotted against the non-dimensional axial distance, i.e.,
the ratio of the axial position from the throat $x$ and the radius of
the throat $R_t$. The radial variation of densities normalized by
the maximum cross-sectional density found at the axis is
studied at the cross-sectional position $x/R_t=13.7$ from
the throat for case i.I in <ref> (a)
and for case i.III in <ref> (b).
The radial density profiles are also studied for other cross-sectional
positions depicted in <ref> in Section <ref>. The density profiles obtained by the hybrid DSMC method show good agreement with the experimental data compared to the continuum method. It should be noted that the error limits in the measured experimental densities reported by Rothe [Rothe, 1970] are $\pm 10~\%$ along the centerline and $\pm 5~\%$ for the relative densities along the cross-sections.

shows the centerline temperature profiles of cases i.I and i.III,
respectively. The equilibrium temperature was calculated by the
continuum method and the translational and rotational temperatures
were predicted by the hybrid DSMC method. Similar to the densities,
all temperatures were normalized by the stagnation temperature $T_o$
and plotted against $x/R_t$. For case i.I, where $Re_N > 500$, the computed
temperatures from all three models agree well with the experimental data,
while for $Re_N < 300$ there are significant differences: Unlike the
case i.I ($Re_N > 500$) where the temperatures decrease monotonically from
the throat to the exit of the nozzle
(<ref> (a)), in the cases
i.II and i.III, $Re_N < 300$, the temperatures reduce to a minimum
at $ x/R_t \approx 6$ and increase toward the exit of the nozzle
(see <ref> (b) and
<ref> (a)). This is a result of stronger
rarefaction effects where the flow gets thermalized due to viscous dissipation, i.e., due to more molecule-surface collisions than molecule-molecule collisions. However, it can be observed that the rotational temperatures predicted by DSMC match best with the measurements as the experimental data [Rothe, 1970] are rotational temperatures.

<ref> shows the normalized
temperature variation against the radial distance normalized by
the cross-sectional radius at
$x/R_t=13.7$ from the throat. The rotational temperatures are again
in good agreement with the experimental data. The rotational
temperature is always greater than the translational temperature
due to flow expansion [Rothe, 1970]. The equilibrium temperature obtained by the
continuum method can be compared with the translational temperatures from the
hybrid DSMC method. The differences between these two temperatures,
particularly for $Re_N < 300$, and the differences in densities observed
above are due to the inaccuracy of the continuum method in the rarefied flow

(a)   (b)

Comparison of density variation along the nozzle
centerline: (a) Case i.I, (b) Case i.III.

(a)   (b)

Comparison of density variation along the radial direction at the cross-section $x/R_t=13.7$ from the throat for the cases (a) i.I and (b) i.III.

(a)   (b)

Comparison of temperature variation along the nozzle centerline
for the cases: (a) i.I and (b) i.III.

(a)   (b)

Comparison of temperature profiles in the radial direction at a non-dimensional distance $x/R_t = 13.7$ from the throat for the cases (a) i.I
and (b) i.III.

§.§ External Flow: DSMC, CFD

External flow simulations are carried out for the test
cases listed in <ref> using pure CFD
and DSMC methods. The flow is assumed to be laminar, which is
in line with the fact that the Reynolds number is relatively
small. As described earlier, the results for the continuum
approach were used to estimate the breakdown Knudsen number $\Knud_B$
(Eq. (<ref>)).
By determining $\Knud_B$ throughout the simulation domain, it is
found that near the surface of the external bodies, $\Knud_B$ is
relatively small ($\Knud_B < 0.004$). This is due to a shock wave
formed around the bodies which eventually increases the pressure
and the density inside the shock layer.
The pressure coefficient $C_p$ is used as a quantity for validating the external flow cases. In <ref> the dotted red lines indicate the measurement line, along which the results are evaluated.

<ref> shows the variation of $C_p$ on the surface plotted
against the axial distance $x$ for cases e.i and e.II, respectively.
Due to the near-continuum flow closer to the surface, both continuum
and DSMC methods show similar trends in the results and reasonably
match the experimental data [Wang and Boyd, 2003, Wang et al., 2002, Moss and Bird, 2004]. However, the DSMC results show
a smoother variation in $C_p$. Furthermore, there is a significant
difference in the results of the methods near the stagnation
region ($x \approx 0$ cm to $0.5$ cm). This is because the DSMC
method predicts a much thicker shock layer around the body compared
to the continuum method particularly near the stagnation
region. This phenomenon is subjected to further study.
Sensitivity studies of different simulation parameters are described in Section <ref>.

(a)   (b)

Comparison of the pressure coefficient $ C_p $ over the surface of
the body for the cases (a) e.I and (b) e.II.

§.§ Impact of DSMC simulation parameters

The results presented in the last two subsections show the level of
accuracy of the DSMC method in resolving rarefied gas flows.
However, as mentioned in Section <ref>, several simulation
parameters are responsible for attaining these accurate results.
In this section, the sensitivity of the results with respect to
different DSMC simulation parameters is presented.
Although the simulations for the internal flow cases were
carried out using the hybrid method, only the DSMC domain,
i.e., the region where $\Knud_B > 0.05$, is considered for this
parameter study.

§.§.§ Effect of simulation particles and DSMC grid size

In order to accurately describe the rarefied gas flow,
it is important to have a sufficient number of simulation
particles in the simulation domain and per grid cell.
These represent the distribution of the actual gas
molecules and are necessary to preserve the statistical
accuracy and resolution of molecular collisions in the
simulation. In SPARTA this property is controlled by the keyword .
The parameter  sets the ratio of real molecules to the
simulation particles. Therefore, the smaller the value of , the
greater the number of simulation particles and the simulation accuracy.
Once the number of simulation particles crosses a certain threshold,
the accuracy of the simulation reaches convergence in results.
Using much smaller values of  compared to this threshold
value increases the computational cost as in DSMC the computational
cost scales linearly  with the number of particles [spa, 2022].
Therefore, it is very important to choose a trade-off value of
in order to optimize the computational costs while ensuring the
accuracy of the simulation.

The simulation domain was discretized with regular grids and the grid cell
size, e.g., $\Delta x$, of the simulation is chosen according to the
criterion mentioned in Section <ref>, i.e.,
$\Delta x \leq \frac{1}{3} \lambda_{min}$ with the minimum
value of the mean free path in the  simulation domain $\lambda_{min}$.

The present convergence study varying  has been performed
for a uniform grid of size $\Delta x$. The simulation particles were
created using the  parameter and were distributed such that
each cell has roughly the same number of particles.
<ref> (a) shows the effect of the
parameter  on the rotational temperature prediction for the
internal flow case i.II. A value of  $\leq 5 \times 10^{15}$
is required to reach convergence and agreement with the
experimental data for a 2D axisymmetric
configuration. Likewise, <ref> (b) shows
the effect of  on the predicted
pressure coefficient $C_p$ for the 2D external flow case e.I.
Here a value of  $\leq 1 \times 10^{17}$ is required
for convergence and a reasonable agreement with the
experimental data. Using smaller values of
increases the number of simulation particles and thereby the
computational cost. Hence, the above mentioned values of
are chosen as a trade-off for the simulation
of the flow cases.

(a)   (b)

Effect of  on (a) centerline rotational temperatures for case i.II. (b) Pressure coefficient $C_p$ for case e.I.

A convergence study has also been performed with respect to the
grid resolution of the simulation. The uniform regular grid is
coarsened until the maximum grid size $\Delta x \leq \lambda_{min}$ is reached
keeping the  value constant (trade-off ).
<ref> shows the
results achieved by the different grid sizes for
case i.II and e.I, respectively. In the limit $\Delta x \leq \lambda_{min}$
there are no significant deviations between the calculations.
Based on these cases the optimal value of the grid size is
$\Delta x \approx \lambda_{min}$ with  $\leq 5 \times 10^{15}$
for the internal and  $\leq 1 \times 10^{17}$ for
the external 2D axisymmetric configurations.
Compared to the previous configuration $\Delta x \approx \frac{1}{3}\lambda_{min}$
these parameters increased the simulation efficiency roughly by a
factor of $3.5$ for the
internal case i.II and by a factor of $4.25$ for the external case e.I.
The 3D simulations are also carried out with $\Delta x \approx \lambda_{min}$ and the corresponding trade-off value of  is $\leq 5 \times 10^{20}$ for the
case i.II and $\leq 2 \times 10^{15}$ for the case e.I,
see <ref>.

The number of simulation particles and the computational effort
could be further reduced without compromising the computational
accuracy by using an adaptive mesh refinement technique. However,
this technique is not considered in the current paper.

(a)   (b)

Effect of grid resolution on (a) centerline rotational temperatures for case i.II and
(b) pressure coefficient $C_p$ for case e.I.

§.§.§ Gas-surface interaction

For the internal flow cases, the simulations are performed using various gas-surface interaction models described in Section <ref>. The surface collisions are treated as either fully specular, fully diffusive or a combination of both. <ref> (a) and (b) show the effect of these surface collision models on the centerline densities for cases i.I and i.III, respectively. In both cases the simulation with a completely diffusive gas-surface interaction matches well with the experimental results. The simulations performed with a specular gas-surface interaction model yield a faster expansion of the flows compared with the experiment and the simulation relying on the diffuse surface interaction. For simulations with the surface interactions modeled by the combination of $50~\%$ diffusive and $50~\%$ specular collisions, the centerline density curve lies in between the completely diffusive and completely specular simulations. As the proportion of specular collisions increases, the curve shifts toward the completely specular simulation curve and vice-versa.

The effect of the gas-surface interactions is also studied regarding radial variations of densities at three different cross-sections of the nozzle. <ref> shows these distributions at the non-dimensional distances of $x/R_t =3.7$ and $x/R_t=6.2$ from the throat for cases i.I and i.III, respectively. The gray shaded regions in the plots represent the error margin in the experiments [Rothe, 1970]. Similar trends to the centerline data are observed in the results, i.e., the simulations with the diffusive-surface-interaction model matches well with the experiment. However, for case i.I ($Re_N > 500$) at the cross-section $x/R_t=3.7$ which is close to the nozzle throat, the gas-surface interaction model with $50~\%$ diffusive and $50~\%$ specular collisions agrees better with the experimental measurements than the completely diffusive model.  This discrepancy is due to the presence of a small inviscid core near the downstream region of the throat. As the flow progresses more in the downstream direction this core vanishes due to the presence of viscous forces. This phenomenon is not present for cases where the nozzle Reynolds number is smaller ($Re_N < 300$) as the viscous forces are present throughout the entire nozzle, , in a viscous boundary layer.

Overall, for regions where the viscous forces are evident the fully diffusive gas-surface interaction model is more suitable. As the flow gets more and more inviscid, i.e., the nozzle Reynolds numbers increases, the gas-surface interactions tend more toward specular collisions. Although not shown, for nozzle Reynolds numbers $Re_N > 1000$, the inviscid core extends till the exit of the nozzle which was also reported in the experimental study [Rothe, 1970]. In such cases, the gas-surface interaction should be biased toward specular. For external cases, e.g., case e.I where the Reynolds number $Re \gg 1000$, the gas-surface interaction model with $10~\%$ diffusive and $90~\%$ specular collisions showed the best agreement with the experiments, see <ref>.

Another important parameter in modeling the gas-surface interaction is the thermal accommodation coefficient. As described in Section <ref>, this parameter quantifies the energy exchange between the surface and the gas. For gas-surface interactions which are fully diffusive, this parameter has a negligible effect on the density. However, it has an influence on the temperature. <ref> (a) and (b) show the effect of the thermal accommodation coefficient on the rotational temperature at the cross-section $x/R_t=13.7$ for cases i.I and i.III, respectively. The thermal accommodation coefficient has a very small effect on the centerline temperature. However, near the wall it has a significant influence. It can be seen that the fully diffusive gas-surface interaction model with an accommodation coefficient of $0~\%$ (adiabatic wall) is in close agreement with the experiment. As the value of the thermal accommodation coefficient increases, the temperatures near the wall increase and diverge from experiments.

(a)   (b)

Effect of gas-surface interactions on densities along the nozzle
axis for the cases (a) i.I and (b) i.III.

(a) Case i.I: $x/R_t = 3.7$   (b) Case i.III: $x/R_t = 3.7$

(c) Case i.I: $x/R_t = 6.2$   (d) Case i.III: $x/R_t = 6.2$

Effect of gas-surface interactions on densities at different cross-sections for case i.I (left column) and case i.III (right column).

(a)   (b)

Effect of gas-surface interactions on $C_p$ for the cases
(a) e.I. and (b) e.II.

(a)   (b)

Effect of thermal accommodation coefficient on the rotational temperature at the cross-section $x/R_w = 13.7$: (a) case i.I and (b) case i.III.

§.§.§ Collision and energy exchange models

More sensitivity studies have been conducted considering the choice of the molecular models (i.e., VHS and VSS) and the energy exchange models (constant and variable relaxation).
According to the evaluation in Section <ref>, the fully diffusive gas-surface interaction model with $0~\%$ thermal accommodation is used in this study for the internal flow cases. The Maxwellian gas-surface interaction model with $10~\%$ diffusive and $90~\%$ specular is used for external flow cases.
<ref> shows the comparison of the various parameters obtained using different molecular and energy exchange models with the experimental data for cases i.III and e.I. It can be seen that the results predicted by these different models lie with in the experimental error estimates. However, in order to rank the best DSMC simulation configuration for internal flows, the accuracy is quantified on the basis of the error estimation given by:

\begin{equation}
\text{error} = \frac{\psi_\text{simulation}-\psi_\text{experiment(mean)}}{\psi_\text{experiment(mean)}} \; ,
\label{eq:error}
\end{equation}

where $\psi $ can be any flow quantity such as the rotational temperature. Again, the collision and energy exchange models have a negligible effect on the density.

The simulation configuration with the variable hard-sphere model with a variable rotational relaxation (obtained using parameters from [Parker, 1959]) yield the best result for internal flow cases. For external flow cases the best result is attained with the variable soft-sphere model with variable relaxation (obtained using parameters from [Parker, 1959]). However, in terms of performance the variable hard-sphere model with a constant rotational collision number (i.e., $\phi_{rot} =0.2$) is the best choice for both internal and external cases as the simulation speed is roughly $1.5$ times higher than the latter. More details on the simulation performance are discussed in Section <ref>.

(a)   (b)

(a) Centerline rotational temperature profiles predicted with different collision models for case i.III. (b) Comparison of $ C_p $ distribution for various collision models for case e.I.

§.§ Performance study

The last part deals with the performance and computational costs of the DSMC simulations running on the high-performance computer HSUper at the Helmut-Schmidt University. HSUper consists of 581 nodes in total, 571 compute nodes with 256 GB RAM per node, 5 nodes with 1 TB RAM per node, 5 nodes with 2 NVIDIA A100 GPUs. Each node contains two sockets equipped with Intel Xeon Platinum 8360Y (36 cores, 2.4 GHz) processors, hence 72 physical cores or 144 virtual CPUs via hyperthreading. The memory is provided by $ 16\times 16 $ GB DDR4 RDIMM 3200 MHz ECC-registered modules. HSUper utilizes InfiniBand HDR, 100 Gb/s, non-blocking fat tree networking.
The test cases with the optimized simulation setup (i.III and e.I) were
used to study the performance.

The  energy consumption of the CPU and DRAM reported in the following are based on the running average power limit [Georgiou et al., 2014] as well as the CPU-time of the job reported by the SLURM workload manager.
In order to keep the result comparable for internal and external flow cases with different configurations, the run-time is and memory usage are  normalized in all cases to 1 core/node. SPARTA uses a hierarchical Cartesian grid, which is not a contour matching mesh. Since a big portion of the grid points is therefore, outside the computational domain and without simulation particles, the usage of load balancing methods is necessary to optimize computational efforts.

<ref> compares the strong scaling results on 1 core up to 144 cores equivalent to 2 nodes on HSUper using static and dynamic load balancing for the internal and external cases. The dynamic load balancing had a huge impact on the speedup of the simulation. The effect of the dynamic load balancing is stronger for the internal case than for the external flow because the external case is a hypersonic flow and due to the very high flow velocity the computational domain gets saturated very quickly with sufficient simulation particles. Thus, the static load balancing at the beginning of the simulation already leads to good results contrary to the internal case where the flow velocity is much smaller than in the external case and the simulation needs therefore a much longer time to reach the steady-state solution.
In <ref> (c) and (d), a similar behavior is also observed in the energy consumption results, which is strongly related to the run-time.

Weak scaling is achieved by varying the parameter , e.g., the number of particles is increased by a factor of three when running on three cores, compared to the single-core case. The results also shown in <ref> are in full agreement with the argumentation given for the strong scaling results. Note that the  parameter decreases correspondingly as the number of cores increases, e.g.,  the  value for 1 core is decreased by a factor of 3 for 3 cores. So the number of simulation particles increases for 3 cores and the number of simulation particles per core remains the same as for 1 core. The  parameter is set for 9, 18, 36 and 72 cores, correspondingly. Note that in <ref> (d) for the internal test case, the weak scaling configuration starts from a different number of particles on 1 core  than the strong scaling configuration ($\texttt{fnum}=5e15$) in <ref> (c).

<ref> demonstrates the strong scaling results on 1 to 16 nodes (i.e. up to 1152 cores). It has to be mentioned that all 72 cores of each HSUper node are occupied for the simulation. Unlike in <ref> (c) and (d), the energy consumption on more than one node increases constantly in <ref> (c) and (d). By using a few processors of a node, the other processors of the node are in the idle mode and still consume some energy resources. Therefore, using more cores of one node which leads to a smaller run-time may require less energy, because the reduction of the computation time may compensate the effect of using more cores. <ref> (c) shows that this is not the case when using more than one node. Considering simulations using 1 node and 2 nodes, if the speedup of the simulation from 1 to 2 nodes would be exactly 2, then the energy consumption would be the same, since the speedup of the simulation would fully compensate the additional energy consumption of more nodes. But as can be seen in  <ref> (a), the speedup of the simulation from 1 to 2 nodes is less than two. So the computational time reduction can not compensate the energy consumption of the additional node and therefore the total energy consumption increases (by a factor of 2). It can be summarized that because an entire node is always allocated exclusively. The faster the simulation is finalized  on it, the less energy is used and this
rationalizes the decrease in energy consumption up to 72 cores. For more cores the
energy consumption increases, because more nodes are allocated. Based on this explanation, in the SPARTA cases it is necessary in terms of energy consumption to always fully use all cores of a node as long as runtime decreases are observed. The weak scaling results in <ref> are in full agreement with the argumentation given for the strong scaling.

(a)   (b)

(c)   (d)

Scaling results on up to 144 cores. Left column: strong scaling.
Right column: weak scaling.
(a) & (b) normalized run-time.  (c) & (d) memory usage. (c) & (d) energy consumption for one time step.

(a)   (b)

(c)   (d)

Scaling results on 1 to 16 nodes. Left column: strong scaling.
Right column: weak scaling.
(a) & (b) normalized run-time. (c) & (d) memory usage. (c) & (d) energy consumption for one time step.

§ CONCLUSIONS

The rarefied nitrogen gas flow through a convergent-divergent nozzle (internal flow) and over a conical body (external flow) were numerically investigated using a continuum-based Navier-Stokes solver and a stochastic Boltzmann solver (DSMC) for different Knudsen- and Reynolds-number regimes. The numerical results were validated against the available experimental data. For the higher Knudsen number range ($\Knud > 0.1$), the DSMC method showed a better accuracy in predicting the state properties such as density, pressure, rotational temperature, etc. However, in the continuum regime ($\Knud < 0.1$) the Navier-Stokes solver showed a higher accuracy and efficiency than the DSMC method. Therefore, for internal-flow test cases having variable Knudsen regimes, a one-way-coupled hybrid approach combining the Navier-Stokes and the DSMC solver was used to accurately resolve the flow in the entire domain and reduce the computational costs.

The sensitivity study of the DSMC simulations revealed several essential points, taking into consideration parameters such as the ratio of real molecules to simulation particles , the grid size, and the surface and molecular collision models. Threshold values for  or the number of particles to approach optimal simulation accuracy were investigated and are reported for both internal and external cases in 2D and 3D configurations.The threshold value of  has been examined for the internal and external cases in 2D and 3D configurations. It was found that different grid sizes $\Delta x$ of the DSMC domain in the range of the minimum mean free path $\lambda_{min}$ have negligible effects on the results. Gas-surface interactions play a major role in obtaining accurate results. For the cases where the Reynolds number is very low ($Re < 500$), viscous effects are dominant and the gas-surface interactions are prone to be diffusive, , the scattering distribution is isotropic. In case of inviscid flow regions, the scattering must be biased toward the specular direction. Furthermore, the adiabatic nature of the wall in the internal-flow case necessitates the use of a diffusive model with an incomplete thermal accommodation, which is achieved using the CLL model. For inviscid regions, a combination of diffuse and specular surface collisions are used according to the Maxwellian model. The sensitivity study on molecular models and energy exchange models did not show major significance on the results and all predictions lie within the experimental margin. However, from the performance point of view, the VHS collision model with constant relaxation is found to have a higher computational speed for all cases.

Performance studies regarding strong and weak scaling
were carried out analyzing the computational speedup
and the energy consumption of the DSMC simulations.
This investigation showed that the dynamic load balancing option of SPARTA provided the most efficient solution. The iterative study shown exploits an ideal parametric setup for hybrid DSMC/CFD simulations without compromising the computational efficiency.
This benchmarked hybrid methodology has a huge
potential for simulating molecular beam experiments,
in which the gas flow passes through a wide
Knudsen number range. At Deutsches Elektronen-Synchrotron
(DESY), molecular beams are generated using an aerosol
injector or a cryogenic buffer gas cell for
nano-particle imaging experiments [Roth et al., 2018, Samanta et al., 2020].
The hybrid CFD/DSMC methodology is currently being
integrated into the simulation framework which
simulates such molecular beam experiments. Furthermore,
the simulation framework can be used to study the
gas-particle interactions in the experiments and also
to optimize the experimental design.

§ ACKNOWLEDGEMENTS

This work was supported by the
HELMHOLTZ Data Science Graduate School for the
Structure of Matter (DASHH, HIDSS-0002), by the Helmut-Schmidt
University, University of the Armed Forces Hamburg, and by
Deutsches Elektronen-Synchtrotron DESY, a member of the
Helmholtz Association (HGF).

We also acknowledge dtec.bw –
Digitalization and Technology Research Center of the
Bundeswehr (project hpc.bw) for provision of computational
resources on HSUper. dtec.bw is funded by the European
Union – NextGenerationEU. Furthermore, we acknowledge
the computational resource on Maxwell operated at DESY.

[Boyd, 2003]
authorI. Boyd, titlePredicting breakdown of the
continuum equations under rarefied flow conditions, journalAIP
Conference Proceedings volume663 (year2003)
pages899–906, doi10.1063/1.1581636.

[Alexander and Garcia, 1997]
authorF. Alexander, authorA. Garcia,
titleThe direct simulation Monte Carlo method,
journalComputers in Physics volume11
(year1997) pages588–593,

[Bird, 1994]
authorG. A. Bird, titleMolecular Gas Dynamics and Direct
Simulation of Gas Flows, publisherOxford, year1994.

[Shen, 2005]
authorC. Shen, titleRarefied Gas Dynamics: Fundamentals,
Simulations and Micro Flows, publisherSpringer,

[La Torre et al., 2009]
authorF. La Torre, authorS. Kenjeres,
authorC. R. Kleijn, authorJ.-L. P. A. Moerel,
titleEvaluation of Micronozzle Performance through DSMC,
Navier-Stokes and Coupled DSMC/Navier-Stokes Approaches, in:
booktitleComputational Science–ICCS 2009: 9th Int. Conf. Baton
Rouge, LA, USA, May 25-27, 2009 Proc., Part I 9,
organizationSpringer, ISBN isbn978-3-642-01969-2,
pages675–684, doi10.1007/978-3-642-01970-8_67,

[Wang and Boyd, 2003]
authorW.-L. Wang, authorI. Boyd, titleHybrid
DSMC-CFD Simulations of Hypersonic Flow over Sharp and Blunted Bodies, in:
booktitle36th AIAA Thermophysics Conference, Orlando, Florida,
USA, pages1–13, doi10.2514/6.2003-3644,

[Glass and Horvath, 2003]
authorC. E. Glass, authorT. J. Horvath,
titleComparison of a 3-D CFD-DSMC solution methodology with a
wind tunnel experiment, journalAIP Conference Proceedings
volume663 (number1) (year2003)
pages441–448, doi10.1063/1.1581580,

[Ghazanfari et al., 2022]
authorV. Ghazanfari, authorM. Shademan,
authorF. Mansourzadeh, titleInvestigation of feed flow
effect using CFD-DSMC method in a gas centrifuge,
journalJournal of Nuclear Research and Applications
volume2 (year2022) pages7–14,

[Virgile et al., 2022]
authorC. Virgile, authorA. Albert,
authorL. Julien, titleOptimisation of a hybrid
NS–DSMC methodology for continuous–rarefied jet flows,
journalActa Astronautica volume195
(year2022) pages295–308, ISSN

[Malaikannan et al., 2018]
authorG. Malaikannan, authorR. Kumar,
authorA. K. Chinnappan, titleA novel efficient hybrid
DSMC–dynamic collision limiter algorithm for multiscale transitional
flows, journalInternational Journal for Numerical Methods in
Fluids volume86 (number9) (year2018)
pages565–581, doi10.1002/fld.4466,

[Espinoza et al., 2016]
authorD. E. R. Espinoza, authorV. Casseau,
authorT. J. Scanlon, authorR. E. Brown,
titleAn open-source hybrid CFD-DSMC solver for high speed
flows, journalAIP Conference Proceedings
volume1786 (number1) (year2016)
pages050007, doi10.1063/1.4967557,

[Farber et al., 2016]
authorK. Farber, authorP. Farber,
authorJ. Gräbel, authorS. Krick,
authorJ. Reitz, authorP. Ueberholz,
titleDevelopment and validation of a coupled
Navier-Stokes/DSMC simulation for rarefied gas flow in the production
process for OLEDs, journalAppl. Math. Comput.
volume272 (numberP3) (year2016)
pages648–656, ISSN issn0096-3003,

[Abbate et al., 2008]
authorG. Abbate, authorC. R. Kleijn,
authorB. J. Thijsse, titleValidation of a hybrid
Navier-Stokes/DSMC method for multiscale transient and steady-state gas
flows, journalInternational Journal for Multiscale Computational
Engineering volume6 (year2008)

[Auweter-Kurtz et al., 2005]
authorM. Auweter-Kurtz, authorM. Fertig,
authorD. Petkow, authorT. Stindl,
authorM. Quandt, authorC.-D. Munz,
authorP. Adamidis, authorM. Resch,
authorS. Roller, authorD. D'Andrea,
authorR. Schneider, titleDevelopment of a hybrid
PIC/DSMC Code, in: booktitle29th Int. Electric Propulsion
Conference, pages1–15, year2005.

[Gallis et al., 2014]
authorM. A. Gallis, authorJ. R. Torczynski,
authorS. J. Plimpton, authorD. J. Rader,
authorT. Koehler, titleDirect simulation Monte
Carlo: The quest for speed, journalAIP Conference
Proceedings volume1628 (number1)
(year2014) pages27–36,

[Plimpton et al., 2019]
authorS. J. Plimpton, authorS. G. Moore,
authorA. Borner, authorA. K. Stagg,
authorT. P. Koehler, authorJ. R. Torczynski,
authorM. A. Gallis, titleDirect simulation Monte
Carlo on petaflop supercomputers and beyond, journalPhysics of
Fluids volume31 (number8) (year2019)
pages086101, doi10.1063/1.5108534.

[Li et al., 2022]
authorJ. Li, authorZ. Zhang, authorL. Deng,
authorX. Ding, authorD. Jiang,
authorP. Wang, titleMassively Parallel Acceleration of
Unstructured DSMC Computing, in: booktitleProceedings of the 2022
6th High Performance Computing and Cluster Technologies Conference, HPCCT
'22, publisherAssociation for Computing Machinery,
addressNew York, NY, USA, ISBN isbn9781450396646,
pages1–6, doi10.1145/3560442.3560443,

[Ivanov et al., 1997]
authorM. Ivanov, authorG. Markelov,
authorS. Taylor, authorJ. Watts,
titleParallel DSMC strategies for 3D computations, in:
editorP. Schiano, editorA. Ecer,
editorJ. Periaux, editorN. Satofuka (Eds.),
booktitleParallel Computational Fluid Dynamics 1996,
publisherNorth-Holland, addressAmsterdam, ISBN
isbn978-0-444-82327-4, pages485–492,

[Roohi and Darbandi, 2012]
authorE. Roohi, authorM. Darbandi,
titleRecommendations on performance of parallel DSMC algorithm
in solving subsonic nanoflows, journalApplied Mathematical
Modelling volume36 (number5) (year2012)
pages2314–2321, ISSN issn0307-904X,

[Fedosov et al., 2005]
authorD. A. Fedosov, authorS. V. Rogasinsky,
authorM. I. Zeifman, authorM. S. Ivanov,
authorA. A. ALexeenko, authorD. A. Levin,
titleAnalysis of Numerical Errors in the DSMC Method,
journalAIP Conference Proceedings
volume762 (number1) (year2005)
pages589–594, ISSN issn0094-243X,

[Chung et al., 1995]
authorC.-H. Chung, authorS. C. Kim,
authorR. M. Stubbs, authorK. J. De Witt,
titleLow-density nozzle flow by the direct simulation Monte
Carlo and continuum methods, journalJournal of Propulsion and
Power volume11 (number1) (year1995)
pages64–70, doi10.2514/3.23841.

[Koura and Matsumoto, 1991]
authorK. Koura, authorH. Matsumoto,
titleVariable soft sphere molecular model for inverse power law or
Lennard-Jones potential, journalPhysics of Fluids A: Fluid
Dynamics volume3 (year1991) pages2459,

[Liu et al., 2006]
authorM. Liu, authorX. Zhang,
authorG. Zhang, authorY. Chen, titleStudy
on micronozzle flow and propulsion performance using DSMC and continuum
methods, journalActa Mechanica Sinica/Lixue Xuebao
volume22 (year2006) pages409–416,

[Moss et al., 1997]
authorJ. Moss, authorR. Wilmoth,
authorJ. Price, titleDSMC Simulations of blunt body
flows for Mars entries: Mars pathfinder and mars microprobe capsules,
in: booktitle32nd Thermophysics Conference,
pages2508, doi10.2514/6.1997-2508,

[Braunstein and Cline, 2003]
authorM. Braunstein, authorJ. Cline,
titleProgress on Parallelizing a General Purpose Direct Simulation
Monte Carlo (DSMC) Code for High Performance Computing Applications, in:
booktitleNASA Publications 20030107284, AMOS 2003 Technical
Conference, Maui, HI, United States, 10 Sep. 2003, year2003.

[Khanlarov and Lukianov, 1999]
authorG. O. Khanlarov, authorG. A. Lukianov,
titleDSMC of the inner atmosphere of a comet on shared memory
multiprocessors, in: editorP. Sloot, editorM. Bubak,
editorA. Hoekstra, editorB. Hertzberger (Eds.),
booktitleHigh-Performance Computing and Networking,
publisherSpringer Berlin Heidelberg, addressBerlin,
Heidelberg, ISBN isbn978-3-540-48933-7,
pages1187–1189, year1999.

[Celoria, 2022]
authorM. Celoria, titlePorting of DSMC to multi-GPUs
using OpenACC, in: booktitleMaster Thesis,
publisherSISSA, year2022.

[28]
authorESI-OpenCFD, titleOpenFOAM - The Open Source
CFD Toolbox, howpublished<https://www.openfoam.com>,
noteversion v2112, year2004–present.

[Sandia National Laboratories, 2022]
authorSandia National Laboratories, titleSPARTA -
Stochastic PArallel Rarefied-gas Time-accurate Analyzer,
noterelease 18Jul2022, year2022.

[Swaminathan-Gopalan and Stephani, 2016]
authorK. Swaminathan-Gopalan, authorK. A. Stephani,
titleRecommended direct simulation Monte Carlo collision model
parameters for modeling ionized air transport processes,
journalPhysics of Fluids
volume28 (number2) (year2016)
pages027101, doi10.1063/1.4939719.

[Ze, 1705]
authorX. Ze, titleDSMC Simulation on Couette Flow in
Micro-channels, in: booktitleProceedings of the 2017 2nd Int.  Conference on Materials Science, Machinery and Energy Engineering (MSMEE
2017), publisherAtlantis Press, ISBN
isbn978-94-6252-346-3, ISSN issn2352-5401,
pages991–994, doi10.2991/msmee-17.2017.192,

[Weaver and Alexeenko, 2015]
authorA. B. Weaver, authorA. A. Alexeenko,
titleRevised variable soft sphere and Lennard-Jones model
parameters for eight common gases up to 2200 K, journalJournal
of Physical and Chemical Reference Data volume44
(year2015) pages023103,

[Stefanov et al., 2005]
authorS. K. Stefanov, authorR. W. Barber,
authorM. Ota, authorD. R. Emerson,
titleComparison between Navier-Stokes and DSMC Calculations
for Low Reynolds Number Slip Flow Past a Confined Microsphere, in:
editorM. Capitelli (Ed.), booktitleRarefied Gas
Dynamics: 24th International Symposium on Rarefied Gas Dynamics, vol.
volume762 of seriesAmerican Institute of Physics
Conference Series, pages701–706,
doi10.1063/1.1941617, year2005.

[Larsen and Borgnakke, 1974]
authorP. S. Larsen, authorC. Borgnakke,
titleRarefied Gas Dynamics 1: Statistical collision model for
simulating polyatomic gas with restricted energy exchange,
vol. volume1, publisherDFVLR-Press,

[Larsen and Borgnakke, 1975]
authorP. S. Larsen, authorC. Borgnakke,
titleStatistical collision model for Monte Carlo simulation of
polyatomic gas mixture, journalJ. Comput. Phys
volume18 (year1975) pages405.

[Xiao et al., 2014]
authorH. Xiao, authorY. Shang, authorD. Wu,
titleDSMC Simulation and Experimental Validation of Shock
Interaction in Hypersonic Low Density Flow, journalThe Scientific
World Journal volume2014 (year2014)
pages732765, doi10.1155/201/732765.

[Klothakis and Nikolos, 2015]
authorA. Klothakis, authorI. Nikolos,
titleModeling of rarefied hypersonic flows using the massively
parallel DSMC kernel "SPARTA", in: booktitle8th Int.  Congress on Computational Mechanics, vol. volume12,

[Wu et al., 2007]
authorJ. S. Wu, authorY. Y. Lian,
authorG. Cheng, authorY.-S. Chen,
titleParallel hybrid particle-continuum (DSMC-NS) flow
simulations using 3-D unstructured mesh, journalParallel
Computational Fluid Dynamics 2006  (year2007)

[Plotnikov and Shkarupa, 2012]
authorM. Y. Plotnikov, authorE. V. Shkarupa,
titleSelection of sampling numerical parameters for the DSMC
method, journalComputers & Fluids
volume58 (number12) (year2012)

[Hedahl and Wilmoth, 1996]
authorM. Hedahl, authorR. Wilmoth,
titleComparisons of the Maxwell and CLL gas/surface interaction
models using DSMC, typeTech. Rep., institutionNASA
Langley Research Center Hampton, VA, United States,

[Pfeiffer et al., 2019]
authorM. Pfeiffer, authorA. Mirza,
authorP. Nizenkov, titleEvaluation of particle-based
continuum methods for a coupling with the direct simulation Monte Carlo
method based on a nozzle expansion, journalPhysics of Fluids
volume31 (year2019) pages073601,

[Falchi et al., 2017]
authorA. Falchi, authorE. Minisci,
authorM. Vasile, authorM. Kubicek,
titleAero-thermal re-entry sensitivity analysis using DSMC and a
high dimensional model representation-based approach, in:
booktitle7th European Conference on Space Debris,

[Rothe, 1970]
authorD. E. Rothe, titleExperimental Study of Viscous
Low-Density Nozzle Flows, typeTech. Rep.,
institutionCAL-AI-2590-A-2, Cornell Aeronautical Laboratory, Inc.,
Cornell University, Buffalo, New York, year1970.

[Greenshields et al., 2010]
authorC. Greenshields, authorH. Weller,
authorL. Gasparini, authorR. J.,
titleImplementation of semi-discrete, non-staggered central
schemes in a colocated, polyhedral, finite volume framework, for high-speed
viscous flows, journalInternational Journal for Numerical Methods
in Fluids volume62 (number1) (year2010)

[Kurganov and Tadmor, 2000]
authorA. Kurganov, authorE. Tadmor, titleNew
High-Resolution Central Schemes for Nonlinear Conservation Laws and
Convection–Diffusion Equations, journalJournal of Computational
Physics volume160 (year2000)
pages241–282, doi10.1006/jcph.2000.6459.

[Kraposhin et al., 2015]
authorM. Kraposhin, authorA. Bovtrikova,
authorS. Strijhak, titleAdaptation of Kurganov- Tadmor
numerical scheme for applying in combination with the PISO method in
numerical simulation of flows in a wide range of mach numbers,
journalProcedia Computer Science volume56
(year2015) pages43–52.

[Bird, 2013]
authorG. A. Bird, titleThe DSMC Method,
publisherAmazon fulfilment, year2013.

[Corporation, 2022]
authorS. Corporation, titleSparta Users Manual,
organizationSandia National Laboratories, year2022.

[spa, 2022]
titleSparta Benchmarks, noteaccessed: 2010-09-30,

[Parker, 1959]
authorJ. G. Parker, titleRotational and vibrational
relaxation in diatomic gases, journalPhysics of Fluids
volume2 (number4) (year1959)
pages449–462, ISSN issn0031-9171,

[Lordi and Mates, 1970]
authorJ. A. Lordi, authorR. E. Mates,
titleRotational relaxation in nonpolar diatomic gases,
journalPhysics of Fluids
volume13 (number2) (year1970)
pages291–308, ISSN issn0031-9171,

[Hinchen and Foley, 1966]
authorJ. J. Hinchen, authorW. M. Foley,
titleScattering Molecular beams by metallic surfaces, in:
booktitle5th International Symposium on Rarefied Gas Dynamics,
pages505, year1966.

[Lord, 1991]
authorR. G. Lord, titleSome extensions to the
Cercignani-Lampis gas-surface scattering kernel,
journalPhysics of Fluids volumeA3
(year1991) pages706–710.

[Lord, 1995]
authorR. G. Lord, titleSome further extensions to the
Cercignani-Lampis gas-surface interaction model,
journalPhysics of Fluids volume7
(year1995) pages1159–1161.

[Wang et al., 2002]
authorW.-L. Wang, authorQ. Sun,
authorI. Boyd, titleTowards development of a hybrid
DSMC-CFD method for simulating hypersonic interacting flows,
journalAIAA Paper  (year2002)

[Moss and Bird, 2004]
authorJ. Moss, authorG. Bird, titleDSMC
Simulations of Hypersonic Flows With Shock Interactions and Validation With
Experiments, in: booktitle37th AIAA Thermophysics Conference,
pages2585, doi10.2514/6.2004-2585,

[Georgiou et al., 2014]
authorY. Georgiou, authorT. Cadeau,
authorD. Glesser, authorD. Auble,
authorM. Jette, authorM. Hautreux,
titleEnergy Accounting and Control with SLURM Resource and Job
Management System, in: booktitleDistributed Computing and
Networking, publisherSpringer Berlin Heidelberg,
addressBerlin, Heidelberg, pages96–118,
doi10.1007/978-3-642-45249-9_7, year2014.

[Roth et al., 2018]
authorN. Roth, authorS. Awel, authorD. A.
Horke, authorJ. Küpper, titleOptimizing aerodynamic
lenses for single-particle imaging, journalJournal of Aerosol
Science volume124 (year2018)
pages17–29, ISSN issn0021-8502,

[Samanta et al., 2020]
authorA. K. Samanta, authorM. Amin,
authorA. D. Estillore, authorN. Roth,
authorL. Worbs, authorD. A. Horke,
authorJ. Küpper, titleControlled beams of
shock-frozen, isolated, biological and artificial nanoparticles,
journalStructural Dynamics
volume7 (number2) (year2020)
pages024304, ISSN issn2329-7778,