Simulation and analysis of physical characteristics of a molecular beam in a collimator using methods of kinetic theory

The paper is dedicated to the study of gas outflow from a vessel for the case of a simple collimator and to the determination of molecular beam and collimator characteristics enabling a maximum reduction in the beam width. A software system was developed for studying gas flow in various geometries using the GP-GPU.


Introduction
Currently, the following two methods are used for simulation of rarefied gas flows -the method of direct statistical simulation (Monte Carlo) [1] and of the finite-difference solution of the Boltzmann equation. In order to reduce computing requirements, the Boltzmann equation is sometimes replaced with its relaxation model having close physical meaning [2]. Relaxation processes are considered in various fields of science and technology [3][4][5][6][7]. Similar relaxation processes can occur in different systems [7][8][9][10][11][12][13] and surfaces [14][15][16].
The first method dominated over the other one for a long time and is still quite popular due to its low computing resources requirements and mathematical simplicity. However, it has essential drawbacks -low accuracy and difficulties in calculating low-density regions and flows at low Knudsen numbers and low Mach numbers.
Rapid development of computational capacities and mathematical tools gave a strong impetus to evolution of the second method. A conservative method of calculating the collision integral was developed [17][18][19], which allowed solving the Boltzmann equations efficiently [20] and removed the last obstacles to the use of this method in practice. Development of methods and technologies of distributed and parallel computing on clusters and graphical processors significantly increased the available capacities of computers and computing systems. In view of the above, the direct solution method was used in this work for the gas simulation.
Gas outflow from a vessel of fixed volume is considered in this work. Behavior of a molecular jet with low Knudsen number for the case of a monatomic gas, and the influence of the simple collimator on the jet width are simulated and studied. The calculations are performed in a graphical processing unit using the Nvidia CUDA technology [21,22].

Statement of the problem
To study the influence of the collimator on the gas jet we will consider a simple installation -a slotted rectangular vessel and a simple collimator. The installation is schematically depicted in Figure 1.
The installation is effectively two-dimensional since its parameters are independent of position on the axis perpendicular to the plane of the schematic. This simplification is made in order to reduce the computational capacity required for solution of the problem.
Initially, the vessel is uniformly filled with gas. There is no gas in the rest of the computational region. Upon collision with the vessel walls or collimator, gas molecules scatter in random directions. When the molecules get on the boundary of the computational region, they leave it forever.
As the system evolves, gas flows out of the vessel and forms a jet. This jet encounters the collimator in its path, and geometrical characteristics of the jet change.
However, it is not so easy to determine these changes numerically. Indeed, the jet has no clear boundary, and thus the definition of the jet width shall be refined. Let us use the fact that the problem has a symmetry axis and define half-width of the jet. Definition 1. Half-width of the jet is the distance from the symmetry axis at which the concentration of particles is half the concentration on the axis. Denoting the half-width by R and the coordinate along the symmetry axis by X, we can naturally define the jet spread.
Definition 2. Jet spread is dR/dX = tan α--slope ratio R(X), wheeα is the jet slope. It is the spread that will be a measured geometrical characteristic of the jet. For its calculation, it is necessary to calculate the concentration of particles in all points of the vessel.

Mathematical basis
An unknown function derived from the Boltzmann equation is a distribution function defined such that The concentration of particles is determined from the distribution function To find the distribution function, it is necessary to solve Boltzmann equations with specified initial and boundary conditions [2]: where the right-hand side is the collision integral: For efficient numerical solution of the problem, we should pass to dimensionless variables. To this end let us define the gas atomic mass as m 0 , and the atomic diameter as d 0 . We will denote the initial temperature and concentration of particles by T 0 and n 0 , respectively. From these values we can obtain the characteristic molecular velocity, free path length and time: With their help we can pass to dimensionless variables: In all calculations that follow the tilde will be omitted to avoid overloading the symbols. Using these new notations, we can set initial and boundary conditions. The initial conditions described when stating the problem are formalized as n 0 =1 inside the vessel and n 0 = 10 −6 in the rest of the computational region. Such concentration is set for efficient simulation of the absence of gas in that region without using the physically incorrect zero concentration value. The value of the distribution function is obtained assuming the Maxwellian distribution in the gas: Boundary conditions on the vessel walls are formalized as diffuse ones, which means that the mass flows coming to and going from the wall are equal, with the reflected molecules distributed according to the Maxwell law. If the wall is perpendicular to the Ox axis, this condition takes the form: When molecules get on the other boundaries of the computational region, they leave it forever, which means that the distribution function does not vary in the direction "from the computational region". If we consider, for example, the right-hand boundary, this condition will take the form: To speed up the calculations, it is reasonable to make use of the problem symmetry and perform calculations in the upper half of the computational region only. Reflective boundary conditions are set on the symmetry axis:

