應用非結構性網格之通用平行化三維DSMC程式(PDSC) 的研究與發展

全文

(1)國立交通大學機械工程學系博士論文. 應用非結構性網格之通用平行化三維 DSMC 程式(PDSC) 的研究與發展 Development of a General-Purpose Parallel Three-Dimensional DSMC Code (PDSC) Using Unstructured Tetrahedral Mesh. 研究生：曾坤樟指導教授：吳宗信. 博士. 中華民國九十四年七月.

(2) 應用非結構性網格之通用平行化三維 DSMC 程式(PDSC) 的研究與發展 Development of a General-Purpose Parallel Three-Dimensional DSMC Code (PDSC) Using Unstructured Tetrahedral Mesh. 研究生：曾坤樟. Student：Kun-Chang Tseng. 指導教授：吳宗信. Advisor：Dr. Jong-Shinn Wu. 國立交通大學機械工程學系博士論文. A Thesis Submitted to Department of Mechanical Engineering College of Engineering National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in. Mechanical Engineering July 2005 Hsinchu, Taiwan. 中華民國九十四年七月.

(3) Development of a General-Purpose Parallel Three-Dimensional DSMC Code (PDSC) Using Unstructured Tetrahedral Mesh. Student: Kun-Chang Tseng. Advisor: Dr. Jong-Shinn Wu. Department of Mechanical Engineering National Chiao-Tung University. Abstract A general-purpose parallel three-dimensional direct simulation M onte Carlo code (PDSC) using unstructured tetrahedral mesh is developed and validated in this thesis. Important features of this PDSC include parallel processing with dynamic domain decomposition, combination of variable time-step scheme with solution-based adaptive mesh refinement, conservative weighting scheme for treating trace species and chemical reaction functions for hypersonic air flows. A multi-level graph-partitioning technique is employed to adaptively decompose the computational domain according to the workload distribution among processors during runtime, which can alleviate the unbalancing loading before it becomes a problem. A three-dimensional h-refined unstructured adaptive mesh with simple mesh quality control, based on a preliminary DSM C simulation, is used to obtain suitable mesh resolution to increase the accuracy of the DSM C solution. A variable time-step method using the concept of fluxes (mass, momentum and energy) conservation across the cell interface is implemented to reduce the number of simulated particles and the number of iterations of transient period to reach steady state, without sacrificing the solution accuracy. A conservative weighting scheme, ensuring exact momentum and near exact energy conservation during each particle collision, is incorporated into the PDSC to efficiently treat flows with trace species. This method is validated and shows it can greatly reduce both the number of particles and computational time. Chemical reaction module, including dissociation, recombination and exchange reactions, is incorporated into the PDSC for treating reactive flows. It is verified by comparing probability, degree of dissociation and mole i.

(4) fractions of a single 2-D cell with theoretical data. Completed PDSC is then applied to compute several complicated, challenging flow problems to demonstrate its superior computational capability. The results are also validated with experimental data or previous simulation data wherever available. Organization of this thesis is briefly described as follows. Chapter 1 describes the background, motivation and objectives of the current study. Chapter 2 describes the general D SM C method and overview of the current implementation of the PDSC. Chapter 3 introduces the variable time-step scheme in combination with solution-based adaptive mesh refinement on unstructured tetrahedral mesh. Chapter 4 describes the parallel implementation and performance of the DSM C method using dynamic domain decomposition. Chapter 5 describes the conservative weighting scheme and its superiority in treating flows having trace species. Chapter 6 describes the chemical reaction functions for treating hypersonic air flows along with its validation using single-cell simulation. Chapter 7 describes the results of simulating several challenging flow problems using the PDSC. Chapter 8 concludes and summarizes the important findings of the current study, along with recommended future studies.. Keywords: direct simulation M onte Carlo, unstructured tetrahedral mesh, parallel processing, dynamic domain decomposition, variable-time-step, adaptive mesh refinement, conservative weighting scheme, chemical reaction, graph-partitioning. ii.

(5) 應用非結構性網格之通用平行化三維 DSMC 程式(PDSC)的研究與發展. 學生: 曾坤樟. 指導教授: 吳宗信博士. 機械工程學系博士班國立交通大學. 中文摘要本論文發展並驗證一通用平行化之三維直接模擬蒙地卡羅程式(PDSC)於非結構性四面體網格。這個 PDSC 程式有幾個重要的特色：以動態區域切割之平行化處理、結合變時步方法於可調適網格、可有效率處理稀少氣體之權值守恆法 (conservative weighting scheme)以及處理高速流之化學反應等等。此程式利用多層式圖形切割技術，在模擬的過程中根據每顆處理器的負載做動態區域切割。此外，利用 DSM C 的初步結果，三維非結構可調適網格(h-refined)法以及簡易的網格品質控制可以用來增加 DSM C 的準確性。變時步法則利用通過界面的流量(質量、動量以及能量)守恆之觀念，在避免犧牲結果準確性的前提下，同時減少模擬分子數目以及到達穩態的疊代次數。確保動量及能量在每次碰撞後均能守恆的權值守恆法 (CWS)可以用來有效的模擬具有稀少氣體的流場。這個方法被證實可以大量減少模擬分子的數目及計算時間。分離(Dissociation)、結合(Recombination)以及交換 (Exchange)的化學反應亦被結合在 PDSC 中來處理有化學作用的流場。最後利用發展完全的 PDSC 模擬數個複雜、具有挑戰性的流場來顯示它優越的計算能力。計算結果與實驗結果或前人的模擬結果進行比較，驗證程式的正確性。本論文的結構簡述如下：第一章描述此研究之背景、動機以及目的。第二章為 DSM C 法之簡介及目前 PDSC 軟體結構內容之概要描述。第三章介紹結合變時步法之非結構性可調適四面體網格。第四章為利用動態區域切割之平行化 D SM C 程式。第五章介紹權值守恆法以及其處理稀少氣體流場之優點。第六章為化學反應於高速流以及利用單一網格之驗證。第七章為利用 PDSC 模擬幾個具有挑戰性. iii.

(6) 的應用及結果。第八章為此研究的重要發現，以及對未來發展的建議。. 關鍵字: 直接模擬蒙地卡羅法，非結構四面體網格，平行化處理，動態區域切割，變時步法，可調適網格，權值守恆法，化學反應，圖形切割. iv.

(7) Acknowledge ments I would like to thank my research advisor, Professor Jong-Shinn Wu, for his support and guidance for this work. He gave me continuous encouragement and assistance to my research and life. It lets me has interests and keep studying during this research. I also learn how to figure problems out by my own and have a chance of one year study in the United States. It is real a turning point of my life. My graduate colleagues and research associates are also unforgettable. Our continuous discussions are significantly efficient and helpful to improve the performances of this research. You also make my research span became more colorful with your companies. I am grateful to Professor I. D. Boyd and those colleagues in his group at the University of M ichigan in the United States for giving me chance to study chemical reaction function of the DSM C method. It is a wonderful year to study with them. I also want to thank the National Science Council (NSC) and the National Center for High-Performance Computing (NCHC) in Taiwan. These simulations with immense computing cost cannot be finished without your supports. Finally, a thankyou gives to all my friends and family for all their supports. This thesis would not complete without you.. v.

