Objectives of the Thesis

Based on previous reviews, the ultimate goal of the current study is to develop a general-purpose three-dimensional DSM C code (PDSC) with the following features:

1. Using unstructured tetrahedral mesh for better treatment of complex geometry of boundaries.

2. Parallel processing with dynamic domain decomposition for a fast and load-balancing simulation.

3. Using variable time-step scheme, in combination with solution-based adaptive mesh, to reduce the simulated particles and to increase the accuracy of the numerical solution.

4. Using conservative weighting scheme to efficiently treat gas flows with trace species.

5. Simulating hypersonic air flows considering chemical reactions including dissociation, exchange and recombination reactions.

The thesis begins with descriptions of the conventional DSM C method and overview of the current implementation of the PDSC in Chapter 2. M esh adaptation and variable time-step scheme for an unstructured mesh will be discussed in Chapter 3.

Parallel computing of DSM C are presented in Chapter 4. Chapters 5 and 6 are description and verification of conservative weighting scheme and molecular chemical reactions of PDSC, respectively. Chapter 7 presents some applications of the current DSM C implementation to some realistic flow fields. Finally, conclusions of the current study and the recommended future studies are summarized in Chapter 8.

Chapter 2 An Overview of the Curre nt Implementation of the DSMC Method

2.1 The Boltzmann Equation

As mentioned previously, the Boltzmann equation is valid for all flow regimes. It describes the statistical distribution of particles in a fluid. It is one of the most important equations of non-equilibrium statistical mechanics, the area of statistical mechanics that deals with systems far from thermodynamics equilibrium. The Boltzmann equation is used to study how a fluid transports physical quantities such as heat and current, and thus to derive transport-related properties such as viscosity, thermal conductivity and electrical conductivity. There are some assumptions made in the derivation of the Boltzmann equation which define limits of applicability.

1. M olecular chaos is an essential component 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. Assume 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);

The meaning of particle phase-space distribution function f is the number of particles with center of mass located within a small volume d³r near the point r, and velocity within a range d³u, at time t. F_i is an external force per unit mass, t is the time and u_i is the 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. The 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.

2. The second term gives the spatial variation of the distribution function.

3. The third term describes the effect of a force on the particles.

4. The term at right hand side of the equation is called the collision integral. It is the source of most of the difficulties in obtaining solutions of the Boltzmann equation.

However, the Boltzmann equation is rarely used to numerically solve the practical problems because its higher dimensionality (up to seven) and the difficulties of modeling the integral collision term. Instead, the DSM C method has been used to simulated problems involving rarefied gas dynamics, which is the main topic in the current thesis.

2.2 General Description of the Standard DS MC

Due to the expected rarefaction caused by the very small size of micro-scale devices or the rarefied gas flows, the current research is performed using the DSM C method [10, 11], which is a particle-based method. The basic idea of DSM C is to calculate practical gas flows through the use of a method that has a physical rather than a mathematical foundation, although it has been proved that the DSM C method is equivalent to solving the Boltzmann equation [12, 13]. The assumptions of molecular chaos and a dilute gas are required by both the Boltzmann formulation and the DSM C method [10, 11]. The molecules move in the simulated physical domain so that the physical time is a parameter in the simulation and all flows are computed as unsteady flows. An important feature of DSM C 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 DSM C method are described in the next section, and the consequences of the computational approximations can be found in Bird [10, 11].

In real molecular collision, the force between molecules is strongly repulsive at short distance and weakly attractive at larger distance. M odels for analytical and numerical studies involve some degree of approximation. These models are developed to imitate the real particle collision according to experiment. There are three molecular collision models, which are the Hard Sphere (HS), Variable Hard Sphere (VHS) and Variable Soft Sphere (VSS) molecular models, in the standard DSM C method [11]. The total collision cross section of the hard sphere model is proportional to the square of the constant diameter. It has the advantage of easily calculated collision mechanics because of the isotropic scattering that means all directions are equally possible for the post-collision velocity in the center-mass frame of reference. But the cross-section should vary with relative velocity in reality. The variable hard sphere (VHS) model proposes the collision diameter is a function of relative speed, which can predict the viscosity more accurately. The cross-section is determined from the viscosity coefficient, but the ratio of the momentum to the viscosity cross-section follows the hard sphere value. Thus, the variable soft sphere (VSS) model is developed to predict the correct viscosity and diffusion coefficients, which the scattering of post-collision is not isotropic anymore.

The procedure for the collision is based on the cell, which collision pairs are chosen from the cells. The correct probability of collision between two particles is proportional to the product of their relative speed and total collision cross-section. 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.

Note that the corresponding molecular data including reference diameter (d_ref), reference temperature (T_ref), and the viscosity temperature exponent (ω) for each species are taken from those listed in Ref. [11]. Solid walls for all cases considered in this study are assumed to be fully diffusive (100% thermal accommodation), unless otherwise specified.

2.3 The S tandard DS MC Procedures

Figure 2.1 is a general flowchart of the DSM C method. Important steps of the DSM C 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.

2.3.1 Initialization

The first step to use the DSM C method in simulating flows is to set up the geometry and flow conditions. A physical space is discretized 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 M axwell-Boltzmann distribution according to the free-stream velocities and temperature, and the positions of each particle are randomly allocated within the cells.

2.3.2 Particle 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 DSM C code by Bird [10, 11], 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.

2.3.3 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 identities of the collision partners can be easy determined according to this numbering system.

2.3.4 Gas-Phase Collisions

The other most important phase of the DSM C method is gas phase collision. The current DSM C 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

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. σ_T and c are the cross section and the relative speed, respectively. The _r 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 [72], which assumes in equilibrium.

The collision process will be finished when all the collision pairs are handled for all cells and then progress to the next step.

