The explanation of the DSMC method in Section 5.1 was compiled during my employment as a university educational material and a subcommittee material for academic societies, and similar content is also published in the Computational Mechanics Handbook III (Reference 3) compiled by the Japan Society of Mechanical Engineers. More than 20 years have passed since the first material was distributed, so it is not up-to-date, but it contains many of the points necessary to understand the content that appears on this website (blog). Although some parts have been revised from the original description, it is mainly based on the books (References 1, 2) and papers (References 12, 13, 14) by Bird, the founder of the DSMC method, so this section does not mention U_system or program parallelization.
Basics of the DSMC method
1. Introduction
The direct simulation Monte Carlo method (DSMC method) was developed by Bird in 1963 as a numerical calculation method for rarefied gas flows. The position coordinates and velocity components of many molecules are stored in the computer, and the “spatial movement” and “intermolecular collisions” of these molecules are separated in small time steps (uncoupling of molecular spatial movement and intermolecular collisions), and the flow field is analyzed by repeating the calculations alternately. This method became widely known around the world with the publication of Bird’s book on the DSMC method in 1976 (Reference 1), and a revised and significantly revised version was published in 1994 (Reference 2). During that time, molecular and collision models have made great advances, and coupled with advances in computer technology, it has become possible to use this method in a wide range of fields, from molecular flows to continuum flows, as well as in three-dimensional analysis. Furthermore, there have been attempts to apply this method, taking advantage of its microscopic analysis characteristics, to the analysis of unstable fluid phenomena (see Chapters 3 and 4 of this website). Here, in order to explain the basics of the DSMC method, we will take up a “supersonic rarefied free jet in an axisymmetric field” as an example and look at the details of the DSMC method while referring to its Fortran program (simplified version). Fortran has a very small number of instructions (only seven statements need to be understood: variable types, input/output statements, assignment statements, block IF statements, DO loops, GOTO statements, and subprograms; it is similar to VBA language, but calculations are fast because it is basically a compiler system) and even beginners can learn it in a relatively short time, so it is effective for understanding the DSMC method. In the DSMC method, many parts of the program can be used in common regardless of the problem being addressed, so we believe that the example problem in this section can be used as a core to apply it to various flow fields. It is also possible to develop this example problem from monatomic gases to diatomic or polyatomic gases, from simple gases to mixed gases, and even flows that include chemical reactions. We will look at the details in order, but first we will list the outline of the calculation procedure of the DSMC method and its features.
2. Overview of the calculation procedure of the DSMC method
[1] First, a molecule is given initial coordinates and initial velocity components by random numbers and placed in the target physical space. The physical space is divided into tiny spaces called cells. Ideally, the size of one side of a cell should be smaller than the local mean free path.
[2] The calculation of the movement of molecules in physical space (including collisions with solid walls and entry and exit at inflow and outflow boundaries) and the calculation of intermolecular collisions are separated by an infinitesimal time step
(a) Each molecule moves through space independently of other molecules for a time
(b) Since the position coordinates of each molecule change due to molecular movement, the cell to which the molecule belongs is investigated, and all molecules are rearranged in the “cross-reference array” in preparation for intermolecular collisions. In other words, each cell is given a consecutive number (seat number) equal to the number of molecules present in it as an array subscript, and the molecular number is stored in the array element. In collision calculations, once a seat number belonging to a cell is determined, this array is used to cross-reference the molecular number from the seat number, and information such as the molecular velocity components can be extracted.
(c) Intermolecular collisions are calculated for each cell based on the collision law. The main content is to calculate the number of collisions that occur within
[3]The above steps (a) to (c) are repeated, and the ensemble average of the flow field quantities is calculated for unsteady problems, and the time average after the steady state is achieved is calculated for steady problems. Note that ensemble averaging is performed to compensate for the small number of molecules per cell, but if you average cells in a range with little spatial change, or take a time average of a range with little temporal change, it is possible to perform an unsteady analysis without performing ensemble averaging. The animations shown in Sections 3.1 and 4.2 were obtained in this way.
3. Characteristics of the DSMC method (advantages and disadvantages)
[1] The calculation is based on the motion of individual molecules, so it has the disadvantage of taking a huge amount of time. DSMC is good at rarefied gases and high-speed gases, but when the density is high (even if high it is below 1 atmosphere) or when the speed is low, it takes a long time to reach a steady state, making the analysis more difficult. When the most probable molecular thermal speed is 350 m/s, the flow is analyzed based on that molecular motion, so in order to keep the statistical fluctuation of the results low, a significant number of sample molecules (or calculation repetitions) are required to clarify a phenomenon at a low speed, such as 1 m/s. However, if the analysis target is not the flow speed, it is possible to analyze even a stationary gas.
[2] Boundary conditions are easy to set. For example, it is only necessary to set the molecular reflection law at a solid wall, and in engineering calculations, simple diffuse reflection (scattered reflection) is often sufficient. The calculation of the inflow position and velocity components of molecules at the inflow boundary, or the number of inflow molecules, can also be easily determined according to the macroscopic state quantity. In the DSMC method, the conditions given to the inflow and outflow boundaries are applied only to molecules flowing in from the boundary, and are not directly applied to molecules flowing out of the calculation domain (of course, they are indirectly affected by collisions with molecules flowing in from the boundary). Therefore, even if an extreme state such as infinite distance is given to the boundary condition, the outflowing molecules are not strictly bound by it, and the calculation does not diverge. Of course, it is preferable for the boundary conditions to be such that there is no difference in the macroscopic state of the inflowing and outflowing molecules, so a method for achieving this is important. Note that if there is a difference between the two, distortion will occur in the flow field near the boundary, but this does not mean that the results for the entire domain will be completely meaningless.
[3] The DSMC method is essentially an unsteady calculation because the calculation progresses over time. Therefore, there is no essential difference in the difficulty of calculations between steady and unsteady problems. This means that even if it is a problem for solving a steady condition, a solution cannot be obtained unless an unsteady state has passed.
[4] Cell division (grid cutting) of the flow field is easy. Since the spatial movement of molecules is performed regardless of the cell structure, it is not necessarily necessary to divide the cells along the boundary shape (details will be described later).
[5] The main part of creating a new program is the routine for cell division and molecular movement. Existing routines can be used for the part related to intermolecular collisions. However, creating a program for molecular movement is extremely difficult without experience.
[6] It can be used to analyze not only rarefied gas flows but also micro flows. The DSMC method was developed as a numerical analysis method for rarefied gas flows. As shown in the example, it is now used in regions where the state changes significantly, from molecular flow to continuous flow, and in the future, it will be improved so that it can be applied to flows at atmospheric pressure. On the other hand, in the world of micro flows, the DSMC method can be used as is even at atmospheric pressure. The dimensionless number that represents the degree of rarefaction is the Knudsen number
4. Outline of the example
To explain the basics of the DSMC method, we will take up the example of “supersonic rarefied free jet in an axisymmetric field” and proceed with the explanation while referring to a Fortran program (Click here to see the program list). Note that the one to three digits and one space following each number on the left side of each line in the program indicate the line number (1 to 677) and are not included in the program. In the text, when describing variables in the Fortran program, we will underline the variable symbols to avoid misunderstandings (for example, TEMP). The example is a very simplified program of about 700 lines, but it has almost all the basic features. As shown in Figure 1, in a simple flow field consisting of a zero-thickness wall with a circular orifice of radius EZA (actual radius is 0.00125 m) and the space (computational space) that spreads to the right of it, the flow field is surrounded by inflow and outflow boundaries (boundaries with a density equivalent to the downstream back pressure, hereafter referred to as downstream boundaries) except for the orifice wall (circular aperture and wall surface). The downstream boundary can be divided into a boundary parallel to the central axis (peripheral boundary) and a boundary perpendicular to the central axis (vertical boundary), but when simply called downstream boundary, it refers to both.