(8) Table of Contents Abstract ...........................................................................................i 中文摘要 ........................................................................................iii Acknowledgements .......................................................................... v Table of Contents ............................................................................ vi List of Tables ................................................................................... x List of Figures................................................................................ xii Nomenclature ............................................................................... xxi Chapter 1 Introduction .................................................................... 1 1.1 Background and Motivation of the Thesis............................................................. 1 1.2 Reviews of the DS MC Method................................................................................ 2 1.2.1 General Features ............................................................................................... 2 1.2.2 Existing DS MC S oftware.................................................................................. 3 1.2.2.1 Visual DS MC Program .............................................................................. 3 1.2.2.2 DAC ............................................................................................................. 3 1.2.2.3 MONACO ................................................................................................... 4 1.2.2.4 S MILE......................................................................................................... 4 1.2.3 Structured and Unstructured Mesh in DS MC................................................ 4 1.2.4 Adapti ve Mesh Refinement............................................................................... 5 1.2.5 Variable Time-Step S cheme.............................................................................. 9 1.2.6 The Parallel DS MC Method........................................................................... 10 1.2.6.1 Domain Decomposition............................................................................ 10 1.2.6.2 Static Domain Decomposition (S DD)...................................................... 14 1.2.6.3 Dynamic Domain Decomposition (DDD)............................................... 14 1.2.7 Conservative Weighting S cheme.................................................................... 15 1.2.8 Chemical Reactions......................................................................................... 16 1.3 Objectives of the Thesis.......................................................................................... 17. Chapter 2 An Overview of the Current Implementation of the DSMC Method.......................................................................... 18 2.1 The Boltzmann Equation....................................................................................... 18 2.2 General Description of the Standard DS MC ....................................................... 19 vi.

(9) 2.3 The Standard DS MC Procedures.......................................................................... 20 2.3.1 Initialization..................................................................................................... 21 2.3.2 Particle Movement.......................................................................................... 21 2.3.3 Indexing............................................................................................................ 21 2.3.4 Gas-Phase Collisions ....................................................................................... 22 2.3.5 S ampling........................................................................................................... 22 2.4 Overview of the Current Implementation of PDSC............................................ 23 2.4.1 MuST Visual Preprocessor............................................................................. 24 2.4.2 PDSC................................................................................................................. 24 2.4.3 Postprocessor ................................................................................................... 26 2.5 Concluding Remarks.............................................................................................. 26. Chapter 3 Unstructured Adaptive Mesh Refinement with Variable Time-Step Scheme.......................................................... 27 3.1 Variable Time-Step S cheme................................................................................... 27 3.2 Unstructured Adapti ve Mesh Refinement............................................................ 30 3.2.1 General Features ............................................................................................. 30 3.2.2 Adaptation Parameter and Criteria............................................................... 31 3.2.3 Adaptation Procedures.................................................................................... 31 3.2.3.1 Two-Dimensional Adaptation .................................................................. 33 3.2.3.2 Three-Dimensional Adaptation ............................................................... 34 3.3 Verifications of Unstructured Adaptive Mesh Refinement and Variable Time-Step S cheme ........................................................................................................ 35 3.3.1 Two-Dimensional Flows .................................................................................. 35 3.3.2 Three-Dimensional Flows ............................................................................... 41 3.4 Concluding Remarks.............................................................................................. 44. Chapter 4 Parallel Computing of DSMC ......................................... 46 4.1 The Parallel DS MC Method.................................................................................. 46 4.1.1 Static Domain Decomposition (S DD)............................................................. 49 4.1.2 Dynamic Domain Decomposition (DDD)...................................................... 50 4.1.2.1 Decision Policy for Repartitioning.......................................................... 50 4.1.2.2 Repartition the Domain ........................................................................... 51 4.1.2.3 Cell/Particle Migration ............................................................................ 52 4.2 Parallel Performance of the Parallel DS MC Method ......................................... 54 4.2.1 Flow Conditions of Driven Cavity Flow........................................................ 54 4.2.2 Simulations on IBM-S P2................................................................................. 56 4.2.3 Simulations on IBM-S MP............................................................................... 61 4.3 Verifications of the Parallel DS MC....................................................................... 61 4.3.1 Two-Dimensional Flows .................................................................................. 62 4.3.2 Three-Dimensional Flows ............................................................................... 64 vii.

(10) 4.4 Concluding Remarks.............................................................................................. 65. Chapter 5 Conservative Weighting Scheme ..................................... 66 5.1 Conservative Weighting S cheme (CWS ).............................................................. 66 5.2 Verifications of Conservative Weighting S cheme................................................ 67 5.2.1 Maxwell-Boltzmann Distribution of A Single Cell ....................................... 68 5.2.2 Hypersonic Flow over a Quasi-2-D Cylinder................................................ 69 5.3 Concluding Remarks.............................................................................................. 69. Chapter 6 Chemical Reactions........................................................ 71 6.1 Chemical Reactions................................................................................................ 71 6.1.1 Dissociation Reaction ...................................................................................... 72 6.1.2 Recombination Reaction................................................................................. 74 6.1.3 Exchange Reaction .......................................................................................... 75 6.2 Verifications of Chemical Reactions ..................................................................... 75 6.2.1 A S ingle Cell with One-Directional Reaction................................................ 76 6.2.2 A S ingle Cell with Pure Gas (Dissociation and Recombination)................. 80 6.2.3 A S ingle Cell with 5 S pecies with 34 Reactions............................................. 81 6.3 Concluding Remarks.............................................................................................. 83. Chapter 7 Applications and Examples............................................. 84 7.1 Near-Continuum Underexpanded Twin-jet Interaction..................................... 84 7.2 Hypersonic Flow around a Three-Dimensional Apollo Re-entry Vehicle ......... 86 7.3 A Three-Dimensional Re-entry S phere Flow ....................................................... 87 7.4 Concluding Remarks.............................................................................................. 88. Chapter 8 Conclusions ................................................................... 90 8.1 Summary................................................................................................................. 90 8.2 Recommended Future Studies............................................................................... 91. References ..................................................................................... 92 Appendix A Derivation of the Probability of Dissociation/Exchange.. 99 Appendix B Derivation of the Probability of Recombination .......... 101. viii.

(11) Autobiography............................................................................. 260 Publications ................................................................................. 261. ix.

(12) List of Tables Table 1. 1 Comparison of some well-known DSM C codes. ........................................ 103. Table 3. 1 The complete listing of physical and VHS parameters of a hypersonic flow over a cylinder. ........................................................................................... 104 Table 3. 2 The cell numbers of adaptive mesh with or without cell quality control.... 105 Table 3. 3 The complete listing of physical and VHS parameters of a hypersonic flow o over a 15 -compression ramp..................................................................... 106. Table 3. 4 The cell numbers of adaptive mesh with or without cell quality control.... 107 Table 3. 5 The complete listing of physical and VHS parameters of a hypersonic flow over a sphere............................................................................................... 108 Table 3. 6 The cell numbers of adaptive mesh with or without cell quality control.... 109. Table 4. 1 Comparisons of the parallel machines at NCHC......................................... 110 Table 4. 2 Number of cells and particles for three different problem sizes for the driven cavity flow.................................................................................................. 111. Table 6. 1 The constants of the rate coefficients in chem.dat file. ............................... 112 Table 6. 2 Comparison of the reaction probability and rate constant of N 2 + N 2 → N + N + N 2 reaction ............................................................. 113 Table 6. 3 Comparison of the reaction probability and rate constant of N 2 + N → N + N + N reaction ................................................................. 114 Table 6. 4 Comparison of the reaction probability and rate constant of N + N + N2 → N 2 + N2 reaction ............................................................... 115 Table 6. 5 Comparison of the reaction probability and rate constant of N + N + N → N 2 + N reaction ................................................................ 116 Table 6. 6 Comparison of the reaction probability and rate constant of O2 + O2 → O + O + O2 reaction................................................................ 117 Table 6. 7 Comparison of the reaction probability and rate constant of O2 + O → O + O + O reaction .................................................................. 118 Table 6. 8 Comparison of the reaction probability and rate constant of O + O + O2 → O2 + O2 reaction................................................................ 119 x.