2.3.5 S ampling

After the particle movement and collision process finishes, 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 [10, 11];

=nm

n, m are the number density and molecule mass, receptively. c, c_o, and c’ are the total velocity, mean velocity, and random velocity, respectively. In addition, T_tr, T_rot, T_v and T_tot 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 degrees 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.

2.4 Overview of the Current Implementation of PDS C

So far, there have only several well-developed DSM C codes which mentioned in Section 1.2.2. It is very helpful to develop a general-purpose parallel DSM C code for simulating problems of rarefied gas flows. Thus, development of this numerical solver is the main subject of this thesis. To make the PDSC friendlier to use, an applicable preprocessor and postprocessor are necessary. Figure 2.2 shows the overview of the planed PDSC, which includes a graphic-based preprocessor, main numerical solver and a postprocessor. They are described in the next sections in tern.

2.4.1 MuS T Visual Preprocessor

A graphical user interface (GUI), which names M uST Visual Preprocessor (Fig.

2.3(b)), is developed by Professor Wu’s group [73]. It is used to ease the parameter settings of boundary conditions, initial conditions and parallel processing because these procedures are tedious and complicated as can be shown in Fig. 2.3(a). This preprocessor can be either used for cell-based (e.g. DSM C) or for node-based (Finite Element M ethod) numerical methods. Firstly, the M uST Visual Preprocessor will ask for mesh connectivity information, which is easy to obtain from mesh generators. The mesh types can be quadrilaterals, triangles, hexahedrons, tetrahedrons, prisms and pyramids. Secondly, the preprocessor will transfer mesh format into graph and process the initial partition according to the weighting of each cell. Thirdly, the processor can create the boundary surfaces automatically and we can assign the boundary types very easily to create the input files for the PDSC.

2.4.2 PDSC

Parallel DSM C Code (PDSC) is the main solver developed 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.

Unstructured Tetrahedral M esh

Reasons of PDSC using unstructured tetrahedral mesh are: (a) it can be easily used for flows with complicated boundary conditions, (b) parallel processing can be easier implemented via graph-partitioning technique, which can handle irregular inter-processor boundary of dynamic domain decomposition, (c) it can be coupled with unstructured node-based numerical method (e.g. N-S equations).

According to these advantages of using unstructured mesh, 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. Briefly speaking, the movement of a particle is determined by the velocity and initial position of the particle. If the intersecting face is an I/O boundary, the particle will be removed.

If not, then process the interaction according to the specified wall boundary condition.

The details of particle ray-tracing techniques of two- and three-dimensional domain are described in Ref. [8, 74].

Collision Cross-Section Data

As mentioned in Section 2.2, the variable soft sphere (VSS) model can reproduce the viscosity and diffusion coefficients correctly. The relevant parameters of using VSS

model for the DSM C method can be found in Bird’s book [10, 11]. This reference provides some usual gaseous species. When the flow involves some special species, it has problem to obtain the relevant parameters of the VSS molecular model. To overcome this problem, a quantum chemistry method is proposed to calculate the intermolecular energy surface according to the distance between molecules [75]. Then the simulated intermolecular energy potential is fitted through the Lennard-Jones (L-J) potential to obtain the constants. Based on these constants and gas kinetic theory, the transport coefficients, which are viscosity and diffusion coefficients, are derived. Finally, the parameters of the VSS model are derived by fitting these computed coefficients to those derived from the VSS model.

Pressure Boundary Treatment

In order to perform accurate simulation for inflow/outflow pressure boundaries, general procedure for treating these conditions by using the concept of particles flux conservation is developed in PDSC [76]. This function is useful for applications of micro-manifold, micro-nozzle and slider air bearing.

Unstructured Adaptive M esh with Variable Time-Step Scheme

To obtain accurate simulated results, two- and three-dimensional h-refinement adaptive mesh with variable time-step scheme is developed [39, 40]. Some parameters are used to determine the adaptive level and a simple cell quality control can prevent the creation of high aspect ratio cells. This module is not only valid for PDSC, but also suitable for other numerical simulators. The detail of adaptive mesh refinement and variable time-step scheme can be found in Chapter 3.

Parallel with Dynamic Domain Decomposition

To save the enormous computational cost of the standard DSM C code, a parallel DSM C with dynamic domain decomposition. M essage passing interface (M PI) is used for data communication. This function can automatically repartition the graph domain according to the loading of each processor, which is the particle number of each cell, to achieve the load balancing of the simulation. It also can be used for other particle simulation and equation solvers. The detail of this feature is presented in Chapter 4.

Conservative Weighting Scheme

When the flow involving trace particle species, the simulation needs lots of simulated particles to satisfy the DSM C limit, which will lead to immense computational time. A weighting scheme is developed to deal with this kind of flows [77]. The basic concept is assigning the lower weight for trace particle species to create

more simulated particles. This method does not use particle cloning and destroying to avoid the statistical error. The detail of conservative weighting scheme is introduced in Chapter 5.

M olecular Chemical Reaction

Finally, PDSC also has the function to simulate flows with chemical reactions.

Chemistry is important and needs to be considered when the flow velocity and temperature is very high. The chemistry in PDSC is developed with help of Professor I.

D. Boyd at the University of M ichigan in United States. It has dissociation, exchange and recombination reactions in PDSC. Chapter 6 is a section of detailed introduction and validation of this feature.

2.4.3 Postprocessor

After the DSM C solver (PDSC), it needs a postprocessor to view the result of the simulation. The output data of PDSC can be transferred from cell-based data into node-base data easily and then importing into both Tecplot and Grapher for displaying purpose.

2.5 Concluding Remarks

This chapter has presented an overview of the current implementation of the

This chapter has presented an overview of the current implementation of the