Solution algorithm
The direct solution method used consists of two stages -solution of the Boltzmann equation with the zero right-hand part, and calculation of the collision integral [23].
At the first stage, the Boltzmann equation reduces to a system of equations with parameter v.
Moreover, on the assumption that the distribution is close to the Maxwell one, values of v with the modulus greater than a certain cutting velocity v cut can be dropped. If the cutting velocity is on the order of 5 (in dimensionless units), the error in the case of the Maxwell distribution is less than 10 −6 . Now the Boltzmann equation can be solved numerically for each v such that |v| < v cut using standard methods. A discrete mesh is introduced with the time increment τ , the coordinate meshes h x , h y , and the velocity meshes ∆v x , ∆v y , ∆v z . Using the mesh, the gread function f k To calculate the value in the next temporal layer, there are used the difference secondorder TVD-scheme with the minmod limiter and the alternating direction method [24]. The boundary conditions are also written in terms of mesh functions, and are taken into account in each temporal layer. This method gives a nonnegative stable solution of the Boltzmann equation and allows going to calculation of the collision integral.
The second stage is implemented using a third-party library [25] that employs the two-point projection method described in [17]. This library changes the value of the distribution function in each node x i so that the influence of the collision integral is taken into account.
The solution algorithm described above is presented in more detail in the scheme below ( Figure 2).
Since the calculations of the distribution function values in each spatial point are performed independently, the algorithm can be naturally parallelized. For this purpose the computational region shall be broken down into sectors, with calculations in each of the sectors to be carried out by a separate processor core.
In this work this was implemented with the Nvidia CUDA, the technology wherein all computer cores of the video card have read and write access to the entire RAM. This adds convenience to solution of this problem, since it is not necessary to explicitly pass data at the sector boundaries from one core to another as would be the case with the MPI technology. Notice that the approached used in this work is more difficult to learn and requires profound knowledge of programming [21,22].
Based on the above-described algorithm, the library for calculation of the collision integral, and CUDA technology, there was developed an application for simulation of the gas flow in this problem. The application has a modular structure; it allows easily changing the problem geometry and can readily be used for other problems as well.

Simulation results
The simulation was carried out for four different configurations presented in Figure 3.
Density fields ( Figure 4) and pressure fields ( Figure 5) were obtained for each configuration for the quasi-steady flow.
Graphics of the Mach number versus distance at the symmetry line were also calculated. The graphics for t = 100 and t = 1600 are presented in Figure 7 and Figure 7, respectively.
Values of spread k, its corresponding angle and time derivative of the spread were calculated. The results are given in Table 1 below.
Thus it can be seen that the use of the GPU for solution of such problems is fully justified. Notice that in configuration (c) the distance to the screen is close to the distance at which the Mach number exceeds unity.      Figure 3 A comparative analysis of the speed of the application operation on the CPU and GPU with meshes of various sizes was carried out. The results are shown in Table 2.   Thus it can be seen that the use of the GPU for solution of such problems is fully justified.

Conclusions
The paper presents the study of gas outflow from a vessel in the case of a simple collimator. Various configurations are considered, such as without a collimator, with a collimator at short MRSU 2020 Journal of Physics: Conference Series 1560 (2020) 012062 IOP Publishing doi:10.1088/1742-6596/1560/1/012062 8 distance, with a collimator at long distance, and with a collimator with a smaller hole. The optimum collimator configuration was found to reduce the spread of the molecular flow.
The application was developed that allows simulating gas flows in different problems using the CUDA technology. It has been demonstrated that this tool is an efficient means of simulating flows of simple gases, including those with low Knudsen numbers. The use of the GPU allowed accelerating the execution of the code 3-15 times compared to the CPU.