(13) Table 6. 9 Comparison of the reaction probability and rate constant of O + O + O → O2 + O reaction .................................................................. 120 Table 6. 10 Comparison of the reaction probability and rate constant of NO + N 2 → N + O + N 2 reaction ............................................................. 121 Table 6. 11 Comparison of the reaction probability and rate constant of NO + O2 → N + O + O2 reaction.............................................................. 122 Table 6. 12 Comparison of the reaction probability and rate constant of NO + NO → N + O + NO reaction........................................................... 123 Table 6. 13 Comparison of the reaction probability and rate constant of NO + N → N + O + N reaction................................................................ 124 Table 6. 14 Comparison of the reaction probability and rate constant of NO + O → N + O + O reaction................................................................. 125 Table 6. 15 Comparison of the reaction probability and rate constant of N + O + N 2 → NO + N 2 reaction ............................................................. 126 Table 6. 16 Comparison of the reaction probability and rate constant of N + O + O2 → NO + O2 reaction.............................................................. 127 Table 6. 17 Comparison of the reaction probability and rate constant of N + O + NO → NO + NO reaction........................................................... 128 Table 6. 18 Comparison of the reaction probability and rate constant of N + O + N → NO + N reaction................................................................ 129 Table 6. 19 Comparison of the reaction probability and rate constant of N + O + O → NO + O reaction................................................................. 130 Table 6. 20 Comparison of the reaction probability and rate constant of N 2 + O → NO + N reaction ..................................................................... 131 Table 6. 21 Comparison of the reaction probability and rate constant of NO + N → N 2 + O reaction ..................................................................... 132 Table 6. 22 Comparison of the reaction probability and rate constant of NO + O → O2 + N reaction...................................................................... 133 Table 6. 23 Comparison of the reaction probability and rate constant of O2 + N → NO + O reaction...................................................................... 134 Table 6. 24 M ole percentage of each species by using original chem.dat.................... 135. xi.

(14) List of Figures Figure 1. 1 Effective limits of major approximations in the DSM C method [10, 11]. 136 Figure 1. 2 Sketch of graph and mesh. ......................................................................... 137. Figure 2. 1 The flowchart of the standard DSM C method........................................... 138 Figure 2. 2 Schematic diagram of the PDSC................................................................ 139 Figure 2. 3 M uST Visual Preprocessor for numerical simulations [73]....................... 140 Figure 2. 4 Important features of the PDSC (a) future PDSC; (b) planned study of PDSC. .................................................................................................................... 141. Figure 3. 1 Sketch of the concept of variable-time-step scheme.................................. 142 Figure 3. 2 The flowchart of the DSM C method with mesh adaptation...................... 143 Figure 3. 3 M esh refinement scheme for unstructured quadrilateral mesh [39]. ......... 144 Figure 3. 4 M esh refinement scheme for unstructured triangular mesh [39]. .............. 145 Figure 3. 5 The diagram of mesh adaptation of tetrahedral cell [40]. .......................... 146 Figure 3. 6 The tree diagram of removing hanging nodes of tetrahedral mesh adaptation [40]. ............................................................................................................ 147 Figure 3. 7 Removal of hanging nodes in mesh adaptation procedures for two-dimensional triangular and three-dimensional tetrahedral meshes [40]. .................................................................................................................... 148 Figure 3. 8 Sketch of the computational domain of a nitrogen hypersonic flow over a cylinder (N2 gas, Kn∞=λ∞/D=0.025, M∞=20, T∞=20 K, n∞=5.1775E19 particles/m3, D= 1m). ................................................................................. 149 Figure 3. 9 Evolution of 2-D, unstructured triangular cells for a hypersonic cylinder flow with cell quality control (a) initial (7,025); (b) level-2 (33,773); (c) level-4 (75,099). ......................................................................................... 150 Figure 3. 10 Normalized density contour of a hypersonic cylinder flow with different meshes. ....................................................................................................... 151 Figure 3. 11 Normalized temperature contours of a hypersonic cylinder flow with different meshes (a) translational; (b) rotational; (b) total. ........................ 152 Figure 3. 12 Normalized number density and temperatures along the stagnation line with different meshes. ................................................................................ 153 Figure 3. 13 Comparison of the adaptive mesh with cell quality control for a hypersonic. xii.

(15) cylinder flow............................................................................................... 154 Figure 3. 14 Normalized number density and temperatures along the stagnation line with cell quality control.............................................................................. 155 Figure 3. 15 Particle distribution of a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes. .......................................... 156 Figure 3. 16 Comparison of particle count of a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes. .......................... 157 Figure 3. 17 Comparison of normalized density contour of a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes (x-y plane).......................................................................................................... 158 Figure 3. 18 Comparison of normalized temperature contours of a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes (x-y plane). ................................................................................................. 159 Figure 3. 19 Comparison of normalized density and temperatures of a hypersonic cylinder flow along the stagnation line using variable time-step (VTS) and constant time-step (CTS) schemes. ............................................................ 160 Figure 3. 20 Schematic diagram of flow features in a typical hypersonic flow over a o hypersonic 15 -compression ramp (N2 gas, Kn∞=λ∞/D=2.E-4, M∞=14.36,. ρ ∞=5.221E-4 kg/m3, T∞=84.83 K, Xc=43.891 cm, Xr=36.86 cm). ............ 161 Figure 3. 21 Evolution of 2-D, unstructured triangular cells for a hypersonic flow over a o hypersonic 15 -compression ramp with adaptive refinement (a) initial. (15,063); (b) level-1 (30,219); (b) level-2 (83,776). .................................. 162 o Figure 3. 22 Normalized density contour over a hypersonic 15 -compression ramp with. different meshes (a) initial (15,063); (b) level 2 (83,776).......................... 163 Figure 3. 23 Pressure, shear and heat transfer coefficient distributions along the solid o wall for a hypersonic 15 -compression ramp............................................. 164. Figure 3. 24 Comparison of the adaptive mesh with cell quality control for a hypersonic o 15 -compression ramp................................................................................ 165. Figure 3. 25 Comparison of normalized density contour with cell quality control for a hypersonic 15o-compression ramp (a) with CQC; (b) without CQC......... 166 Figure 3. 26 Sketch of a hypersonic flow over a one to sixteenth sphere (N2 gas, Kn∞=0.01035, D=1.28 cm, Tw=300 K, T∞=66.25 K, M∞=4.2)................... 167 Figure 3. 27 Normalized density contours of a 3-D, unstructured tetrahedral cell for a. xiii.