Since the field is axisymmetric, the calculation space is represented only above the central axis (the central axis of the jet). The cell division is two-dimensional, but the flow field has a three-dimensional spread. The flow velocity is assumed to be uniform and equal to the speed of sound (
5. Generation of random numbers from a prescribed distribution
Here, we will introduce three methods for creating random numbers from a prescribed distribution required for DSMC calculations using real random numbers (hereafter referred to as uniform random numbers
[1] Inverse transform method (Direct method)
When the value
<Assignment> When molecules are uniformly distributed inside a spherical shell with an inner sphere of radius
[2] Acceptance-rejection method
The inverse transform method mentioned above cannot be used unless the equation connecting the cumulative distribution function

[3] A method to simultaneously obtain two normally distributed values
There is an efficient method to simultaneously obtain two normally distributed values, such as thermal velocity components in an equilibrium gas. Let
Now, if we replace the variables
Therefore,
However, because
That is, the velocity components
6. Function subprogram
In addition to the main program, the example program has three subprograms at the beginning. The first, RANF0(), is a function that generates uniform random numbers using the congruential method, and has a period of
The second subprogram, GAM(X), calculates the gamma function, and performs an approximation using Chebyshev’s series. This is used in the calculation of the VHS and VSS molecular models. The third subprogram, ERF(X), calculates the error function, and is used to calculate the number of inflowing molecules when there is a macroscopic flow rate at the inlet of the molecules. In the example, it is used to calculate the number of inflowing molecules from the orifice aperture. The formula for this error function uses the formula by Abramowitz et al.
7. Array variables used in DSMC calculations
The real variables used in the DSMC method are basically double-precision real numbers (see IMPLICIT statement). This is to minimize errors in the spatial movement of molecules. However, since the arrays for storing information on molecules or cells are large, single-precision real numbers are used for these large arrays in the example (if the PC has sufficient memory, it is better to use double-precision arrays as they are). In other words, the molecular velocity component P(3,MOL), the position coordinates PZ1(MOL) and PZ2(MOL), and the cell volume CTA(NCELL). Note that MOL is the maximum number of molecules, and NCELL is the maximum number of cells. By changing the values in the PARAMETER statement according to the capacity of the computer, it is possible to handle even larger numbers of molecules and cells. LCR stores the cross-reference array used in intermolecular collision calculations, and IP stores the cell number to which the molecule belongs. Note that since the example is an axisymmetric calculation, the
8. Starting the example (supersonic free jet) program
IERRx is a counter to check the number of errors during calculation, RVRALL and RVRGT are variables to check the maximum number of collisions in intermolecular collisions, IRAN is a random integer number (odd number), and PI is the circular constant.
The upstream stagnation point temperature TEMP and the orifice wall temperature TEMPW are 15[°C] (288[K]), the upstream stagnation point pressure PTOR is 50[Torr], PPAS is the pressure in [Pa], ANUM is the number density
In the DSMC method, molecular movement and intermolecular collisions are separated (uncoupled) at a time step
9. Cell division of the flow field
The cells that divide the flow field have two roles. One is that intermolecular collisions occur between molecules in the same cell, and the other is that the cells are used as infinitesimal elements to sample flow field information such as density, velocity, and temperature. Since molecules in the same cell are subjected to collision calculations regardless of their relative positions, ideally the length of one side of a cell should be smaller than the local mean free path. However, in axisymmetric or three-dimensional problems, it is difficult to make the cell length less than the mean free path when the size of the flow field and its density become large, and it is sometimes unavoidable to make the cells larger due to limitations in computer capacity. However, even if the cells are made larger, the results obtained do not become completely unreliable; the change rates of physical quantities such as density and temperature simply appear somewhat more relaxed than they should be. Therefore, cells can be made larger in places where the state changes are gradual, while if cells are made larger in places where the change is steep, results that are somewhat more relaxed than the original changes can be obtained. Keep this in mind when using the DSMC method. The simplest method of cell division is to divide the space into rectangular cells. Cell division can be done rectangularly even if the boundary is an arbitrary curved surface. This is because the cells are not involved in the spatial movement of molecules at all, and even in a rectangular cell structure, molecular movement is consistent with the boundary shape (the volume of the space in which the gas actually exists must be calculated). In the case of a complex boundary shape, some cells may become extremely small in volume, causing concern that there will be almost no molecules in the cell and causing problems in collision calculations. However, this effect will not be apparent unless a considerable number of nearby cells become small collectively.
Bird (Reference 12) devised the subcell division method in 1987 when comparing calculation methods with the molecular dynamics method (MD method) for vortex generation. The time counter method (TC method), which was the mainstream intermolecular collision methods at the time, required more than 20 molecules to be contained in a cell, so in order to ensure this number and keep the distance between colliding molecules below the mean free path, the normal cell grid was divided into even finer cells. However, in the currently commonly used maximum collision number method (null collision method, no-time-counter method), the theoretical collision number is satisfied even if the number of molecules in the cell is about one (Reference 15), so the role of the subcell method is not as important as it was originally. In addition, the program code of the subcell method shown in Bird’s book (Reference 2) was created for one-dimensional calculations, so it does not necessarily provide accurate calculations when applied to two- or three-dimensional problems. Recently, Usami et al. (References 3, 4, 5) developed a method called the U-system to calculate collisions taking into account the molecular positions within the cell. The details of this method, which gives very good results even if the cells are much larger than the mean free path, are explained in Section 5.2.
The density of the free jet in the example drops rapidly as it moves downstream, so the cell structure (axial cell division) that subdivides the flow field is designed to increase the cell length in a geometric progression as it moves downstream, so that there is no large difference in the number of molecules in each cell. Now, on the central axis, the axial length of the cell immediately after the orifice is
In this case, the number of cells in the axial direction KS is 164, and the common ratio ES0 is 1.02. On the other hand, the radial cell division is basically first divided equally into 40 cells (I2) of length TKN for the orifice radius EZA, and the entire radial cell structure of the flow field is also divided equally, with the number of divisions KO being 240. AL1 and AL2 in lines 82-83 are
If the radial direction is divided equally in an axisymmetric field, the cell volume on the central axis becomes extremely small and a sufficient number of molecules cannot be secured. Therefore, in the example, only the cells on the central axis are operated to combine three adjacent cells in the radial direction to form one cell.
The cell volume CTA is calculated in the double DO loop in lines 87-104. Because of the axisymmetric flow field, the cells other than those on the central axis have a ring structure. The three cells closest to the central axis are combined to form one, making this a somewhat complicated routine. Additionally, the total number of cells, NC, is KS
<Assignment> A plate with zero thickness is placed diagonally in a flow field divided into a grid of rectangular cells, and some cells are divided into the front and back of the plate. In such a case, consider how to handle the cells divided by the plate.
10. Number of molecules in a cell and weighting factor of simulated molecules
The number of molecules in a gas is
In this way, regardless of the number of molecules in the real gas, the number of molecules in the DSMC method can be determined quite freely, but it is necessary to remember how many real molecules each simulation molecule represents. This is called the weighting factor WEI. In other words, in calculating the number of intermolecular collisions per time step
In the example, the real molecular number density at the upstream stagnation point is ANUM [1/m
<Assignment> If the number of molecules in the cell is one in time average, a question arises as to whether a collision partner cannot be found. However, the simulation is valid. Consider this point.
11. Molecular model (VHS molecule and VSS molecule) – collision calculation part 1
One of the significant improvements in the DSMC method is the molecular model. In the early days (Reference 1), hard sphere molecules were used if calculation speed was a priority, and inverse power law molecules were used if faithfulness to real molecules was a priority. Later, Bird analyzed many calculation results and concluded that “the reason inverse power law molecules produce good results is not because they faithfully reproduce the deflection angle in collision calculations, but because the collision cross section is a function of the relative translational energy.” In a 1981 paper (Reference 13), he published the VHS molecular model, which maintains the hard sphere scattering law (in the collision of hard sphere molecules, there is no probabilistically dominant direction for the relative velocity after the collision) and makes the collision cross section a function of the relative velocity (i.e., a molecular model in which the molecular diameter is a function of the relative velocity). This molecular model has excellent calculation speed because it has collision dynamics almost the same as hard spheres, and it also produces calculation results similar to the inverse power law molecule. Now, assuming that the temperature
Here,
Table 1 Temperature exponent of viscosity coefficient ω

Here,
Furthermore, using the diameter
The Knudsen number
That is,
The VHS molecule was subsequently improved by Koura and Matsumoto (Reference 17) into the VSS molecular model, which also takes into account the diffusion coefficient.
In the VSS molecule, the collision cross section is considered to be the same as in the VHS molecule, but the deflection angle
Here,
Table 2 α for various gases

The reference diameter
Now, as in equation (5) obtained in the explanation of the inverse transform method,
Of course, when
For a single gas, there is almost no difference in the results between VHS and VSS molecules, but it has been reported that VSS molecules are more effective when calculating mixed gases. However, VHS molecules are more advantageous in terms of calculation time.
In the example problem, we deal with a single gas, argon, and whether to use VHS or VSS molecules is determined by substituting 0 or 1 for IVSS (lines 113-123). AMOL is the molecular weight, RGAS is the gas constant, AMAS is the mass of one molecule, BOL is the Boltzmann constant, AMR is the reduced mass, HNH is the specific heat ratio, VISC is the viscosity coefficient, DIRM is the reference molecular diameter, ACXS is the reference collision cross section, VM is the most probable molecular thermal speed (most probable velocity), BETA is
12. Calculation of intermolecular collisions – Collision calculation part 2
The most discussed issue in the development of the DSMC method is the model of intermolecular collisions. To calculate the number of collisions per
Now, consider that there are
Therefore, the usual method is to check whether or not a collision is possible for all combinations of
In other words, the selection of molecular pairs is repeated the number of times shown in equation (34), and the probability of the molecular pairs actually colliding is judged by
In this way, a method for efficiently calculating collisions without trying all combinations has been completed. Note that this method requires that
In the example, the maximum relative speed VRMAX is set to three times the average relative speed (CMUL
Another notable intermolecular collision calculation method is the Nambu method (Nambu method, modified Nambu method) (Reference 20). This method is highly regarded as a method that faithfully incorporates the Boltzmann equation into the collision calculation. However, since only the velocity of one of the molecules changes during the collision, the calculation speed is somewhat slow, and there are reports that the statistical fluctuations in the results are greater than in other methods, so it is used less frequently than the maximum collision number method.
For molecular pairs that are determined to be collidable according to the probability of equation (35), the components of the relative velocity after the collision
Here, the subscripts 1 and 2 distinguish the two molecules. Finally, the velocity components
Diatomic or polyatomic molecules store energy not only through translational motion (translational energy mode), but also through rotational motion around the center of mass (rotational energy mode). In addition, at high temperatures, vibrational motion is also excited (vibrational energy mode). At present, the Larsen-Borgnakke phenomenological molecular model (LB model) (Reference 21) is often used to handle this with the DSMC method. This assumes the conservation of total energy and calculates the value of each energy mode based on the equilibrium energy distribution. However, in order to use this phenomenological model, it is necessary to determine in advance, for example, the number of collisions (or the ratio of elastic and inelastic collisions) required for the translational and rotational modes to reach equilibrium.
13. Inflow of molecules from the inlet boundary
For a stationary gas in equilibrium, the number of molecules incident on a unit area per unit time is
Here,
The number of molecules flowing in from the boundary is calculated from these equations. Next, the velocity components of the inflowing molecules are considered to be separated into components parallel to the boundary and components perpendicular to it. The velocity component parallel to the boundary is first calculated using the same equations (8), (10), and (12) as those for the velocity generation of molecules in equilibrium, and then the macroscopic flow velocity in that direction is added to it to complete the equation. On the other hand, if the molecular velocity component perpendicular to the boundary is
The acceptance-rejection method is used to generate the molecular velocity from this distribution. In the acceptance-rejection method, it is first necessary to determine the “range of
By substituting this for
[1] Find the uniform random number
[2] Substitute the above
[3] Find the uniform random number
Note that
Now, when
Considering that the range of
As in equation (11),
Therefore,
The vertical velocity component of the molecules diffusely reflected from the wall can be calculated in the same way using equation (46). However,
As mentioned earlier, the state at the orifice entrance in the example problem is calculated assuming that the gas reaches the orifice entrance from the stagnation point with an isentropic change, where the flow velocity becomes the speed of sound
DTM is the time step
FIN is the number of simulated molecules flowing in from the orifice inlet, DFIN is the number of simulated molecules flowing in from the downstream vertical boundary perpendicular to the central axis, and DFINS is the number of simulated molecules flowing in from the downstream peripheral boundary parallel to the axis. NM is the total number of simulated molecules, and its initial value is the number of molecules initially placed in the flow field. In addition, the maximum number of molecules is stored in NMM.
Lines 196 to 201 initialize (zero) the array variables for “number of passing molecules” and “cell sampling information”. In line 199,
In lines 205-207, it is determined whether to “start simulation” or “process data”, and branches to the respective routine.
In lines 209-221, the position coordinates
DO 100 in line 224 is the start of the main iterative routine (four routines for molecular movement, etc., and cell information sampling).
Lines 225-229 output error information every 10 RUNs. DO 24 in line 232 is the start of a loop that includes four routines: molecular movement, survei of belonging cell for molecules, rearrangement of molecules in cross-reference array, and intermolecular collision. However, note that there is only one loop (NIS=1) in the present program.
In line 235, IN is
Line 239 starts the molecular movement in physical space, but we will postpone the explanation to the next section 14. First, we will explain the molecular inflow program from the inflow boundary (lines 307-371).
Lines 309-327 are the molecular inflow from the orifice aperture. Lines 309-311 convert the real number of inflowing molecules into an integer using a random number. SS(5) in line 312 is the cumulative of FINBAS, where FINBAS (line 190) is the number of gas molecules at the stagnation point density entering the same cross-sectional area as the orifice aperture in time
Lines 336 to 349 are the inflow of molecules from the downstream peripheral boundary parallel to the central axis, and lines 356 to 371 are the inflow of molecules from the downstream vertical boundary perpendicular to the axis. The calculation process is almost the same except for whether the velocity component perpendicular to the boundary is in the
<Assignment> Perform calculations using a program that removes the addition of U0(KU) in line 367, and compare the magnitude of distortion near the vertical boundary.
14. Movement of molecules in physical space (flow field)
Each molecule moves in physical space independently of other molecules for a time of
As described above, when creating a new program by the DSMC method, most of the effort is usually devoted to creating this molecular movement routine.
Collisions between molecules and walls are generally found by simultaneously solving the molecular trajectory equation and the wall equation. For example, assume that there is a wall in two-dimensional space that can be expressed by the following equation
If the molecular velocity is
Therefore, the time
If
Note that when the wall equation is
Next, we will introduce an equation to find the collision time for a simple wall equation.
[1] When a molecule moves inside a circle (
Of course, if
[2]When a molecule moves outside the circle (
In this case, if
[3]When a molecule moves in the double circle space created by the above two circles, first check for a collision with the inner circle. If there is no collision with the inner circle, check for a collision with the outer circle. Note that there is a possibility of collision with the outer circle more than twice in a row.
[4]When a molecule moves inside a sphere (
[5]When a molecule moves outside a sphere (
[6]When a molecule moves in the space inside the double sphere formed by the two spheres in [4][5] above, first check for collision with the inner sphere. If there is no collision with the inner sphere, check for collision with the outer sphere.
<Assignment> Find the collision time when a molecule moves inside or outside the surface of a cone.
When a molecule moves in an axisymmetric space, it is necessary to perform a coordinate rotation operation after the movement. The molecule moves in space according to three velocity components, but since only two position coordinates are held, a process is required to constantly set the remaining coordinate to zero. Now, consider the
When a molecule collides with a wall during its movement, diffuse reflection (scattered reflection) is calculated for normal walls (engineering walls). Diffuse reflection is when all the momentum and energy that the molecule had at the time of incidence is absorbed by the wall, and the molecule is assumed to be completely accommodated to the wall, and the new reflection velocity is calculated from the equilibrium distribution at the wall temperature. In other words, first, the time until the molecule collides with the wall is calculated, and this is subtracted from the total movement time to calculate the flight time after reflection. Next, calculate the velocity components after reflection as follows, and let the molecule fly at this new velocity.
Now, assume that in two-dimensional space, the normal vector of the wall is tilted at an angle of

First, we calculate the velocity component
Here,
Next, we convert these velocities to velocity components
Another type of wall reflection is specular reflection. The angle of incidence and the angle of reflection are equal, and the molecular velocity does not change except for the reversal of the velocity component perpendicular to the wall, and there is no exchange of energy between the molecule and the wall. In general, such a surface does not exist, but this reflection occurs when a symmetrical surface is considered. There is also a method of simulating an intermediate reflection between the two by predetermining the ratio of diffuse reflection and specular reflection (accommodation coefficient). This method can be used in two ways: one is to determine whether each reflection is diffuse reflection or specular reflection using a random number, and the other is to combine the velocity component in each reflection from both diffuse reflection and specular reflection. The former is more common, but the latter method has the advantage that a lobular reflectance distribution can be easily achieved even with a single incident angle.
Ideally, the starting point of the reflected molecule should be a point on the wall, but in that case, calculation errors may cause the molecule to start from inside the wall, which can have a negative effect on the simulation. To prevent this, it is better to fly the molecule from a point slightly above the wall (a point a little into the flow field). How far the molecule should be above the wall depends on the shape of the wall, but about
The routine for molecular movement starts from Line 239. After advancing the molecule number by one, the position coordinates and velocity components stored in the single precision array are temporarily replaced with double precision variables to reduce calculation errors. Now, for molecules that have been in the flow field from the beginning, the time step
Lines 272 to 279 process molecules that do not collide with the orifice wall or molecules that depart from the wall. In line 272, the outflow counter is advanced for molecules that have crossed the
Here, we will explain the molecule discarding process. Since the information of molecules that have jumped out of the calculation area is not needed, they will be discarded from the array that stores the information, but the molecular information is stored randomly in the array. Therefore, simply discarding them will leave the array in a worm-eaten state. In other words, in order to store the information of new molecules that will flow in, it is necessary to organize the array. Therefore, as shown in lines 284 to 290, the information stored at the end of the array at that time is transferred to the element that stored the molecular information to be discarded, and the number at the end of the array (total number of molecules) is decremented by one.
In lines 300 to 303, after all the calculations for the movement of the molecules that were in the flow field from the beginning are completed, the process moves on to the process of starting the calculation of the movement of the molecules flowing in from the boundary. Once that calculation is also completed, check the update of the maximum number of molecules, and then move on to the routine of “calculation to determine which cell the molecule is in”.
<Assignment> Draw a flow chart for the “spatial movement routine of molecules” in the example program and track the movement of the molecules.
15. Calculating the cell to which the molecule belongs and sorting the molecules on the cross-reference array
Before proceeding with the calculation of intermolecular collisions, a procedure is required as a preparation. Since the position coordinates change due to the spatial movement of the molecules, first survey which cell each molecule belongs to. In the example, in lines 392 to 394, the sequence number in the
On the other hand, in the
Now, to make it easier to calculate intermolecular collisions, we need to process the rearrangement of the molecules in the cross-reference array LCR(K) in the order of the cells to which they belong. Figure 4 shows an example of the array created by this process, where K represents the seat number (address), and the molecule number (molecule identification number) is stored in LCR(K). In Figure 4, the first cell has 6 molecules, reserves 6 seats from 1 to 6, and the seat number just before the first seat (starting_address

To create this cross-reference array, first, the number of molecules in the cell is investigated for all cells. That is, in lines 414 to 418, the number of molecules in the cell, IC1(N), is set to zero, and then the number of molecules in the cell is counted while checking the cells to which all molecules belong. Next, in lines 419 to 423, the number of molecules in the cell is accumulated in the order of cells, and the first seat of each cell on the LCR(K) array is investigated. In fact, the seat number just before the first seat of each cell, starting_address
16. Program for intermolecular collisions – collision calculation part 3
The calculation of intermolecular collisions is performed for each cell in a DO loop (lines 433 to 509) that spans all cells. In line 435, if there are not two or more molecules in a cell, the collision calculation for that cell is skipped. Lines 436 and 437 calculate the number of trials (number of collision attempts) for selecting a colliding molecule pair using equation (34). As mentioned earlier, the collision cross section
Lines 444 to 447 are the calculation to randomly select one of the molecules in the cell using the number of molecules in the cell IC1(N), the starting_address
Lines 468-473 calculate the relative velocity components of VHS molecules after collision using equations (20)-(24), and lines 475-492 calculate the relative velocity components of VSS molecules after collision using equations (26), (28)-(32). Line 495 counts the number of collisions and is not used directly in this example. Lines 497-505 calculate the velocity components of both molecules after collision using equations (36)-(38). Here, VCCM is the center of mass velocity, but since the molecules have the same mass,
17. Sampling of flow field data
The above-mentioned molecular movement and intermolecular collisions are repeated endlessly to form a flow field, but the only things that change in the computer are the position coordinates and velocity components of the many molecules. The information we actually want from the flow field is density (number density), flow velocity, and temperature, so we need to calculate them. The number density can be calculated from the number of molecules in each cell,
When calculating pressure, if the equation of state holds, it is calculated from the density and temperature. If the equation of state cannot be used due to a strong non-equilibrium condition, a virtual wall is placed at the position where the pressure is to be calculated (in the program) and the pressure is calculated from the momentum of the molecules passing through it.
In the example, the DO loop starting from line 514 is a sampling process for the information in the cell. For each molecule, the cell to which it belongs is investigated, and in line 517 the number of molecules is counted using equation (67), in lines 518 to 520
Intermediate summaries of the in-cell sampling data are output to a file every K100 times.
Lines 529 and 530 determine whether the number of times for data output has been reached. If the number of times has been reached, lines 531 to 591 are executed, and if the number of times has not been reached yet, the program returns to line 224.
The file name output every K100 times is different each time. That is, they are named ‘sjetc01’, ‘sjetc02’, …. Lines 531-542 determine the name, and lines 544-553 output the cell information to a file.
The other output file is the net number of passing molecules at the orifice inlet and downstream boundary. In this example, it is non-dimensionalized by the number of incident molecules at the stagnation point
18. Final Data Processing
In the example program, the final data processing routine is created separately from the main routine. In other words, once the simulation is finished, the number of output data to average is determined as described above, and the program is executed again. Entering ‘2’ in lines 205 to 207 will move to the data processing routine. XLIM in line 601 inputs how many times the orifice diameter the flow field length should be processed in the
Lines 636 to 671 are the routines that calculate the final flow field data and output it to the file ‘OUTD’. The control variable I in DO 302 changes in the order of the cells arranged in the
The output data is arranged in the following order for each cell:

Figure 5 shows the density distribution calculated in the example, processed by gnuplot (v3.5). It is expanded in the height direction, so the Mach disk and barrel shock can be clearly observed. However, in the example, the calculation area is not set upstream of the orifice, so the flow rate is larger than it actually is, and the jet itself is larger, and the Mach disk is shifted downstream from its original position.
<Assignment> When considering the effect of back pressure on jet formation, depending on the flow field conditions, the effect of the downstream peripheral boundary parallel to the jet axis may be dominant, and the downstream vertical boundary perpendicular to the jet axis may have almost no effect. In the example program, perform a calculation by setting only the downstream vertical boundary as a complete vacuum, and compare the results.
19. Development from the example program
The trick to making a good program is to create many programs yourself. We hope that you will use this example program as a basis to challenge yourself to more complex problems. Finally, some development tasks are added.
[1] In the example, the free jet is large because there is no computational domain upstream of the orifice, and the shape of the jet immediately after the orifice exit is slightly different from the actual shape. Therefore, create a program with a computational domain upstream and perform a DSMC calculation. In that case, how should the boundary conditions at the upstream boundary be set?
[2] Give the orifice a thickness. Also, how can you make it a conical nozzle?
[3] In the example, the cell length is changed depending on the location. If the time step could also be changed depending on the location, efficient calculations would be possible in situations with a large pressure ratio. How can this be achieved? (Difficult problem)
[4] The example shows an axisymmetric flow field, but if the flow field could be made three-dimensional, a free jet passing through a rectangular orifice could also be calculated. How should the example program be modified?
[5] What part of the example program should be modified to realize a mixed gas flow?
To see the example program, please click here. (Please note that the one to three digits and the space following them on the left side of each line in the program indicate the line number (1 to 677) and are not included in the program itself.)