Concluding Remarks

This chapter has presented an overview of the current implementation of the DSM C method in brief. The first part of this chapter was concerned with the general description of the standard DSM C method, which is proven that it does provide solutions that are consistent with the Boltzmann equation. The detail of the D SM C method can be found in Ref. [10, 11]. The second part describes an overview of the current implementation of PDSC. A graphical user interface, which names M uST Visual Preprocessor, is developed to deal with complicated and tedious procedures. This preprocessor not only can be used for cell-based simulators (e.g. DSM C or PIC), but also for some node-based solvers (e.g. CFD or FEM ). Then, several important features of the PDSC are also introduced to deal with different problems of flows. The next chapter will introduce the unstructured adaptive mesh with variable time-step scheme in DSM C method to improve the accuracy and efficiency of the computation.

Chapter 3

Unstructured Adaptive Mesh Refinement with Variable Time-Step Scheme

The accuracy of the DSM C method depends on the number of simulated particles and a suitable mesh. It is necessary to adapt the existing mesh according to the flow filed to achieve higher resolution of flow properties and apply the variable time-step scheme to obtain a more uniform particle distribution and efficient computation. Thus, in the current section, the general features of the proposed PDSC with variable time-step scheme combining mesh adaptation will be described in detail.

3.1 Variable Time-S tep S cheme

The number of simulated particles per cell is related to the number density, cell volume and particle weight by the following relation;

P

P W

N = nV (3.1)

N_p is the number of simulated molecules of p^th cell. W_p is the particle weight which is defined as ratio of the number of real particles to the number of simulated molecules, and n and V are the number density and the volume of the computational cell, respectively. If the number density is assumed to be a constant, the simulated particle count decreases by decreasing the cell volume or increasing the particle weight. As mentioned previously, mesh refinement can help to obtain a better cell-size distribution, ideally on the order of the local mean free path. The volume of a cell can then be related to density by the fact that the mean free path is inversely proportional to the number density. Thus, the relationship between these variables is shown as Eq. (3.2),

−1

∝

∝

∆x λ n (3.2)

For two-dimensional flow, the cell volume is given by

2

For three-dimensional flow, the cell volume is given by

3

Substituting Eqs. (3.2)~(3.4) into Eq. (3.1), give the following relation between number of simulated particles and flow density (assuming constant particle weight):

D According to Eq. (3.5), the number of simulated particle is inversely linear and square proportional to the number density with respect to the two- and three-dimensional flows.

That is, the lower density regions have larger simulated particle numbers and the higher density regions have fewer simulated molecular numbers. This effect is more obvious in three-dimensional simulation. This will lead to computational waste and incorrect results at lower and higher density regions, respectively. To avoid this problem, a variable time-step scheme is proposed to obtain a more uniform particle distribution as follows:

From Eq. (3.5), the density distribution is inversely proportional to the dimension of the cell. Thus, the first step of variable time-step scheme is to find out the cell, which has the minimum cell volume (V_min), and to calculate the local time-step of the cell as Eq. (3.6). The time-step is also proportional to x∆ and inversely proportional to the number density,

Umean and 2kT m are the mean and thermal velocity, respectively. Then, each local time-step∆ of each cell can be assigned based on t_j ∆t_min and the cell volume V_j as Eq. of the variable time-step scheme for a simulated particle moves across the cell interface is illustrated in Fig. 3.1.

Basic idea of variable time-step method in PDSC is to enforce the flux

conservation (mass, momentum and energy) of moving simulated particle when crossing the interface between two neighboring cells. If we scale the local cell time-step to the local cell size (or local mean free path), then the best way to enforce flux conservation is to change the particle weight factor without destroying or cloning the particles during particle movement across the cell interface. The cloning of particle can generally induce unpredictable random-walk effects in a statistical simulation like DSM C. One of the advantages in implementing the variable time-step scheme is to reduce both the simulated particle numbers and transient time-step to steady state, when the sampling normally starts in DSM C. This will result in appreciable time saving for the steady DSM C simulation. The net flux of the physical particles, including mass, momentum and kinetic energy, should be enforced conservation when a simulated particle crosses the cell interface form the cell 1 and to the cell 2. Thus,