(16) hypersonic sphere flow with three different cross sections........................ 168 Figure 3. 28 Evolution of a 3-D, unstructured tetrahedral cells for a hypersonic sphere flow (a) initial (5,353); (b) level-1 (22,510); (c) level-2 (164,276). .......... 169 Figure 3. 29 Comparison of normalized density of a hypersonic sphere flow with different meshes (a) mesh distribution; (b) normalized contour. ............... 170 Figure 3. 30 Comparison of normalized temperatures of a hypersonic sphere flow with different meshes (a) translational; (b) rotational; (c) total.......................... 171 Figure 3. 31 Comparison of normalized density of a hypersonic sphere flow along the stagnation line with different meshes (a) initial; (b) level-2. ..................... 172 Figure 3. 32 Comparisons of the adaptive mesh and normalized density contour with or without cell quality control at x-y plane for a hypersonic sphere flow...... 173 Figure 3. 33 Density distribution of a hypersonic sphere flow along the stagnation line with or without cell quality control............................................................ 174 Figure 3. 34 Particle distribution of a hypersonic sphere flow using variable time-step (VTS) and constant time-step (CTS) schemes. .......................................... 175 Figure 3. 35 Comparison of particle count of a hypersonic sphere flow using variable time-step (VTS) and constant time-step (CTS) schemes. .......................... 176 Figure 3. 36 Comparison of normalized density contour over a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes.... 177 Figure 3. 37 Comparison of normalized temperature contour over a hypersonic cylinder flow using variable time-step (VTS) and constant time-step (CTS) schemes (a) translational; (b) rotational; (c) total..................................................... 178 Figure 3. 38 Comparison of normalized density of a hypersonic sphere flow along the stagnation line using variable time-step (VTS) and constant time-step (CTS) schemes. ..................................................................................................... 179. Figure 4. 1 The flowchart of the parallel D SM C method............................................. 180 Figure 4. 2 Schematic diagram of the proposed cell numbering scheme (mn is the th number of cells in the n processor, where the starting and ending cell n −1. number are. ∑ mi + 1 i= 1. n. and. ∑ mi , respectively. np is the total number of i= 1. processors) [84].......................................................................................... 181 Figure 4. 3 The flowchart of the parallel D SM C method with dynamic domain decomposition method. .............................................................................. 182 xiv.

(17) Figure 4. 4 Sketch of the dynamic domain decomposition method............................. 183 Figure 4. 5 Sketch of a two-dimensional high-speed driven cavity flow (Ar gas, Vp=8*Cmp, Tw=300 K, L/H=1, L=0.32 m, Kn=0.04)................................... 184 Figure 4. 6 Parallel speedup and efficiency as a function of number of processors for high-speed driven cavity flow at different problem sizes on IBM -SP2 machine (maximum 64 processors)............................................................ 185 Figure 4. 7 Normalized computational time per particle on a single IBM -SP2 processor. .................................................................................................................... 186 Figure 4. 8 Evolution of domain decomposition for large problem size using 64 processors, when activating SAR scheme at intervals of 2∆t, during the simulation for a bottom, lid-driven cavity flow (a) initial; (b) intermediate; (c) final............................................................................................................. 187 Figure 4. 9 Evolution of domain decomposition for large problem size using 64 processors, when activating SAR scheme at intervals of 2∆t, during the simulation for a bottom, lid-driven cavity flow. (a) initial; (b) intermediate; (c) final. ...................................................................................................... 188 Figure 4. 10 Number of particles in each processor and the number of repartitions as a function of the number of simulation time-steps for the large problem size using 16 processors when activating SAR at intervals of 2∆t.................... 189 Figure 4. 11 Number of repartitions as a function of simulation time-steps at 16 processors................................................................................................... 190 Figure 4. 12 Number of repartitions as a function of simulation time-steps at 64 processors................................................................................................... 191 Figure 4. 13 Final normalized particle numbers on each processor for three problem sizes at 64 processors................................................................................. 192 Figure 4. 14 Final normalized particle numbers on each processor for three problem sizes at 64 processors................................................................................. 193 Figure 4. 15 Fraction of time for the DSM C computation and repartition as a function of number of processors.................................................................................. 194 Figure 4. 16 Real CPU running time of doing useful DSM C work and repartition per time with different problem size................................................................. 195 Figure 4. 17 Relative costs within DSM C with different problem size on IBM -SP2 (a) small problem; (b) medium problem; (c) large problem. ........................... 196. xv.

(18) Figure 4. 18 Degree of imbalance as a function of number of processors for static and dynamic domain decomposition................................................................. 197 Figure 4. 19 Parallel speedup and efficiency as a function of number of processors for high-speed driven cavity flow at different problem sizes on IBM -SMP machine (maximum 128 processors).......................................................... 198 Figure 4. 20 Partition development for a hypersonic cylinder flow (64 CPUs) (a) initial; (b) medium; (c) final. ................................................................................. 199 Figure 4. 21 Comparison of the normalized density contour of a hypersonic cylinder flow using static and dynamic domain decomposition methods................ 200 Figure 4. 22 Comparison of the normalized temperature contours of a hypersonic cylinder flow using static and dynamic domain decomposition methods.. 201 Figure 4. 23 Normalized density and temperatures of a hypersonic cylinder flow along the stagnation line using static and dynamic domain decomposition methods. .................................................................................................................... 202 o Figure 4. 24 Partition development for a 15 -compression ramp flow (64 CPUs) (a). initial; (b) medium; (c) final. ...................................................................... 203 Figure 4. 25 Comparison of normalized density contour of a hypersonic o 15 -compression ramp flow using static and dynamic domain decomposition. methods. ..................................................................................................... 204 Figure 4. 26 Comparison of normalized temperature contours of a hypersonic o 15 -compression ramp flow using static and dynamic domain decomposition. methods. ..................................................................................................... 205 Figure 4. 27 Comparison of the coefficients along the solid wall of a hypersonic o 15 -compression ramp flow using static and dynamic domain decomposition. methods (a) pressure; (b) shear stress; (c) heat transfer. ............................ 206 Figure 4. 28 Development of partition of a supersonic sphere flow using dynamic domain decomposition method on the surface of simulation domain........ 207 Figure 4. 29 Comparison of normalized density contour of a hypersonic sphere flow by using static and dynamic domain decomposition methods........................ 208 Figure 4. 30 Comparison of normalized temperature contours of a hypersonic sphere flow by using static and dynamic domain decomposition methods........... 209 Figure 4. 31 Comparison of normalized density along the stagnation line of a hypersonic sphere flow by using static and dynamic domain decomposition methods. ..................................................................................................... 210 xvi.

(19) Figure 5. 1 Schematic diagram of CWS for non-reactive flow. ................................... 211 Figure 5. 2 Simulation of two-component with weight ratio 1:9 (a) velocity distribution; (b) relative error as a function of the number of simulation time-steps..... 212 Figure 5. 3 Simulation of two-component with weight ratio 1:19 (a) velocity distribution; (b) relative error as a function of the number of simulation time-steps..... 213 Figure 5. 4 Simulation of two-component with weight ratio 1:99 (a) velocity distribution; (b) relative error as a function of the number of simulation time-steps..... 214 Figure 5. 5 Simulation of three-component with weight ratio 1:99:9900 (a) velocity distribution; (b) relative error as a function of the number of simulation time-steps.................................................................................................... 215 Figure 5. 6 A 3-D, unstructured tetrahedral cells for a hypersonic cylinder flow. ....... 216 Figure 5. 7 Contours of number density of Ar-He mixture gas with mole fraction 1:99 and 99:1 (a) Ar:He=1:99; (b) Ar:He=99:1.................................................. 217 Figure 5. 8 Contours of number density of Ar-Ne mixture gas with mole fraction 1:99 and 99:1 (a) Ar:Ne=1:99; (b) Ar:Ne=99:1.................................................. 218 Figure 5. 9 Contours of number density of O2-N2 mixture gas with mole fraction 1:99 and 99:1 (a) O2:N2=1:99; (b) O2:N2=99:1. ................................................. 219 Figure 6. 1 The flow chart of chemical reaction in the PDSC. .................................... 220 Figure 6. 2 Comparison of simulation and theoretical data of N 2 + N 2 → N + N + N 2 reaction (a) reaction probability; (b) rate constant..................................... 221 Figure 6. 3 Comparison of simulation and theoretical data of N 2 + N → N + N + N reaction (a) reaction probability; (b) rate constant..................................... 222 Figure 6. 4 Comparison of simulation and theoretical data of N + N + N2 → N 2 + N2 reaction (a) reaction probability; (b) rate constant..................................... 223 Figure 6. 5 Comparison of simulation and theoretical data of N + N + N → N 2 + N reaction (a) reaction probability; (b) rate constant..................................... 224 Figure 6. 6 Comparison of simulation and theoretical data of N 2 + N 2 → N + N + N 2 reaction (a) reaction probability; (b) rate constant..................................... 225 Figure 6. 7 Comparison of simulation and theoretical data of O2 + O → O + O + O reaction (a) reaction probability; (b) rate constant..................................... 226 Figure 6. 8 Comparison of simulation and theoretical data of O + O + O2 → O2 + O2 xvii.

