2. 1 The Boltzmann Equation
The Knudsen number (Kn) is used to indicate the degree of rarefaction. In Figure 2.1, flows are divided into four regimes and three solutions. We have found the Boltzmann equation is valid for all flow regimes which from 10 to 0.0001. It is one of the most important transport equations in non-equilibrium statistical mechanics, which deals with systems far from thermodynamics equilibrium. There are some assumptions made in the derivation of the Boltzmann equation which defines limits of applicability. They are summarized as follows:
1. Molecular chaos is assumed which is valid when the intermolecular forces are short range. It allows the representation of the two particles distribution function as a product of the two single particle distribution functions.
2. Distribution functions do not change before particle collision. This implies that the encounter is of short time duration in comparison to the mean free collision time.
3. All collisions are binary collisions.
4. Particles are uninfluenced by intermolecular potentials external to an interaction.
According to these assumptions, the Boltzmann equation is derived and shown as Eq. (2.1)
molecular velocity.
σ
is the differential cross section and dΩ is an element of solid angle.The prime denotes the post-collision quantities and the subscript 1 denotes the collision partner. Meaning of each term in Eq. (2.1) is described in the following;
1. The first term on the left hand side of the equation represents the time variation of the distribution function of the particles (unsteady term).
2. The second term gives the spatial variation of the distribution function (flux term).
3. The third term describes the effect of a force on the particles (force term).
4. The term at right hand side of the equation is called the collision integral (collision term). It is the source of most of the difficulties in obtaining solutions of the Boltzmann equation.
In general, it is difficult to solve the Boltzmann equation directly using numerical method because the difficulties of correctly modeling the integral collision term. Instead, the DSMC method was used to simulated problems involving rarefied gas dynamics, which is the main topic in the current thesis.
2. 2 General Description of the Standard DSMC
In order to the expected rarefaction caused by the rarefied gas flows, the direct simulation Monte Carlo (DSMC) method which is a particle-based method developed by Bird during the 1960s and is widely used an efficient technique to simulate rarefied gas regime [Bird, 1976 and Bird, 1994]. In the DSMC method, a large number of particles are generated in the flow field to represent real physical molecules rather than a mathematical foundation and it has been proved that the DSMC method is equivalent to solving the Boltzmann equation [Nanbu, 1986 and Wagner, 1992]. The assumptions of molecular chaos and a dilute gas are required by both the Boltzmann formulation and the DSMC method [Bird, 1976 and Bird, 1994]. An important feature of DSMC is that the molecular motion and the intermolecular collisions are uncoupled over the time intervals that are much smaller than the mean collision time. Both the collision between molecules and the interaction between molecules and solid boundaries are computed on a probabilistic basis and, hence, this method makes extensive use of random numbers. In most practical applications, the number of
simulated molecules is extremely small compared with the number of real molecules. The general procedures of the DSMC method are described in the next section, and the consequences of the computational approximations can be found in Bird [Bird, 1976 and Bird, 1994].
In DSMC, there are three molecular collision models for real physical behavior and imitate the real particle collision, which are the Hard Sphere (HS), Variable Hard Sphere (VHS) and Variable Soft Sphere (VSS) molecular models, in the standard DSMC method [Bird, 1994]. The collision pairs then are chosen by the acceptance-rejection method. The no time counter (NTC) method is an efficient method for molecular collision. This method yield the exact collision rate in both simple gases and gas mixtures, and under either equilibrium or non-equilibrium conditions.
Figure 2.2 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;
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. An important feature has to be noted is the size of the computational cell 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 if it did not collide with solid surface. For the standard DSMC code by Bird [Bird, 1976 and Bird, 1994], the particles are moved in a structured mesh. There are two possible conditions of the particle movement. First is the particle movement without interacting with solid wall. The particle location can be easy located according to the velocity and initial locations of the particle. Second is the case that the particle collides with solid boundary. The velocity of the particle is determined by the boundary type. Then, the particle continues its journey from the intersection point on the cell surface with its new absolute velocity until it stops. Although it is easier to implement by using structured mesh, it is difficult for those flows with complex geometry.
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
V
C over a time interval ∆ is calculated by the following equation;t
c r
T
N
c t V
F N
N
( ) /12
max∆
σ
(2.2)N
andN
are fluctuating and average number of simulated particles, respectively.F
Nis the particle weight, which is the number of real particles that a simulated particle represents.
σ
T andc
r are the cross section and the relative speed, respectively. The collision for each pair is computed with probability)max
/(
)
(
σ
Tc
rσ
Tc
r (2.3)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.2), 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.3), 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 [Borgnakke and Larsen, 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, vibration and total temperature, respectively.ε
rotandε
vare the rotational and vibration 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 total temperature. Vibration effect can beneglect 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.
2. 3 General Description of the PDSC
Although the large number of molecules in a real gas is replaced with a reduced number of model particles, there are still a large number of particles must be simulated, leading to tremendous computer power requirements and needing to cost a lot of computational time. As a result, parallel DSMC method is developed to solve the problem. Figure 2.3 illustrates a simplified flow chart of the parallel DSMC method used in the current study. The DSMC algorithm is readily parallelized through physical domain decomposition. The cells of the computational grid are distributed among the processors. Each processor executes the DSMC algorithm in serial for all particles and cells in its domain. Data communication occurs when particles cross the domain (processor) boundaries and are then transferred between processors.
Parallel DSMC Code (PDSC) is the main solver used in this thesis, which utilizes unstructured tetrahedral mesh. Figure 2.4 is the features of PDSC and brief introduction is listed in the following paragraphs.
1.
2D/2D-axisymmetric/3-D unstructured-grid topology: PDSC can accept either
2D/2D-axisymmetric (triangular, quadrilateral or hybrid triangular-quadrilateral) or 3D (tetrahedral, hexahedral or hybrid tetrahedral-hexahedral) mesh [Wu et al.’s JCP paper, 2006]. Computational cost of particle tracking for the unstructured mesh is generally higher than that for the structured mesh. However, the use of the unstructured mesh, which provides excellent flexibility of handling boundary conditions with complicated geometry and of parallel computing using dynamic domain decomposition based on load balancing, is highly justified.2.
Parallel computing using dynamic domain decomposition: Load balancing of PDSC
is achieved by repeatedly repartitioning the computational domain using a multi-level graph-partitioning tool, PMETIS [Wu and Tseng, 2005] by takingadvantage of the unstructured mesh topology employed in the code. A decision policy for repartition with a concept of Stop-At-Rise (SAR) [Wu and Tseng, 2005]
or constant period of time (fixed number of time steps) can be used to decide when to repartition the domain. Capability of repartitioning of the domain at constant or variable time interval is also provided in PDSC. Resulting parallel performance is excellent if the problem size is comparably large. Details can be found in Wu and Tseng [Wu and Tseng, 2005].
3.
Spatial variable time-step scheme: PDSC employs a spatial variable time-step
scheme (or equivalently a variable cell-weighting scheme), based on particle flux (mass, momentum, energy) conservation when particles pass interface between cells.This strategy can greatly reduce both the number of iterations towards the steady state, and the required number of simulated particles for an acceptable statistical uncertainty. Past experience shows this scheme is very effective when coupled with an adaptive mesh refinement technique [Wu et al.’s CPC paper, 2004].
4.
Unsteady flow simulation: An unsteady sampling routine is implemented in PDSC,
allowing the simulation of time-dependent flow problems in the near continuum range [JCP paper submitted in June 2007]. A post-processing procedure called DSMC Rapid Ensemble Averaging Method (DREAM) is developed to improve the statistical scatter in the results while minimizing both memory and simulation time.In addition, a temporal variable time-step (TVTS) scheme is also developed to speed up the unsteady flow simulation using PDSC. More details can be found in [JCP paper submitted in June 2007]. Details of the idea and implementation are described next.
5.
Transient Sub-cells: Recently, transient sub-cells are implemented in PDSC directly
on the unstructured grid, in which the nearest-neighbor collision can be enforced, whilst maintaining minimal computational overhead [JFM paper under preparation, 2007].2. 4 General Description of Unsteady Sampling Method in DSMC [JCP paper in March 2008]
In section 2.3, the PDSC code has been specifically designed for simulating steady flows,
therefore some modification is need for unsteady sampling. The unsteady sampling method has been described in detail in the paper [Cave, et al., 2008].
There are two methods for unsteady sampling, the differences illustrated in Figure 2.5 and the details will be described in the following;
7. The “Ensemble-Average” method, require multiple simulation runs. The flow flied is sampled at the suitable sampling times during the run. The sampling simulation outputs from each run are averaged over the runs. There the results are vary precise, but the method is very computational expensive. Because a large of runs is required to reduce the statistical scatter to smooth data and a large amount of memory is needed to record the sampling data for each simulation.
8. The “Time-Average” method, require one simulation run. It averages a number of time steps over an interval before the sampling time. However it suffers a potential disadvantage in that the results will be “smeared” over the time over which samples are taken. Hence the sample time must be sufficiently short to minimizes time
“smearing” and yet long enough to obtain a good statistical sample. This method of time averaging has been used previously by Auld to model shock tube flow [Auld, 1992]
In PDSC, the method of time-averaging was implemented [Cave, et al., 2008]. Here a technique called the temporal variable time step (TVTS) method was used to reduce the simulation time by increasing the time step between sampling. The code has an option for the user to choose specific output flow times or for output at regular intervals. Figure 2.6 shows the flow chart of the PDSC method with the unsteady sampling procedures implemented.
Here M is the output matrix for sampling interval M. Most parts of the procedure are the same as the steady simulation except the sampling data must be reset after completing each
simulation interval.
2. 5 DSMC Rapid Ensemble Averaging Method (DREAM) [JCP paper March 2008]
Because reducing the statistical scatter greatly, in time-average data requires a very large number of simulation particles with consequent large computational times. In the thesis, we have adapted DREAM code which has been described in detail in this paper by [Cave, et al., 2008]. The illustrated in Figure 2.7
First, we select a raw data set X-n produced by PDSC n sampling intervals prior to the sampling interval of interest X. New particle data is generated from the macroscopic properties in data set X-n by assuming a Maxwellian distribution of velocities. The standard PDSC algorithm is then used to simulate forward in time until the sampling period of interest
X is reached. The time steps close to the sampling point are time-averaged in the same way as
in PDSC and this process is repeated a number of times, thus building up a combination of ensemble-averaged and time-averaged data without having to simulate from zero flow time for each run. This process reduces the statistical scatter in the results by adding to the number of particles in the sample, rather than by some artificial smoothing process. Because only a short period of the flow is processed in this way, the scheme has significant memory and computational advantages over both ensemble-averaging and using a greater number of particles in the time-averaging scheme.