t u weight, conserved flux quantity and time-step, respectively, and the numbers at subscript represents cell numbers. Note that A represents the area of cell interface between cell 1 and 2. N₂ is number of the simulated particle in cell 2, which originated from cell 1. There are several choices of the corresponding parameters to satisfy Eq.

(3.8), with which we can play. The best choice is to set N₂=1 (without particle cloning or destroy) and to keep Φ₁ =Φ₂ without changing the velocity across the cell interface, Eq. (3.8) can be rewritten as Eq. (3.9)

In other words,

i i

t W

∆ will be the same for all cells throughout the computational domain.

Inserting Eqs. (3.7)~(3.9) into Eq. (3.1), the number of simulated particles per cell is,

D

Using this approach, resulting number of simulated particles per cell for the

three-dimensional flow scales with ∆x (~³V , V_{c c} is the cell volume) [43] if cell size is proportional to the local mean free path, which otherwise scales with (∆x)². In doing so, the simulated particle will only have to adapt its weight that is proportional to the size of time-step, which is approximately commensurable to the local mean free path if solution-based adaptive mesh is used. Of course, the remaining time for a simulated particle, when crossing cell interface, should be rescaled according to the ratio of time-steps in original and destination cells. In the PDSC, the procedure of variable time-step scheme is listed in the following;

1. Chose a minimum cell volume to calculate the reference time-step, Eq. (3.6).

2. Assign the time-step for each cell based on the cell size, Eq. (3.7).

3. Determine the particle weight for each cell by Eq. (3.9).

4. The time-step has to be modified if the particle crosses the cell interface.

By using this variable time-step scheme, the simulated particle number and transient time will be reduced to speed up the computing.

3.2 Unstructured Adapti ve Mesh Refinement 3.2.1 General Features

Based on the reviews in previous chapter and considering the applications of DSM C, the general features of an adaptive mesh generation scheme are proposed as follows: (1) unstructured mesh (triangular, quadrilateral and hybrid in two-dimensional mesh and tetrahedral in three-dimensional mesh); (2) h-refinement with mesh embedding; (3) local adaptation criteria and free-stream parameter (relative local density ratio to free-stream value) as the mesh adaptation parameters; (4) upper limit on maximum number of levels of mesh adaptation. And the variable time-step scheme is incorporated into the unstructured adaptive mesh to reduce the number of simulators.

All mesh adaptation methods need some means to detect the requirement of local mesh refinement to better resolve the features in the flow fields and hence to achieve more accurate numerical solutions. This also applies to DSM C. It is important for the adaptation parameters to detect a variety of flow features but does not cost too much computational time. Often gradients of properties, such as pressure, density or velocity, are used as the adaptation parameters to detect rapid changes of the flow-field solution in traditional CFD. However, by considering the statistical nature of the DSM C method, density is adopted instead as the adaptation parameter. Using density as the adaptation

parameter in DSM C is justified since it is generally required that the mesh size be much smaller than the local mean free path to better resolve the flow features, as mentioned previously. Thus, density is used to determine the mesh adaptation in the current study.

3.2.2 Adaptation Parameter and Criteria

There are two adaptation criteria to use if the cell should be refined or not. One is the local cell Knudsen number, Kn_c, based on density. The other is the free-stream parameter φ_i. The details are described in the following.

Firstly, to use the density as an adaptation parameter, a local cell Knudsen number is defined as magnitude of local cell area and local cell volume of two- and three-dimensional domain, respectively. When the mesh adaptation module is initiated, the local cell Knudsen number for each cell is computed and compared with a preset value, Kn_cc. If this value is less than the preset value, then mesh refinement is required. If not, check the next cell until all cells are checked. This adaptation parameter is expected to be most stringent on mesh refinement (more cells are added); hence, the impact to PDSC computational cost might be high, but is required to obtain an accurate solution.

Secondly, considering the practical applications of mesh adaptation in external flows, we have added another constraint, φ_i ≥φ₀, where φ_i, free-stream parameter of each cell, is defined as

∞

= ρ

φ_i ρⁱ (3.12)

and φ₀ is a preset value. Not only the above constraint helps to reduce the total refined cell numbers to an acceptable level by reducing the cell numbers in the free-stream region a great deal, but also reduces the total computational time up to 30% as can be shown later.

3.2.3 Adaptation Procedures

In our simulations, the h-refinement has been developed to obtain a reasonable adaptive mesh. The concepts of the adaptation procedures of two- and