(20) reaction (a) reaction probability; (b) rate constant..................................... 227 Figure 6. 9 Comparison of simulation and theoretical data of O + O + O → O2 + O reaction (a) reaction probability; (b) rate constant..................................... 228 Figure 6. 10 Comparison of simulation and theoretical data of NO + N 2 → N + O + N 2 reaction (a) reaction probability; (b) rate constant..................................... 229 Figure 6. 11 Comparison of simulation and theoretical data of NO + O2 → N + O + O2 reaction (a) reaction probability; (b) rate constant..................................... 230 Figure 6. 12 Comparison of simulation and theoretical data of NO + NO → N + O + NO reaction (a) reaction probability; (b) rate constant..................................... 231 Figure 6. 13 Comparison of simulation and theoretical data of NO + N → N + O + N reaction (a) reaction probability; (b) rate constant..................................... 232 Figure 6. 14 Comparison of simulation and theoretical data of NO + O → N + O + O reaction (a) reaction probability; (b) rate constant..................................... 233 Figure 6. 15 Comparison of simulation and theoretical data of N + O + N 2 → NO + N 2 reaction (a) reaction probability; (b) rate constant..................................... 234 Figure 6. 16 Comparison of simulation and theoretical data of N + O + O2 → NO + O2 reaction (a) reaction probability; (b) rate constant..................................... 235 Figure 6. 17 Comparison of simulation and theoretical data of N + O + NO → NO + NO reaction (a) reaction probability; (b) rate constant..................................... 236 Figure 6. 18 Comparison of simulation and theoretical data of N + O + N → NO + N reaction (a) reaction probability; (b) rate constant..................................... 237 Figure 6. 19 Comparison of simulation and theoretical data of N + O + O → NO + O reaction (a) reaction probability; (b) rate constant..................................... 238 Figure 6. 20 Comparison of simulation and theoretical data of N 2 + O → NO + N reaction (a) reaction probability; (b) rate constant..................................... 239 Figure 6. 21 Comparison of simulation and theoretical data of NO + N → N 2 + O reaction (a) reaction probability; (b) rate constant..................................... 240 Figure 6. 22 Comparison of simulation and theoretical data of NO + O → O2 + N reaction (a) reaction probability; (b) rate constant..................................... 241 Figure 6. 23 Comparison of simulation and theoretical data of O2 + N → NO + O reaction (a) reaction probability; (b) rate constant..................................... 242 Figure 6. 24 Degree of dissociation for idea dissociating gas (a) nitrogen; (b) oxygen. xviii.

(21) .................................................................................................................... 243 Figure 6. 25 Equilibrium composition of air at density of 10-2 atm (a) original; (b) new fitting chem.dat. (lines= simulation; symbols=rate equation analysis)...... 244. Figure 7. 1 Sketch of a flow of twin-jet interaction (N2 gas, Po=870 Pa, To=285 K, L=W=10D, L=20D, D=3 mm, Knth=0.00385)............................................ 245 Figure 7. 2 The mesh of level-2 adaptation along x-y and y-z planes for the twin-jet interaction (657,624). ................................................................................. 246 Figure 7. 3 Surface local cell Knudsen number distribution on initial and level-2 adaptive meshes for the twin-jet interaction in a near-vacuum environment. .................................................................................................................... 247 Figure 7. 4 Initial and final domain decomposition for 32 processors for twin-jet interaction (a) initial; (b) final. ................................................................... 248 Figure 7. 5 Particle distribution for the twin-jet interaction by using constant and variable time-step schemes......................................................................... 249 Figure 7. 6 Normalized density contours (with respect to the density at the sonic orifice) on x-y and y-z planes using vacuum outflow boundary for twin-jet interaction................................................................................................... 250 Figure 7. 7 Normalized density contours (labels in Kelvins) on x-y and y-z planes using vacuum outflow boundary for twin-jet interaction (a) translational; (b) rotational..................................................................................................... 251 Figure 7. 8 Density and rotational temperature distribution along the symmetric centerline (y-axis) using vacuum outflow boundary for twin-jet interaction. .................................................................................................................... 252 Figure 7. 9 M esh and particle distribution of the 3-D Apollo case on X-Z plane (a) mesh; (b) particle distribution............................................................................... 253 Figure 7. 10 Number density contours of the 3-D Apollo case (a) N2; (b) O2; (c) NO; (d) N; (e) O....................................................................................................... 254 Figure 7. 11 Temperature and Velocity contours of the 3-D Apollo case (a) translational; (b) rotational; (c) vibrational; (d) U-velocity; (e) V-velocity; (f) W-velocity. .................................................................................................................... 255 Figure 7. 12 Evolution of adaptive mesh of re-entry sphere (a) initial (100,476); (b) level-2 (669,072). ....................................................................................... 256 Figure 7. 13 Normalized property contours of the re-entry sphere (a) density; (b) overall xix.

(22) temperature; (c) M ach number. .................................................................. 257 Figure 7. 14 Chemical composition along the stagnation line of the re-entry sphere.. 258 Figure 7. 15 Surface coefficients of the re-entry sphere (a) pressure; (b) skin-friction; (c) heat transfer. ............................................................................................... 259. xx.

(23) Nomenclature λ. ：mean free path. ρ. ：density. δ. ：mean molecular spacing. α. ：is a penalty parameter ：symmetric factor ：degree of dissociation ：attack angle. σ. ：the differential cross section. ω. ：viscosity temperature exponent. τ. ：surface skin-friction force. φ. 2 ：m, mv, mv /2 or other internal energy. ：free-stream parameter ：weighting ratio between cells. θ. ：characteristic temperature of dissociation ：circumferential angle. φ0. ：preset value of free-stream parameter for adaptive mesh refinement. ρd. ：characteristic density for dissociation. εrot. ：rotational energy. εv. ：vibrational energy. ζrot. ：rotational degree of freedom. ζv. ：vibrational degree of freedom. Ut ：time-step Ux ：cell size σT. ：the total cross section. Ac. ：cell area for 2-D. B. ：balance factor. c. ：the total velocity. xxi.

