Chapter 2 Numerical Methods
2.1 Direct Simulation Monte Carlo
2.1.2 The DSMC Procedures
Fig. 2.1 is a general flowchart of the DSMC method. Important steps of the DSMC method include setting up the initial conditions, moving all the simulated particles, indexing (or sorting) all the particles, colliding between particles and sampling the molecules within cells to determine the macroscopic quantities. The details of each step will be described in the following subsections.
Initialization
The first step to use the DSMC method in simulating flows is to set up the geometry and flow conditions. A physical space is discredited into a network of cells and the domain boundaries have to be assigned according to the flow conditions. A point has to be noted is the cell dimension should be smaller than the mean free path, and the distance of the molecular movement per time step should be smaller than the cell dimension. After the data of geometry and flow conditions have been read in the code, the numbers of each cell is calculated according to the free-stream number density and the current cell volume. The initial particle velocities are assigned to each particle based on the Maxwell-Boltzmann distribution according to the free-stream velocities and temperature, and the positions of each particle are randomly allocated within the cells.
Molecular Movement
After initialization process, the molecules begin move one by one, and the molecules move in a straight line over the time step. A special particle ray-tracing technique has to be designed to efficiently track the particle movement for the special grid system, unstructured grid, which we use in the current study. The particle ray-tracing technique of three-dimensional domain is described in the following, respectively.
3-D Particle Ray-Tracing in Unstructured Mesh
Fig 2.2 is the sketch of the particle movement in three-dimensional unstructured mesh.
The details are described in the following:
Without considering the external force effects, free-flying position of traced particle at t
On the other hand, cell face can simply be represented as a planar equation as
= 0
)
By computing in turn Eq. (2.3) of each face in the current cell, the correct intersecting face number can be identified by finding the minimum positive ∆ 't. The intersecting coordinate can be found by substituting ∆ 't into Eq. (2.1). If the intersecting face is a normal face between cells, then the particle will continue its trajectory until it stops.
If the intersecting face is a solid face, the particle will be reflected in a special way (e.g., diffusively or secularly) according to the specified boundary condition. These are related by the coordinate transformation matrix between the local coordinate system (on the face) and
the absolute coordinate system for both types of conditions. First, a unit vector x' along
the face is chosen, then y' is the cross product of x' and z' (the normal unit vector of the
The coordination transformation matrix H , is
⎥⎥
transformation H−1 can be easily written as
−1
H =HT (2.6)
where HT is the transpose matrix of H .
Now, the particle velocity can be transformed from the velocity in absolute coordinate system (Vabs) to the velocity in local coordinate system (Vloc) before the reflection by using
H .
Vloc=
H V
abs (2.7) After the reflection of the particle, the new local coordinate system velocity (Vloc') can be written as=
'
Vloc F(Vloc, wall condition) (2.8)
where F(Vloc, wall condition) is a kernel function, depending upon the wall condition and velocity before reflection.
Finally, the absolute velocity after the reflection (Vabs') will be obtained by using the inverse transformer as
Then, the particle continues its journey with its new absolute velocity until it stops.
Indexing
The location of the particle after movement with respect to the cell is important information for particle collisions. The relations between particles and cells are reordered according to the order of the number of particles and cells. Before the collision process, the collision partner will be chosen by a random method in the current cell. And the number of the collision partner can be easy determined according to this numbering system.
Gas Phase Collisions
The other most important phase of the DSMC method is gas phase collision. The current DSMC method uses the no time counter (NTC) method to determine the correct collision rate in the collision cells. The number of collision pairs within a cell of volume Vc over a time interval ∆t is calculated by the following equation;
c
N and N are fluctuating and average number of simulated particles, respectively.
F is the particle weight, which is the number of real particles that a simulated particle N
represents. σ and T cr are the cross section and the relative speed, respectively. The collision for each pair is computed with probability
)max
/(
)
(σTcr σTcr (2.11)
The collision is accepted if the above value for the pair is greater than a random fraction.
Each cell is treated independently and the collision partners for interactions are chosen at random, regardless of their positions within the cells. The collision process is described sequentially as follows:
1. The number of collision pairs is calculated according to the NTC method, Eq. (2.10), for each cell.
2. The first particle is chosen randomly from the list of particles within a collision cell.
3. The other collision partner is also chosen at random within the same cell.
4. The collision is accepted if the computed probability, Eq. (2.11), is greater than a random number.
5. If the collision pair is accepted then the post-collision velocities are calculated using the mechanics of elastic collision. If the collision pair is not to collide, continue choosing the next collision pair.
6. If the collision pair is polyatomic gas, the translational and internal energy can be redistributed by the Larsen-Borgnakke model [1975], which assumes in equilibrium.
The collision process will be finished until all the collision pairs are handled for all cells and then progress to the next step.
Sampling
After the particle movement and collision process finish, the particle has updated positions and velocities. The macroscopic flow properties in each cell are assumed to be constant over the cell volume and are sampled from the microscopic properties of each particle within the cell. The macroscopic properties, including density, velocities and temperatures, are calculated in the following equations [Bird, 1976 and Bird, 1994];
=nm velocity, mean velocity, and random velocity, respectively. In addition, Ttr, Trot, Tv and Ttot
are translational, rotational, vibrational and total temperature, respectively. εrot and εv are
the rotational and vibrational energy, respectively. ζrot and ζv are the number of degree of freedom of rotation and vibration, respectively. If the simulated particle is monatomic gas, the translational temperature is regarded simply as the total temperature. Vibrational effect can be neglect if the temperature of the flow is low enough.
The flow will be monitored if steady state is reached. If the flow is under unsteady situation, the sampling of the properties should be reset until the flow reaches steady state.
As a rule of thumb, the sampling of particles starts when the number of molecules in the calculation domain becomes approximately constant.