three-dimensional simulation are similar, but the details are complicated and non-trivial.

Before outlining the procedures of mesh adaptation, three general rules are described as follows:

(1) Isotropic mesh refinement is employed for those cells, which flag for mesh refinement. A new node is added on each edge (face) of a parent cell and connecting them to form child cells. In general, this will create hanging nodes in the interfacial cell, which is not refined, next to the isotropically refined cell. Existence of hanging node(s) not only complicates the particle movement, but also increases the cost of the cell-by-cell particle tracing due to the increase of face numbers. Hence, a remedy is proposed as follows in item (2).

(2) An-isotropic mesh refinement is utilized in the (interfacial) cells next to those cells have just been isotropically refined. The child cells are formed no matter what type of the interfacial cell is, considering the generality of the practical programming. However, some special treatment is required. The removal of hanging node(s) in the interfacial cells does increase the computational cost; it is, however, trivial as compared with the disadvantages caused by the hanging node(s).

(3) Isotropic mesh refinement is used to remove the cells of high aspect ratio. In the interfacial regions between refined and unrefined cells, bad quality cells (large aspect ratio) do occur. In DSM C, this will introduce large errors when particles are collided in the cell if it is treated the same as other normal cells. The remedy to remove these high aspect ratio cells is recorded these cells, which has been refined by an-isotropic adaptation. If these labeled cells suffering the second an-isotropic refinement, these cells will be forced to refine by isotropic adaptation.

The mesh adaptation procedures are performed after enough samples of data at each original cell are gathered. As a rule of thumb, about 50,000 particles sampled in a cell are considered enough for the mesh adaptation purpose. The mesh adaptation module is initiated and checks through all the cells to determine if mesh enrichment is required based on the specific adaptation parameter, which was explained previously. If mesh enrichment is conducted, associated neighbor identifying arrays are updated or created, coordinates and number of face for new cells are recorded, and sampled data on the coarse parent cell are redistributed (based on the magnitude of cell size) to the finer child cells accordingly. The above procedures are repeated until the prescribed maximum number of adaptation levels has been reached or no mesh enrichment is required for all the cells in the computational domain. Before preceding the D SM C

computation using the most updated mesh, all sampled data are reset to zero. Finally, the DSM C computation is then conducted on the final refined mesh to accumulate enough samples for obtaining the macroscopic properties in the cells.

In summary, the following steps summarize the procedures for mesh refinement:

1. Set up initial grids and input data.

2. Process computation until enough sampled data are gathered at each cell.

3. Compute the adaptation parameters of each cell by using Eq. (3.11) and Eq.

(3.12).

4. Refine all the cells in which both the Kn_c is less than the preset Kn_cc and φ_i is larger than φ₀ by conducting isotropic mesh refinement. If the adaptation criteria are not met, go to step (8).

5. Create and update the neighbor identifying arrays, coordinates, face numbers, and distribute sampled data to child cells, respectively. Reduce the simulation time-step to half.

6. Check if there are any hanging nodes in the interfacial cells. If it does, then conduct an-isotropic mesh refinement. Also create and update associated cell data as described in step (5).

7. Return to (2) if both the accumulated adaptation levels are less than the preset maximum value and mesh refinement is required.

8. If the accumulated adaptation levels are greater than the preset value or no mesh refinement is required, then reset all sampled data to zero and precede the DSM C computation as normal.

The corresponding flowchart is illustrated in Fig. 3.2 Note that the proposed mesh adaptation is capable of refining the mesh close to the body surface following the real surface geometry if the surface contour can be cast into parametric function format. The details of the adaptation procedure are present as follows.

3.2.3.1 Two-Dimensional Adaptation

These rules applied to quadrilateral and triangular mesh in two-dimensional simulation, as illustrated in Fig. 3.3 and Fig. 3.4, respectively. The mesh adaptation can be categorized into isotropic and an-isotropic.

1. Isotropic adaptation:

If the cell is suggested to be isotropic adaptation, quadrilateral and triangular child cells are formed with respect to the parent cell. In general, this will create one to three hanging nodes in the interfacial cell, which is not refined,

next to the isotropically refined cell. For example, for a quadrilateral interfacial cell with three or four hanging nodes (Fig. 3.3), an isotropic cell refinement is conducted and forms four quadrilateral parent cells.