(24) c’. ：random velocity. co. ：mean velocity. cr. ：relative speed. C. ：repartitioning cost. Cf. ：surface skin-friction coefficient. Ch. ：surface heat transfer coefficient. Cp. ：surface pressure coefficient. d. ：molecular diameter. D. ：the throat diameter of the twin-jet interaction. dΩ. ：element of solid angle. dref ：reference diameter E. ：edge cut of graph ：energy. Ea. ：activation energy. Ec. ：cutting edge. f. ：particle phase-space distribution function. F. ：external force per unit mass. FN. ：particle weight. G. ：graph. I. ：degree of imbalance. k. ：the Boltzmann constant. kf. ：forward rate constant. kr. ：reverse rate constant. Ke. ：equilibrium constant. Kn. ：Knudsen number. Knc ：local cell Knudsen number Kncc ：preset value of cell Knudsen number for adaptive mesh refinement L. ：characteristic length; ：degree of load imbalance. xxii.

(25) m M∞ n. ：molecule mass ：free-stream mach number ：number of processor; ：speedup of parallel computing; ：number density. N. ：fluctuating number of simulated particles. Ñ. ： average number of simulated particles. Np. ：the number of simulated molecules of pth cell. Nx. ：the number of sub-domains along x-direction. Ny. ：the number of sub-domains along y-direction. Nz. ：the number of sub-domains along z-direction. p. ：surface pressure. P. ：momentum ：reaction probability. q. ：surface heat flux. r. ：a distance. Re. ：Reynolds number. S. ：weight of a sub-domain. t. ：time. T*. ：potential well-depth temperature. To. ：stagnation temperature. Tref ：reference temperature Trot ：rotational temperature Ttot ：total temperature Ttr. ：translational temperature. Tv. ：vibrational temperature. Tw. ：wall temperature. u. ：molecular velocity. xxiii.

(26) V. ：vertex of graph ：cell volume for 3-D. W. ：degradation function ：particle number of each sub-domain. Wp. ：the particle weight. Zr∞ ：limiting rotational collision number. xxiv.

(27) Chapter 1 Introduction 1.1 Background and Motivation of the Thesis Gas flows are often described correctly using the Navier-Stokes equations. However, in some flow regimes, the Navier-Stokes equations fail to approximate the gas dynamics behavior and the particle nature of the matter must be taken into account. One of these is the rarefied gas flow (Kn≥0.01, see Bird [11]), which the mean free path becomes comparable with, or even larger than, the characteristic length of flows. These high Knudsen number flows are now of practical scientific and engineering importance. For example, the examples include the pumping characteristics of turbo-molecular drag vacuum pump [1-3], the low-pressure plasma etching and chemical vapor deposition (LPCVD) [4], the micro-filter [5], and the micro-electro-mechanical-system (M EM S), [6-8], etc. Due to their importance in practical applications, accurate and efficient numerical modeling of these phenomena, rather than solving the Navier-Stokes equations, becomes necessary for understanding the underlying physics. It is well known that the Boltzmann equation is more appropriate for all flow regimes; it is, however, rarely used to numerically solve the practical problems because of two major difficulties. They include higher dimensionality (up to seven) of the Boltzmann equation and the difficulties of correctly modeling the integral collision term. An alternative method, known as Direct Simulation M onte Carlo (DSM C) method, was proposed by Bird to solve the Boltzmann equation using direct simulation of particle collision kinetics, and the associated monograph was published in 1976 [10] and 1994 [11]. Later on, both Nanbu [12] and Wagner [13] were able to demonstrate mathematically that the DSM C method is equivalent to solving the Boltzmann equation as the simulated particle numbers become large. This method has become a widely used computational tool for the simulation of gas flows in which molecular effects become important. The advantage of using particle method under these circumstances is that molecular model can be implemented directly to the calculation of particle collisions. It has been applied very successfully to compute rarefied hypersonic flows [8, 14], and other fundamental scientific problems, such as flow instabilities [15, 16]. In addition to the space science applications, it has also been utilized in the analysis of ultra-high 1.

(28) vacuum technology [1-3]. Very recently, it was applied to rarefied internal gas flow problems such as channel, pipe, ducted slider air bearing flows and the results generally agree very well with experiments [8, 14, 17]. Figure 1.1 [11] illustrates the effective limits of major approximations in the DSM C method. The dilute gas assumption requires that δ/d>>1 and δ/d=7 has been chosen as the limit, as shown as a vertical solid line. δ and d are mean molecular spacing and molecular diameter, respectively. The longer tilted dash line represents Kn=0.1 which is the demarcation between the continuum approach and the particulate approach. On the right-hand side of this line, continuum approach, such as the Navier-Stokes equations, is valid; while it is necessary to consider the particle nature of the flow on the left-hand side of this line. The shorter tilted dash line has been chosen as L/δ=100 as the criterion for the onset of significant statistical fluctuations. With the increasing importance of the rarefied gas dynamics, an efficient simulation tool becomes necessary to understand and engineer rarefied gas dynamics. Thus, in the current study we intend to develop a general-purpose DSM C code, which can be used efficiently for either fundamental or practical understanding of general rarefied gas dynamics in the modern technological development.. 1.2 Reviews of the DS MC Method 1.2.1 General Features The direct simulation M onte Carlo (DSM C) method [10, 11] is a particle method for solving the Boltzmann equation describing gas flows. The gas is modeled at the microscopic level using simulated particles, which each represents a large number of physical molecules or atoms. The physics of the gas are modeled through the motion of particles and collisions between them. M ass, momentum and energy transports between particles are considered at the particle level. The method is statistical in nature and depends heavily upon pseudo-random number sequences for simulation. Physical events such as collisions are handled probabilistically using largely phenomenological models, which are designed to reproduce real fluid behavior when examined at the macroscopic level. General procedures of the DSM C method consist of four major steps: moving, indexing, collision and sampling. In the current study, we either use Variable Hard Sphere (VH S) or Variable Soft Sphere (VSS) molecular models [10, 11] to reproduce real fluid behavior as well as No Time Counter (NTC) method [11] for the collision. 2.

(29) mechanics. Details of the procedures and the consequences of the computational approximations regarding DSM C can be found in Bird [10, 11]. In the following section, we will introduce several well-known DSM C codes, which include Visual DSMC Program, DAC, MONACO, SMILE and PDSC. The details of the standard DSM C method along with the current implementation will be introduced in Chapter 2. 1.2.2 Existing DS MC S oftware There are relatively few general-purpose DSM C codes available in the public domain. Table 1.1 summarizes the list of the main features of the four D SM C codes, including Visual DSMC Program, DAC, MONACO and SMILE. They are briefly introduced in the following in turn. Interested readers can refer to the listed references for details. 1.2.2.1 Visual DS MC Program Professor G. Bird, who invented the DSM C method back in 1958, at the University of Sydney in Australia, designed the Visual DSMC Program. It is a stand-alone Window version that can be used for two-dimensional, axis-symmetric and three-dimensional flows. It is equipped with the real-time animation of moving particles and displays of macroscopic property contours. It is an excellent tool for understanding the basics of the DSM C method and the physics of rarefied gas dynamics. However, it is not practical to utilize it for practical engineering problem due to its speed of simulation. The related information can be found in http://www.gab.com.au/. 1.2.2.2 DAC For simulating rarefied gas dynamics environment, the DSM C Analysis Code (DAC) [18] is a well-known DSM C software designed by Gerald J. LeBeau at the National Aeronautics and Space Administration (NASA). It was awarded by NASA as the 2002 Software of the Year Award and has been used extensively to support numerous space related programs, projects and missions. For example, it has been used to simulate plume impingement of space shuttle, Russian M ir Space Station, aerodynamics of re-entry space shuttle and Hubble Space Telescope (HST) servicing mission. It is a 2-D/Axi-symmetirc/3-D simulator having features such as parallel computing using dynamic load balancing, adaptive meshing and chemical reaction functions. Both distributed-memory and shared-memory can be used to run the DAC. Structured background grid coupled with unstructured triangular surface grid is used to treat complicated geometry of objects. 3.

(30) 1.2.2.3 MONACO MONACO is a well-developed DSM C code by Professor I. D. Boyd at the University of M ichigan in the United States [19]. It can be run on both workstation- and PC-class clusters. It has been applied to compute plume flows, hypersonic flows and materials processing. It is currently under development of coupling CFD and DSM C methods. The MONACO code also can be used for 2-D/Axi-symmetirc/3-D simulation. It can be run by serial or parallel computation with unstructured mesh. Important features include chemical reactions for hypersonic air flows, variable time-step scheme and “manual” dynamic domain decomposition. Some preprocessors were developed to ease the preprocessing tasks. 1.2.2.4 S MILE SMILE, stands for Statistical Modeling In the Low-density Environment, is a software system based on the DSM C method, which is developed by Professor M . Ivanov at the Institute of Theoretical and Applied M echanics, Novosibirsk, Russia [20]. It is a parallel version with dynamic domain decomposition and several numerical techniques to save the computational time. By now SMILE has been used to compute high-altitude aerodynamics problems, nozzle and plume interactions, to name a few. The detailed features are listed in Table 1.1. In the following, we will introduce and review some important aspects of implementing DSM C in practice, which more or less explain our preference to some specific choices than others. 1.2.3 S tructured and Unstructured Mesh in DS MC M ost applications of the DSM C have used structured grids [10, 11] to discretize the physical domain. For problems with complicated geometry, multi-block meshing techniques were developed first by Bird [11], which involved two steps: dividing the flow field into several blocks followed by discretizing each block into quadrilateral (2-D) or hexahedral (3-D) cells. Subsequent related research has been directed to develop alternative meshing techniques such as the coordinate transformation method by M erkle [21], the body-fitted coordinate system by Shimada and Abe [22] and the transfinite interpolation method by Olynick et al. [23]. However, all of these still used structured grids. It is much easier to program the code using structured grids; however, it often requires tremendous problem specific modification. To alleviate such restriction, an unstructured grid system is an alternative choice. M any physical problems involve very complicated geometry of objects; thus, unstructured mesh has been recommended to 4.

(31) take advantage of the flexibility of handling this situation, although it might be computationally more expensive. In addition, using unstructured mesh has the flexibility of applying graph-partitioning technique for parallel implementation of the DSM C method [24-26]. Boyd's group [27-29] has applied such technique to compute thruster plumes produced by spacecraft and found that the results are very satisfactory. Wilmoth et al. [30] have used two types of grid (unstructured tetrahedral and structured Cartesian grids) to compute the low-density, hypersonic flows about reusable launch vehicle. Both methods were shown to give comparable results. Wu's group [8, 31] has also developed 2-D and 3-D codes to compute nozzle plume and vacuum pump, respectively. It was concluded that unstructured grid has certain advantages in grid refinement as compared with structured grid. 1.2.4 Adaptive Mesh Refinement The success of the DSM C method relies on the proper distribution of simulated particles and cells. The collision partners of each cell are selected from each cell without considering their relative positions. This assumption is reasonable as long as the cell dimension is less than a mean free path. Ideally, the cell size has to be small enough, e.g., one third of the mean free path [10, 11]. Thus, solution-based mesh adaptation becomes critical step in obtaining accurate flow solution using DSM C, especially for flows having highly non-uniform density variations. The development of mesh adaptation in CFD and DSM C are described in the following for the purpose of comparison. CFD: For the past decade, the development of CFD using unstructured adaptive meshes has greatly extended the capability of predicting complex flow fields. Several adaptive mesh techniques have been developed to increase the resolution of “important” region and decrease the resolution of “unimportant” region within the flow field, as reviewed by Powell et al. [32]. In general, mesh adaptation can be categorized into three methods [33]: (1) re-meshing (mesh generation), (2) mesh movement, and (3) mesh enrichment (or h-refinement). For the first method, a solution based on the initial mesh is obtained, and then the mesh is regenerated, which the mesh points are more concentrated on where resolution of the solution is needed. This new mesh may contain more or fewer mesh points than the original mesh. For the second method, the total mesh points remain the same in the computational domain. It is common to use a spring analogy, in which the nodes of the mesh are connected by springs whose stiffness is proportional to certain 5.

(32) measure of solution activity over the spring. The mesh points are moved closer into the region where solution gradients are relatively large. This is often applied to the spatial adaptation of a structured mesh. For the third method, mesh enrichment, mesh points are added or embedded into the regions where relatively large solution gradients are detected, while the global mesh topology remains intact. It is generally regarded that mesh enrichment method has certain advantages over the first two methods [33, 34]. One of the most important advantages is that the mesh enrichment technique is in general many times faster and robust than the re-meshing technique [34]. In Ref. [33], it is mentioned that the disadvantage, however, is that the implementation of mesh enrichment involves a significant modification to existing numerical schemes due to the appearance of hanging nodes [33]. This can be easily overcome, however, by some simple methods through the elimination of hanging nodes, as proposed by Kallinderis and Vijayan [35]. DSMC: The corresponding development and the application of the adaptive mesh technique in particle method, such as the DSM C method, has been largely ignored. Applying adaptive mesh technique in the DSM C method, as in CFD, not only improves the flow field resolution without increasing the computational cost much, but also more or less equalizes the statistical uncertainties in the averaging process of obtaining the macroscopic quantities. Among the very few studies about this subject, Wang and Harvey [36] have first applied a solution-based, re-meshing adaptive grid technique (mesh regeneration using the advancing front method) in unstructured mesh to study the hypersonic flow field with highly non-uniform density variations involving shocks. Later on, in the same group, Robinson [25] has applied a similar technique combining a parallel D SM C method to compute a hypersonic flow over compression ramp at different Knudsen numbers. However, some unexpected results such as lower accuracy for a refined mesh, as compared with a coarse mesh, arose due to smaller particle-per-cell caused by too many cells. Cybyk et al. [37] have developed a technique using the M onotonic Lagrangian Grid (MLG) in the DSM C method, which provides a time-varying grid system that automatically adapts to local number densities within the flow field. However, the application of this MLG technique to external gas flows is not promising due to the particle sorting problems inhered in the scheme. Additionally, this technique highly restricts the time-step size as compared with the traditional DSM C method, which 6.