2. An-isotropic adaptation:

For the an-isotropic adaptation, triangular child cells are formed no matter what type of the interfacial cell is, considering the generality of the practical programming. Typical methods of interfacial mesh refinement in the quadrilateral and triangular cells are also shown in Figs. 3.3 and 3.4, respectively. However, some special treatment is required. For example, for a quadrilateral interfacial cell with one or two hanging nodes (Fig. 3.3), an isotropic cell refinement is conducted to form three or four triangular cells, respectively.

3.2.3.2 Three-Dimensional Adaptation

The basic ideas of refining the 2-D unstructured mesh, first by isotropic refinement, then by an-isotropic refinement, remain the same for the three-dimensional tetrahedral unstructured mesh, although more complicated conditions will be encountered. The adaptation has two kinds of refinement, which are isotropic and an-isotropic refinements as 2-D adaptation. When the cell is determined to be refined, it will process the isotropic refinement, that is, the original cell will be divided to eight subcells equally.

This process will create hanging nodes for the neighboring cells. In order to remove those hanging node, an-isotropic mesh is processed to over come the difficulty of this problem. The tree diagram of mesh adaptation is shown as Fig. 3.5.

No matter the parent cell is adapted by isotropic or an-isotropic mesh refinement.

The child cells are all formed as tetrahedral cell. As Fig. 3.6 illustrates, there are six types with different hanging node numbers. The tetrahedral interfacial cell having one to five hanging node always forms tetrahedral child cells by adding appropriate node(s) on the edge(s) of the interfacial cell. The number of refined cells in the interfacial tetrahedral cell is either two, four or eight depending upon the distribution of the hanging nodes. Rules of removing hanging nodes can be briefly summarized as follows:

1. An-isotropic adaptation:

a. One hanging node (Type 1): Forming two tetrahedral child cells by directly connecting the hanging node to the other two nodes in the interfacial cell.

b. Two and three coplanar hanging nodes (Type 2a, 3b): In the case of the two coplanar hanging nodes, one node is added to the same coplanar face. The

cell, with these two types with different hanging nodes, will be refined as four tetrahedral child cells by an-isotropic adaptation.

2. Isotropic adaptation:

a. Others (Type 2b, 3a, 4, 5, 6): Forming eight tetrahedral child cells by adding appropriate nodes to fill up the six edges of the interfacial cell and then form the child cells as that in isotropic mesh refinement described earlier.

Usually, no matter for the two- or three-dimensional adaptation, the an-isotropic adaptation will form two or four child cells with high aspect ratio. These cells with bad quality maybe lead to poor simulated results. Thus, some modifications must be operated to improve the cell quality. Isotropic mesh refinement is again used to remove the cells having high aspect ratio, as shown in Fig 3.7, which represents the simplest approach of mesh quality we can think of. The method is first marked the cell with an-isotropic adaptation. Then, these cells will be forced to refine by isotropic adaptation if these labeled cells suffering second an-isotropic adaptation. It is effective in improving the mesh quality, although it is simple. Besides, one other important step during the mesh adaptation procedures is to modify the old cell neighbor-identifying arrays and add the new cell neighbor-identifying arrays accordingly.

3.3 Verifications of Unstructured Adaptive Mesh Refinement and Variable Time-S tep S cheme

The test cases include two-dimensional and three-dimensional simulations. A hypersonic cylinder flow and a 15^o-compression ramp flow are simulated in two-dimensional domain. A three-dimensional supersonic flow over a sphere is also presented as a benchmark test. Related issues about the advantages of using variable time-step method and unstructured adaptive mesh will be discussed along with the presentation of simulated results. Results are then compared with experimental data and previous simulation wherever available. In addition, physics of the flow field related to these problems will be described as brief as possible, since we are interested in

The test cases include two-dimensional and three-dimensional simulations. A hypersonic cylinder flow and a 15^o-compression ramp flow are simulated in two-dimensional domain. A three-dimensional supersonic flow over a sphere is also presented as a benchmark test. Related issues about the advantages of using variable time-step method and unstructured adaptive mesh will be discussed along with the presentation of simulated results. Results are then compared with experimental data and previous simulation wherever available. In addition, physics of the flow field related to these problems will be described as brief as possible, since we are interested in