(33) makes the cost of obtaining the steady-state solution comparably high. Garcia and Bell [38] have developed an adaptive mesh and algorithm refinement (AMAR) embedding the D SM C method within a continuum method (N-S equations solver) at the finest level of an adaptive mesh refinement (AMR) hierarchy. This method can cope with problems involving physics in several orders of magnitude of length scale. Recently, Wu’s group [39, 40] also proposed 2-D and 3-D DSM C methods with adaptive mesh using h-refinement technique. Density was used as the adaptation parameter considering the statistical nature of the DSM C method. However, highly skewed cells often appear since no quality control of the mesh has been implemented. Thus, some policy during mesh adaptation should be added to improve the mesh quality. Concerns Related to M esh Adaptation Before implementing the adaptive mesh techniques, several concerns need to be considered as pointed out by the excellent review article by Powell et al. [32] and the references cited therein. These concerns are mainly applied to CFD; however, most of these are true to the DSM C method as well. These include the data structure, the initial mesh generation, the mesh adaptation procedure, the adaptation parameters and criteria, and the effect of mesh adaptation on computational algorithm. These are briefly described in turn in the following from the perspectives of the DSM C method. Data structure There exists strong relationship between the selected mesh adaptation and the data structure to be used. The majority of the DSM C codes apply structured mesh as mentioned previously. Structured mesh allows the particle tracking relatively easy and accesses the mesh information in memory more directly; however, it lacks the flexibility on mesh adaptation. There are two ways of adapting the structured mesh. One way, called r-refinement, is to “distort” the mesh distribution, so that mesh redistributes more crowded in the region where it needs mesh refinement. Another way, called h-refinement, is to add mesh in the localized region in both x- and y-directions; however, it increases unexpectedly the mesh population in other unexpected regions as well, where they do not require at all. To allow for the addition (or deletion) of mesh in the computational domain avoiding the problem outlined in the above, more sophisticated mesh data structure, rather than the simple structured mesh, has to be adopted. Unstructured mesh cannot be mapped onto a computational space with structured (i,j) indexing. Instead, the connectivity information of mesh has to be stored, which makes 7.

(34) the mesh data access indirectly. However, the spawning of the mesh in the regions of interest is much easier as compared to structured mesh due to the un-ordered data structure. It was thus concluded that unstructured mesh is superior to structured mesh considering the advantages and implementation of mesh adaptation. Initial Mesh Generation The initial un-adapted mesh required for computation is generated via either advancing front method or Delaunay triangulation if unstructured mesh is used. The detailed descriptions of these two methods can be found in Lohnern and Parikh [41] and Baker [42], respectively, and are not repeated here. Mesh Adaptation Procedure We have to decide how to adapt the mesh. For the unstructured mesh, there are generally two ways of adapting the mesh: re-meshing and embedding [36]. For the re-meshing procedure, generation of connectivity information for all the mesh is required at each adaptation step. This is expensive in general as mentioned previously. For the embedding technique, local h-refinement is used to introduce new mesh points and only the mesh in the immediate vicinity of new mesh need to be connected. Hence, it requires less computational effort. In addition, we also have to decide if isotropic or an-isotropic refinement is used. Generally, the quality of the mesh using isotropic refinement is superior to that using an-isotropic refinement. Thus, some mesh quality control policies may have to be implemented to ensure the good quality of cells. Adaptation parameters and criteria The decision of where to refine or coarsen the mesh is one of the very critical issues in mesh adaptation scheme. For the DSM C method, it often requires that the computational mesh size is much smaller (at least 1/3~1/2) than the local mean free path [10, 11], which is inversely proportional to the local number density. However, density is naturally a parameter to be considered. In addition, the choice of density as the adaptation parameter helps equalize the statistical uncertainties throughout the cells due to the sampling process for obtaining the macroscopic quantities. Effect on Computational Algorithm The effects of mesh adaptation on DSM C are trivial since the cell is used mainly for selecting collision partners and sampling the particles. Thus, the computational procedure is exactly the same except the number of the cells increases after mesh adaptation.. 8.

(35) 1.2.5 Variable Time-S tep S cheme The accuracy of a DSM C simulation is directly related to the number of simulated particles per cell throughout the cells. As the number of simulated particles increases, the statistical uncertainties of the macroscopic properties reduce due to better collision condition. The number of simulated particle per cell is shown to inversely proportional linear and square of gas density for two- and three-dimensional flows, respectively. That is, the simulated molecules are fewer in higher density regions, while lower density regions are over resolved. M ore computational time is spent calculating the lower density regions than is needed. A strategy to increase the computational speed without jeopardizing the accuracy of the solution is to reduce the number of simulated particles by using cell/particle weighting, but maintaining near-uniform particle distribution per cell, e.g., cell weighting for axisymmetric flow [10, 11], particle weighting and variable time-step [43, 44]. Kannenberg and Boyd [43] presented strategies for efficient particle resolution in DSM C. The authors manipulated variations of particle weight, variations of time-step and grid arrangement to obtain a more uniform particle count throughout the flow field. It was shown that careless use of cell/particle weighting often introduces some detrimental effects to the statistical accuracy, which is caused by repeatedly cloning the particles in the flow field [43]. Nevertheless, variable time-step method represents one of the simplest and most efficient ways of particle weighting that avoids the problem of particle cloning, if careful grid manipulation is done [43]. To obtain a more uniform distribution of model particles per cell throughout the computational domain, a variable time-step scheme is highly recommended. M arkelov and Ivanov [44] proposed a method with zoned variable time-steps in the DSM C method to simulate an axisymmetric shock wave/laminar boundary layer interaction by the DSM C method and to analyze the influence of flare length on the separation region extent. The authors employed the SMILE code to obtain adaptive mesh. The simulation domain was divided into sub-domains with different time-steps. This technique allows controlling the number of particle in each sub-domain. A relatively uniform particle distribution and 30% decrease in the total number of molecules was reported, correspondingly, to a higher computational efficiency. M oss et al. [45] have developed a local time-step scheme in the DSM C method and apply to compute a 5-deg wedge flow. A grid generation and adaptation procedure is also incorporated to insure the cell size requirement of the DSM C method. In this method, only one time-step is used throughout the flow field for unsteady flow. If the 9.

(36) flow is steady, then the computational effort can be facilitated by subdividing computational domain into an arbitrary number of regions with different local time-steps. Thus, the number of particles and overall computational time for steady state are reduced using the local time-step for each sub-domain. Usami [46] also has developed a different time-step scheme to simulate supersonic jet expansion at a very large pressure ratio. The flow filed has been divided into 12 blocks and classified into seven classes, where the further block from the orifice has a longer time-step. If the particle passes through the interface of each domain, the time-step of the further block from the orifice becomes twice as the closer block. 1.2.6 The Parallel DS MC Method The DSM C method has become a widely used computational tool for the simulation of gas flows in which molecular effects become important. The advantage of using a particle method under these circumstances is that molecular model can be applied directly to the calculation of particle collisions, while the continuum methods use macroscopic averages to account for such effects. Therefore, particle methods can predict these effects with much higher accuracy. Also, it is the only viable tool for analyzing the gas flows in the transitional regime. Nevertheless, the main drawback of such direct physical method is its high computational cost. That is why the DSM C method was only used for analyzing high Knudsen number flows (or transitional flows). For lower Knudsen number gas flows near the continuum regime, the computational cost is prohibitively high, even with the most advanced computer nowadays. Hence, it is important to increase the computational speed to extend the application range of the DSM C method. 1.2.6.1 Domain Decomposition Generally, these are two methods to partition the simulated domain, which are geometry-based and graph-based domain decomposition. Geometry-based method is usually faster but provided poor edge cut (Ec ) since they pay no heed to the connections of the point in the mesh. M any physical problems can be expressed within the framework of graph theory, such as discrete optimization problems and matrix reordering. Sketch of graph and mesh is shown in Fig. 1.2. A small portion of a triangular grid, which is usually made by commercial mesh software, is shown as the thinner lines. The bold solid circles and bold lines represent the vertices and edge cuts of the graph, respectively. The graph G(V, E) is the collection of these vertices (V) and edge cuts (E) on the basis of the connectivity between the cells. One of the advantages 10.