以蒙地卡羅法模擬次音速流通過垂直平板的現象

(1)

國立交通大學

機

械

工

程

機

械

工

程

機

械

工

程

機

械

工

程

學

系

DSMC Simulation of the Subsonic Flow Past a Vertical

Plate

以蒙地卡羅法模擬次音速流通過垂直平板的現象

研

研究

究

究生

究

生

生：

：

：林士傑

林士傑

指導教授

指導教授：

：

：吳宗信

：

吳宗信

吳宗信教授

吳宗信

教授

中華民國九十

中華民國九十七

七

七年

年

年七

七

七月

七

月

(2)

以蒙地卡羅法模擬次音速流通過垂直平板的現象

DSMC Simulation of the Subsonic Flow Past a Vertical Plate

研究生：林士傑

Student: Shin-Chieh Lin

指導教授：吳宗信博士 Advisor: Dr. Jong-Shinn Wu

國立交通大學

機械工程學系

碩士論文

A Thesis

Submitted to Institute of Mechanical Engineering

College of Engineering

National Chiao Tung University

in partial Fulfillment of the Requirements

for the Degree of

Master of Science

in

Mechanical Engineering

July 2008

Hsinchu, Taiwan

中華民國九十七年七月

(3)

致

謝

在交大求學的兩年時間，特別感謝吳宗信老師的照顧，讓我在學習過程及生活方面都得到不少獲益。在研究及做學問方面，吳老師悉心指導與督促，使我從中學習到許多，處理事情以及解決問題變得更有效率，讓我順利的完成研究並且不斷地成長。還有老師愛台灣精神及推廣本土文化不遺餘力也在我心中根深蒂固。同時也感謝口試委員在口試時提供的寶貴意見，使得本論文更加充實完備。另外特別感謝曾坤樟學長長期發展關於此篇論文中所用到的程式之辛勞，且在不吝於付出時間與心力來教導我於研究中所遭遇的問題，在此一併致謝。實驗室的氣氛融洽，使我在良好的環境下學習，學習過程中更有效率。由衷地感謝洪維呈、謝昇汎、盧勁全、林宗漢、王柏勝等以畢業的學長，以及 APPL 實驗室的成員，李允民、周欣芸、李富利、許哲維、鄭凱文、胡孟樺、邱沅明、江明鴻、林雅茹、林昆模學長姐的指導，呂政霖、鄭丞志、劉育宗、吳玟琪、柳志良、蘇正勤，與你們努力奮鬥相處的時光將是我最美好的回憶，穎志、俊傑、逸民、必任等學弟妹的協助與鼓勵，以及來自紐西蘭的訪問學者 Hadley M. Cave (瓦片)和法國的 Postdoc 學者 ALIAT Abdelaziz，使我這兩年過得相當充實且溫馨，並能順利完成學業。此外，特別感謝我的大學老師牛仰堯導師，感謝他對我的期許鼓勵及愛護，以及陪伴著我的家人以及好友，有你們的鼓勵與支持，使得我更能堅持下去。在這離別的季節，大家各奔前程，希望大家追求自己的夢想前進，擁有光明的未來和生活。 林士傑謹誌 九十七年七月于風城

(4)

學生

學生:

:

: 林士傑

:

林士傑

指導教授

指導教授:

:

: 吳宗信

吳宗信

國立交通大學機械工程學系

摘要

Vortex-shedding 屬於流體力學中的外流場問題，產生的原因是當流體通過不同形狀的物體時，使物體的尾流產生週期性的剝離的渦流現象，此現象就稱 vortex-shedding。常發現在鳥在空中飛行、車子在路面上行走、橋的橋墩以及氣流受到島嶼影響等。在過去也有很多科學家做過相關的研究，但是大多數的研究 vortex-shedding 都是在連續流及不可壓縮流流場的範圍，而少數針對稀薄流體區域做研究，主要由於在稀薄流體區做實驗以及在非穩態流場模擬也較為困難。本文的目的是使用直接模擬蒙地卡羅法及結構性格網來模擬次音速流體通過垂直平板，研究 vortex-shedding 現象。我們使用 time-averaging with temporal variable time step 平均取樣時間方法的模擬，這種方法稱為 TVTS。

利用不同 Unsteady time average with temporal variable time step (TVTS)、particle per cell、number of temporal node、domain size 以及 Reynolds number 等這些參數，觀察垂直平板尾流層產生 vortex-shedding 的變化情形。

由結果顯示 TVTS=100 和 150 設定條件下，尾流層都會發生擺動的現象，而 TVTS=100 的時候，尾流層產生明顯 vortex-shedding。當固定 TVTS=100，模擬不同 Reynolds number，則會發生 Strouhal number 和 aerodynamics coefficient 會隨著 Reynolds number 增加。

(5)

DSMC Simulation of the Subsonic Flow Past a Vertical Plate

Student:

Shin-Chieh Lin

Advisor: Dr. Jong-Shinn Wu

Department of Mechanical Engineering

National Chiao-Tung University

Abstract

The phenomena of vortex shedding associated with the subsonic external flow problems in different length scales are visible everywhere in fluid dynamics. For example, aviation of fruit flies and birds, driving car in the wind, flowing river through piers under a bridge, and the air current interaction with an island and so forth. A large number of experimental and numerical studies have been reported on the vortex-shedding flows in the continuum limit, while there have been very few studies focusing on similar flows in the rarefied gas regime. Major obstacle of the investigation in rarefied regime mostly came from the difficulties of experiments and also numerical simulations for unsteady flows in this regime.

In the present paper, a general-purpose Parallel Direct Simulation Monte Carlo Code, named PDSC, is used to simulate the subsonic flow pasts a 2D vertical plate for studying the vortex-shedding phenomena. An unsteady time-averaging with temporal variable time step sampling method, called TVTS. Parametric studies, including temporal variable time step (TVTS) factor, particles per cell, number of temporal nodes, domain size and Reynolds number, are conducted to obtain the Strouhal number and aerodynamics coefficients. The results are compared to experimental data in the continuum region and simulations from the literature wherever they are available. Results of TVTS=100 and 150 has oscillation phenomenon, but results of TVTS=100 has results clear vortex shedding. Both the Strouhal number (0.174, 0.188, and 0.21) and the average drag coefficients (1.05, 1.14, 1.35, and 1.4) are increased with respect to Re=73, 126, 287 and 412 respectively, expect that the Strouhal

(6)

(7)

List of Contents

致謝...I 摘要... II Abstract ... III List of Contents ... V List of Tables ...VII List of Figures... VIII Nomenclature ... XVII

Chapter 1 Introduction ... 1

1. 1 Motivation and Background... 1

1.1. 1 Importance of Flow past a Vertical Flat Plate ... 1

1.1. 2 Classification of Flow Rarefaction ... 2

1.1.3 Direct Simulation Monto Carlo Method ... 3

1. 2 Literature Survey ... 4

1. 3 Specific Objectives of the Thesis ... 5

Chapter 2 Numerical Method ... 7

2. 1 The Boltzmann Equation ... 7

2. 2 General Description of the Standard DSMC... 8

2. 3 General Description of the PDSC ... 13

2. 4 General Description of Unsteady Sampling Method in DSMC [JCP paper in March 2008] ... 14

2. 5 DSMC Rapid Ensemble Averaging Method (DREAM) [JCP paper in March 2008]……….16

Chapter 3 Results and Discussion... 17

3. 1 Problem Description and Test Conditions ... 17

3.1 1 Test Flow past a Vertical Flat Plate with Different TVTS Factors ... 17

3.1 2 Test Flow past a Vertical Flat Plate with Different Particles per cell ... 18

3.1 3 Test Flow past a Vertical Flat Plate with Different Number of Temporal Nodes ... 19

3.1 4 Test Flow past a Vertical Flat Plate with Different Domains Sizes ... 20

3.1 5 Test Flow past a Vertical Flat Plate with Different Reynolds Numbers ... 20

3. 2 Effects of TVTS Factor... 21

3.2 1 General Simulation Results ... 23

3.2 2 Property Distributions of Vertical Flat Plate ... 24

3.2 3 Stagnation Point From Flow past a Vertical Flat Plate ... 25

3. 3 Effects of Particle per cell... 25

3.3 1 General Simulation Results ... 26

(8)

3.3 3 Stagnation Point From Flow past a Vertical Flat Plate ... 28

3. 4 Effects of Number of Temporal Node ... 28

3.4.1 General Simulation Results ... 29

3.4.2 Property Distributions of Vertical Flat Plate ... 29

3.4.3 Stagnation Point From Flow past a Vertical Flat plate... 30

3. 5 Effects of Domain Size ... 30

3.5.3 Stagnation Point From Flow past a Vertical Flat Plate ... 33

3. 6 Effects of Reynolds Number... 33

3.6.3 Stagnation Point From Flow past a Vertical Flat Plate ... 35

3. 7 Effects of Knudsen Number... 36

Chapter 4 Conclusions and Recommendation of Future Work ... 36

4. 1 Summary ... 36

4. 2 Recommendation of Future Work ... 37

(9)

List of Tables

Table 1 Flow past a vertical plate with different TVTS factors ... 42

Table 2 Flow past a vertical plate with different particles per cell ... 43

Table 3 Flow past a vertical plate with different number of temporal nodes ... 44

Table 4 Flow past a vertical plate with different domain sizes ... 45

(10)

List of Figures

Figure 2. 1 Classifications of Gas Flows... 47

Figure 2. 2 Flow chart of the DSMC method. ... 48

Figure 2. 3 Simplified flow chart of the parallel DSMC method for np processors... 49

Figure 2. 4 The additional schemes in the parallel DSMC code. ... 50

Figure 2. 5 Sampling method in DSMC include (a) steady sampling (b) unsteady ensemble sampling (c) unsteady time averaging.(d) unsteady time averaging with temporal variable time step (TVTS)... 51

Figure 2. 6 Simplified flow chart of the unsteady parallel DSMC method. ... 52

Figure 2. 7 Simplified flow chart of the DSMC Rapid Ensemble Averaging Method (DREAM)... 53

Figure 3.1 Computational domains for the developing vertical flat plat flow. (a) H=3L; (b) H=5L; (c) H=7L ... 54

Figure 3. 2 The mesh for Re=126 and Kn=0.01 vertical flat plat flow. (a) 500 by 150 (H=3L); (b) 500 by 250 (H=5L); (c) 500 by 350 (H=7L) (∆x = ∆y~ 2λλλ) ... 55 λ Figure 3. 3 The mesh for Re=73, Kn=0.017 and Mesh=500 by 250 vertical flat plat. (∆x = ∆y~ λλλλ) ... 56

Figure 3. 4 The mesh for Re=287, Kn=0.0044 and Mesh=1000 by 500 vertical flat plat. (∆x = ∆y~ 2λλλλ) ... 56

Figure 3. 5 The mesh for Re=412, Kn=0.0031 and Mesh=1000 by 500 vertical flat plat. ( ∆x = ∆y~ 3λλλλ )... 56

Figure 3. 6 Contours of U-velocity at different instant times for the 2D vertical flat plate vortex-shedding problem. (Case 1 at Table 1, TVTS factor = 100) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 57

Figure 3. 10 Contours of U-velocity at different instant times for the 2D vertical flat plate vortex-shedding problem. (Case 5 at Table 1, TVTS factor = 300) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 61 Figure 3. 11 Contours of V-velocity at different instant times for the 2D vertical flat plate

(11)

vortex-shedding problem. (Case 1 at Table 1, TVTS factor = 100) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 62 Figure 3. 12 Contours of V-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 5 at Table 1, TVTS factor = 300) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 66 Figure 3. 16 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 1 at Table 1, TVTS factor = 100) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 67 Figure 3. 17 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 5 at Table 1, TVTS factor = 300) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 71 Figure 3. 21 Sketch of the vortex shedding flow after a vertical flat plat flow. (a) H=3L; (b)

H=5L; (c) H=7L ... 72 Figure 3. 22 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in TVTS factor=100. (Case 1 at Table 1) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 73 Figure 3. 23 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in TVTS factor=150. (Case 2 at Table 1) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01;(c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 74

(12)

Figure 3. 24 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate vortex- shedding problem in TVTS factor=198. (Case 3 at Table 1) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 75 Figure 3. 25 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in TVTS factor=220. (Case 4 at Table 1) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 76 Figure 3. 26 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in TVTS factor=300. (Case 5 at Table 1) (a) x=0.03,y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0... 77 Figure 3. 27 Time trace of Drag Coefficient distributions for 2D vertical flat plate vortex-

shedding problem. (a) TVTS factor=100 (Case 1 at Table 1); (b) TVTS factor=150 (Case 2 at Table 1); (c) TVTS factor=198 (Case 3 at Table 1); (d) TVTS factor=220 (Case 4 at Table 1); (e) TVTS factor=300 (Case 5 at Table 1)

... 79 Figure 3. 28 Time trace of Pressure distributions for 2D vertical flat plate vortex-shedding

problem. (a) TVTS factor=100 (Case 1 at Table 1); (b) TVTS factor=150 (Case 2 at Table 1); (c) TVTS factor=198 (Case 3 at Table 1); (d) TVTS factor=220 (Case 4 at Table 1); (e) TVTS factor=300 (Case 5 at Table 1)... 80 Figure 3. 29 The stagnation point for flow past a vertical flat plate with different TVTS

factors at normalized time. (Table 1)... 81 Figure 3. 30 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 1 at Table 2, Particle per cell = 50) (a) 74.8µs; (b) 972.4 µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 82 Figure 3. 31 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 2 at Table 2, Particle per cell = 100) (a) 74.8µs; (b) 972.4 µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs

... 83 Figure 3. 32 Contours of U-velocity at different instant times for the 2D vertical flat plate

... 84 Figure 3. 33 Contours of V-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 1 at Table 2, Particle per cell = 50) (a) 74.8µs; (b) 972.4 µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 85 Figure 3. 34 Contours of V-velocity at different instant times for the 2D vertical flat plate

(13)

74.8µs; (b) 972.4 µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 86 Figure 3. 35 Contours of V-velocity at different instant times for the 2D vertical flat plate

... 87 Figure 3. 36 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 1 at Table 2, Particle per cell = 50) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 88 Figure 3. 37 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 2 at Table 2, Particle per cell = 100) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 89 Figure 3. 38 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 3 at Table 2, Particle per cell = 200) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs ... 90 Figure 3. 39 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in Particle per cell = 50. (Case 1 at Table 2) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0... 91 Figure 3. 40 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in Particle per cell = 100. (Case 2 at Table 2) (a) x=0.03,y=0.01; (b)x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0... 92 Figure 3. 41 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in Particle per cell = 200. (Case 3 at Table 2) (a) x=0.03,y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0... 93 Figure 3. 42 Time trace of Drag Coefficient distributions for 2D vertical flat plate vortex-

shedding problem. (a) Particle per cell = 50 (Case 1 at Table 2); (b) Particle per cell = 100 (Case 2 at Table 2); (c) Particle per cell = 200 (Case 3 at Table 2). 94 Figure 3. 43 Time trace of Pressure distributions for 2D vertical flat plate vortex-shedding

problem. (a) Particle per cell = 50 (Case 1 at Table 2); (b) Particle per cell = 100 (Case 2 at Table 2); (c) Particle per cell = 200 (Case 3 at Table 2) ... 95 Figure 3. 44 The stagnation point for flow past a vertical flat plate with different particles

per cell at normalized time. (Table 2)... 96 Figure 3. 45 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 1 at Table 3, Number of temporal node = 140) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f)

(14)

8751.6µs ... 97 Figure 3. 46 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 2 at Table 3, Number of temporal node = 280) (a) 37.4µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f)

vortex-shedding problem. (Case 3 at Table 3, Number of temporal node = 560) (a) 18.7µs; (b) 729.3µs; (c) 7349.1µs; (d) 7461.3µs; (e) 7554.8µs; (f)

7648.3µs ... 99 Figure 3. 48 Contours of V-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 1 at Table 3, Number of temporal node = 140) (a) 74.8µs; (b) 972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f)

vortex-shedding problem. (Case 2 at Table 3, Number of temporal node = 280) (a) 37.4µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f)

vortex-shedding problem. (Case 3 at Table 3, Number of temporal node = 560) (a) 18.7µs; (b) 729.3µs; (c) 7349.1µs; (d) 7461.3µs;(e) 7554.8µs; (f)

7648.3µs ... 102 Figure 3. 51 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 1 at Table 3, Number of temporal node = 140) (a) 74.8µs; (b)972.4µs; (c) 8527.2µs; (d) 8602µs; (e) 8676.8µs; (f) 8751.6µs

shedding problem. (Case 2 at Table 3, Number of temporal node= 280) (a) 74.8µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs

shedding problem. (Case 3 at Table 3, Number of temporal node = 560) (a) 18.7µs; (b) 729.3µs; (c) 7349.1µs; (d) 7461.3µs; (e) 7554.8µs; (f)

7648.3µs ... 105 Figure 3. 54 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in Number of temporal node = 140. (Case 1 at Table 3) (a) x=0.03, y=0.01; (b) x=0.06, y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 106 Figure 3. 55 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

(15)

Table 3) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 107 Figure 3. 56 Time traces of stream-wise U-velocity and V-velocity for 2D vertical flat plate

vortex- shedding problem in Number of temporal node = 560. (Case 3 at Table 3) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (c) x=0.09, y=0.01; (d) x=0.03, y=0; (e) x=0.06, y=0; (f) x=0.09, y=0 ... 108 Figure 3. 57 Time trace of Drag Coefficient distributions for 2D vertical flat plate vortex-

shedding problem. (a) Number of temporal node = 140 (Case 1 at Table 3); (b) Number of temporal node = 280 (Case 2 at Table 3); (c) Number of temporal node = 560 (Case 3 at Table 3) ... 109 Figure 3. 58 Time trace of Pressure distributions for 2D vertical flat plate vortex-shedding

problem. (a) Number of temporal node = 140 (Case 1 at Table 3); (b) Number of temporal node = 280 (Case 2 at Table 3); (c) Number of temporal node = 560 (Case 3 at Table 3)... 110 Figure 3. 59 The stagnation point for different Number of temporal nodes at normalized

time. (Table 3)...111 Figure 3. 60 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 1 at Table 4, Cell Number = 500 by 150 (H=3L)) (a) 37.4 µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs ... 112 Figure 3. 61 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 2 at Table 4, Cell Number = 500 by 250 (H=5L)) (a) 37.4 µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs ... 113 Figure 3. 62 Contours of U-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 3 at Table 4, Cell Number = 500 by 350 (H=7L)) (a) 37.4 µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs ... 114 Figure 3. 63 Contours of V-velocity at different instant times for the 2D vertical flat plate

(16)

(H=7L)) (a) 37.4 µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs ... 117 Figure 3. 66 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 1 at Table 4, Cell Number = 500 by 150 (H=3L)) (a) 37.4µs; (b) 972.4µs; (c) 8602µs; (d) 8714.2µs; (e) 8826.4µs; (f) 8938.6µs

... 120 Figure 3. 69 Time traces of stream-wise U-velocity and -Velocity for 2D vertical flat plate

vortex- shedding problem in Cell Number = 500 by 150. (Case 1 at Table 4, H=3L) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06,y=0; (f) x=0.09,y=0 ... 121 Figure 3. 70 Time traces of stream-wise U-velocity and V-Velocity for 2D vertical flat plate

vortex- shedding problem in Cell Number = 500 by 250. (Case 2 at Table 4, H=5L) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (c)

x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06,y=0; (f) x=0.09,y=0... 122 Figure 3. 71 Time traces of stream-wise U-velocity and V-Velocity for 2D vertical flat plate

vortex- shedding problem in Cell Number = 500 by 350. (Case 3 at Table 4, H=7L) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (c)

x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06,y=0; (f) x=0.09,y=0... 123 Figure 3. 72 Time trace of Drag Coefficient distributions for 2D vertical flat plate vortex-

shedding problem. (a) Cell Number = 500 by 150 (Case 1 at Table 4, H=3L); (b) Cell Number = 500 by 250 (Case 2 at Table 4, H=5L); (c) Cell Number = 500 by 350 (Case 3 at Table 4, H=7L)... 124 Figure 3. 73 Time trace of Pressure distributions for 2D vertical flat plate vortex-shedding

problem. (a) Cell Number = 500 by 150 (Case 1 at Table 4, H=3L); (b) Cell Number = 500 by 250 (Case 2 at Table 4, H=5L); (c) Cell Number = 500 by 350 (Case 3 at Table 4, H=7L) ... 125 Figure 3. 74 The stagnation point for different domain sizes at normalized time. (Table 4)

... 126 Figure 3. 75 Contours of U-velocity at different instant times for the 2D vertical flat plate

(17)

63.5µs; (b) 1651.0 µs; (c) 9080.5µs; (d) 9207.5µs; (e) 9334.5µs; (f)

vortex-shedding problem. (Case 2 at Table 5, Reynolds number=126) (a) 18.7µs; (b) 729.3µs; (c) 7349.1µs; (d) 7461.3µs; (e) 7554.8µs; (f)

vortex-shedding problem. (Case 3 at Table 5, Reynolds number=287) (a) 18.7µs; (b) 878.9 µs; (c) 4301.0µs; (d) 4357.1µs; (e) 4413.2µs ; (f)

vortex-shedding problem. (Case 4 at Table 5, Reynolds number=412) (a) 11.25µs; (b) 855.0µs; (c) 2452.5µs; (d) 2520.0µs; (e) 2587.5µs; (f)

vortex-shedding problem. (Case 1 at Table 5, Reynolds number=73) (a) 63.5µs; (b) 1651.0 µs; (c) 9080.5µs; (d) 9207.5µs; (e) 9334.5µs; (f)

vortex-shedding problem. (Case 2 at Table 5, Reynolds number=126) (a) 18.7µs; (b)729.3µs; (c) 7349.1µs; (d) 7461.3µs ;(e) 7554.8µs; (f)

vortex-shedding problem. (Case 3 at Table 5, Reynolds number=287) (a)18.7µs; (b) 878.9 µs; (c) 4301.0µs; (d) 4357.1µs; (e) 4413.2µs; (f) 4469.3µs ... 133 Figure 3. 82 Contours of V-velocity at different instant times for the 2D vertical flat plate

vortex-shedding problem. (Case 4 at Table 5, Reynolds number=412) (a) 11.25µs; (b) 855.0 µs; (c) 2452.5µs; (d) 2520.0µs; (e) 2587.5µs; (f)

2655.0µs ... 134 Figure 3. 83 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 1 at Table 5, Reynolds number=73) (a) 63.5µs; (b) 1651.0µs; (c) 9080.5µs; (d) 9207.5µs; (e) 9334.5µs; (f) 9461.5µs... 135 Figure 3. 84 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 2 at Table 5, Reynolds number=126) (a) 18.7µs; (b) 729.3µs; (c) 7349.1µs; (d) 7461.3µs; (e) 7554.8µs; (f) 7648.3µs... 136 Figure 3. 85 Streamline at different instant times for the 2D vertical flat plate vortex-

(18)

878.9µs; (c) 4301.0µs (d) 4357.1µs; (e) 4413.2µs; (f) 4469.3µs... 137 Figure 3. 86 Streamline at different instant times for the 2D vertical flat plate vortex-

shedding problem. (Case 4 at Table 5, Reynolds number=412) (a) 11.25µs; (b) 855.0µs; (c) 2452.5µs; (d) 2520.0µs; (e) 2587.5µs; (f) 2655.0µs... 138 Figure 3. 87 Time traces of stream-wise U-velocity for 2D vertical flat plate vortex-

shedding problem in Reynolds number=73. (Case 1 at Table 5) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06,y=0; (f) x=0.09,y=0 ... 139 Figure 3. 88 Time traces of stream-wise U-velocity for 2D vertical flat plate vortex-

shedding problem in Reynolds number=412. (Case 4 at Table 5) (a) x=0.03, y=0.01; (b) x=0.06,y=0.01; (c) x=0.09,y=0.01; (d) x=0.03, y=0; (e) x=0.06,y=0; (f) x=0.09,y=0 ... 142 Figure 3. 91 Time trace of Drag Coefficient distributions for 2D vertical flat plate vortex-

shedding problem. (a) Reynolds number=73 (Case 1 at Table 5); (b) Reynolds number=126 (Case 2 at Table 5); (c) Reynolds number=287 (Case 3 at Table 5); (d) Reynolds number=412 (Case 4 at Table 5) ... 143 Figure 3. 92 Time trace of Pressure distributions for 2D vertical flat plate vortex-shedding

problem. (a) Reynolds number=73 (Case 1 at Table 5); (b) Reynolds

number=126 (Case 2 at Table 5); (c) Reynolds number=287 (Case 3 at Table 5); (d) Reynolds number=412 (Case 4 at Table 5) ... 144 Figure 3. 93 The Stagnation point from flow past a vertical plate at normalized time. (a)

Reynolds number=73 (Case 1 at Table 5); (b) Reynolds number=126 (Case 2 at Table 5); (c) Reynolds number=287 (Case 3 at Table 5); (d) Reynolds

number=412 (Case 4 at Table 5) ... 145 Figure 3. 94 Strouhal number variation as a function of Reynolds number. (Table 5) ... 146 Figure 3. 95 Strouhal number variation as a function of Knudsen number. (Table 5) ... 146 Figure 3. 96 Dimensionless critical time at which the symmetrical twin vortices to become

(19)

Nomenclature

λ ：mean free path

ρ ：density

σ ：the differential cross section ω ：viscosity temperature exponent Ω : space domain rot ε ：rotational energy v ε ：vibrational energy rot

ζ ：rotational degree of freedom

v

ζ ：vibrational degree of freedom t ：time-step

σT ：the total cross section

c ：the total velocity

c’ ：random velocity co ：mean velocity cr ：relative speed d ：molecular diameter

D ：the throat diameter of the twin-jet interaction

dref ：reference diameter

E ：energy

k ：the Boltzmann constant Kn ：Knudsen number

Knmax : continuum breakdown parameter

. max

Thr

Kn : the threshold value of continuum breakdown parameter

Q

Kn : local Knudsen numbers based on flow property Q L ：characteristic length;

m ：molecule mass

M∞ ：free-stream mach number n ：number density

Tne

(20)

.

Thr Tne

P : the threshold value of thermal non-equilibrium indicator Re ：Reynolds number

To ：stagnation temperature T∞ ：free-stream temperature

Tref ：reference temperature

Trot ：rotational temperature Ttot ：total temperature

Ttr ：translational temperature Tv ：vibrational temperature Tw ：wall temperature S t ：Strouhal number D C ：Drag coefficient

f ：half plate shedding frequency

t ：simulation flow time

x ：the length of the wake bubble

Re ：Reynolds number

U ：U-velocity

(21)

Chapter 1 Introduction

1. 1

Motivation and Background

1.1. 1 Importance of Flow past a Vertical Flat Plate

The phenomena of vortex shedding associated with the subsonic external flow problems in different length scales are visible everywhere in fluid dynamics. The periodic two -dimensional vortex shedding behind bluff bodies exposed to uniform flow has fascinated researchers since the days of Leonardo da Vinci. The occurrence of this flow phenomenon is due to instabilities and depends on the geometry of the bluff body and the Reynolds number; it has often been the cause of failure of flow-exposed structures in various fields of engineering.

The motion of a viscous incompressible fluid past an object results in the generation of a wake containing a double row of vortices, the so-called von Kármán [T. von Kármán, 1956]. The behavior of such flows may be characterized through a single parameter, the Reynolds number. From experimental observations it is held that for small values of Re less than some threshold value Reynolds number the flow is steady, whilst when Reynolds number is exceeded the flow becomes unsteady and vortices are shed, leading eventually to regular periodic shedding. As Re increases still further, eventually a turbulent regime is encountered and the vortex shedding pattern then assumes an irregular structure.

At Re < 40 a stationary, symmetric vertical pattern is formed behind a buff body, and at higher Reynolds a periodic wake is formed. Wake formation is described in detail by Younis [1988] and Roshko [1954]. As Reynolds the flow separates from the plate ends creating free shear layers that continue downstream before rolling into vortices4. These disturbances develop increasingly closer to the plate.

(22)

points on either side of the plate, corresponding to the plate tips. Thus, a vortex will be generated in the region behind the separation point on one side, while a corresponding vortex on the other side will break away from the cylinder and move downstream in the wake. When the attached vortex reaches a particular strength, it will in turn break away and a new vortex will begin to develop again on the second side and so on.

A unique relationship is found to exist between Re and the dimensionless Strouhal number, St. Roshko [1954] noted that the bluffness of a body is related to the ratio of wake width to body size. A bluffer body diverge the flow more extensively producing a wider wake and a lower shedding frequency. Roshko [1954] determined that St for flow over a flat plate remained nearly constant, 0.13≤ St ≤ 0.140, over a wind rang of Re.

1.1. 2 Classification of Flow Rarefaction

Knudsen number (Kn=λ/L) is usually used to indicate the degree of rarefaction. Note that the mean free path λ is the average distance traveled by molecules before collision and L is the flow characteristic length. In general, flows are divided into four regimes and three solutions. When the local Kn number approaches zero, the flow reaches inviscid limit and can be solved by Euler equation. As local Kn increases, molecular nature of the gas becomes dominated. Hence, when the flow is close to the continuum regime (Kn approach 0.01), the well known Navier-Stokes equation may be applied to yield accurate result for engineering purposes. For Kn larger than 0.01, continuum assumption begins to break down and the particle-based method is necessary and a kinetic approach, based on the Boltzmann equation [Cercignani, 1998]. It is important to note that the kinetic approach is valid in the whole range of the gas rarefaction. This is an important advantage when systems with multiscale physics are investigated; however it is 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. The well known

(23)

direct simulation Monte Carlo (DSMC) method [Bird, 1994] is also a powerful computational scheme.

1.1.3 Direct Simulation Monto Carlo Method

Direct Simulation Monte Carlo (DSMC), was proposed by Bird to solve the Boltzmann equation using direct simulation of particle collision kinetics, and the associated monograph was published in 1994 [Bird’s book]. Later on, both Nanbu [1986] and Wagner [1992] were able to demonstrate mathematically that the DSMC method is equivalent to solving the Boltzmann equation as the simulated number of particles becomes large. The DSMC method is a particle method for the simulation of 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. Mass, 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. This method has become a widely used computational tool for the simulation of gas flows in the transitional flow regime, in which molecular effects become important.

The DSMC method becomes very time-consuming as the flow approaches continuum regime since the sampling cell size has to be much smaller than the local mean free path for the solution to be accurate. Several remedies in speeding up the DSMC computation include: (1) parallel computing [Robinson, 1996-1998]; (2) variable time-step scheme for steady flows [Kannenberg, 2000]; (3) sub-cells within each sampling cell [Bird’s book, 1994], and (4) unsteady flows sampling. Details of the “parallel computing”, “variable time-step scheme” and “sub-cell” can be found in those references cited in the above and are not described here

(24)

for brevity. Only “unsteady flows sampling” concept is described here since it was rarely discussed in detail in the literature. Unsteady flow simulations have difficulties in DSMC sampling, Cave, et al. [2008] developed new model for under-expanded jets from sonic nozzles during the start up of rocket nozzles and during the injection phase of the Pulsed Pressure Chemical Vapor Deposition (PP-CVD) process. But sampling over a small time interval requires either a very large number of simulated molecules or the average of a large number of separate simulations (“ensemble-averaging”). The costs high computational expense and large memory. Xu [1993] used one method of decreasing the statistical scatter of the results is to using statistical smoothing procedures in two dimensional unsteady problems. Wagner [1992] proved Bird’s two-dimensional axis-symmetric code which incorporates unsteady sampling techniques in which a number of time intervals close to the sampling point are averaged (“time-averaging”); however this is single processor code. The increased computational capacities of parallel-DSMC techniques have the potential to enable the simulation of time-dependent flow problems at the near-continuum regime.

Accordingly, in this thesis develops an unsteady time-averaging sampling method [Cave,

et al., 2008] for flow past a vertical plate and uses DSMC rapid ensemble averaging method (DREAM) to reduce the statistical scatter with an acceptable runtime for unsteady flow simulation.

1. 2 Literature Survey

Flat plate flows are encountered in systems not in equilibrium. Proto-type flows of this kind are circular cylinder problem and flat plate flow problem in two or three dimension. The circular cylinder problem has been studies [Blasius, 1908]. Recently remarks advance in numerical integration of the Navier-Stokes equations has been made owing to the tremendous of development high-speed computers, and time-dependent solutions for impulsively started flow past a circular cylinder. G.A.Bired et al [1997] calculated the two dimensional flow over

(25)

a vertical plate at Knudsen number of 0.0025 and 0.005 for Mach number of 0.5 using the direct simulation Monte Carlo method. S.Taneda, et al [1971] was investigated the separated flow past a flat plate experimentally. Their measurement were conducted using flow-visualization techniques for the length of wake bubble at Reynolds numbers ranging from 18.1 to 1135.Taneda, et al [1963] was investigated experimentally at Reynolds number ranging from -1

10 to₁₀-5_{.They measured import points resulting from experiments are that the}

drag and life coefficients.

Cave et al., [2008] developed a unsteady sampling procedures for the parallel direction simulation Monte Carlo method. This paper is simulations of a shock tube and the development of Couette flow are then carried out as validation studies. To overcome the large computational expense and memory requirements usually involved in DSMC simulations of unsteady flows. “Time-averaging” method was implemented which has considerable computational advantages over ensemble-averaging a large number of separate runs. Also, in order to reduce the statistical scatter with an acceptable runtime for unsteady flow simulation using DSMC technique. DSMC-DREAM (Rapid Ensembled Average Method) was developed whereby a combination of time- and ensemble-averaged data was build up by regenerating the particle data a short time prior to the sampling point of interesting, assuming a Maxwell-Boltzmann distribution of particle velocities. In this thesis, we used this method and validated this code.

1. 3 Specific Objectives of the Thesis

The current objectives of the thesis are summarized as follows:

1. To simulation of subsonic (Ma=0.77) flow over a two-dimensional vertical flat plate and vortex shedding behind a vertical flat plate using parallelized DSMC code (PDSC). 2. To simulation vertical flat plate, including same vertical flat plate with height of 0.02

(26)

3. To simulation vertical flat plate using an unsteady time-averaging with temporal variable time step sampling method.

4. To discuss the effects of TVTS (Unsteady time average with temporal variable time step) factor, Particle per cell, Number of temporal node, Domain size and Reynolds Number.

The organization of the thesis is stated as follows: Chapter 1 describes the Introduction, Chapter 2 describes the Numerical Method, Chapter 3 describes Simulation Results and finally Chapter 4 recommendation of future work.

(27)

Chapter 2 Numerical Method

2. 1 The Boltzmann Equation

The Knudsen number (Kn) is used to indicate the degree of rarefaction. In Figure 2.1, flows are divided into four regimes and three solutions. We have found the Boltzmann equation is valid for all flow regimes which from 10 to 0.0001. It is one of the most important transport equations in non-equilibrium statistical mechanics, which deals with systems far from thermodynamics equilibrium. There are some assumptions made in the derivation of the Boltzmann equation which defines limits of applicability. They are summarized as follows:

1. Molecular chaos is assumed which is valid when the intermolecular forces are short range. It allows the representation of the two particles distribution function as a product of the two single particle distribution functions.

2. Distribution functions do not change before particle collision. This implies that the encounter is of short time duration in comparison to the mean free collision time. 3. All collisions are binary collisions.

4. Particles are uninfluenced by intermolecular potentials external to an interaction. According to these assumptions, the Boltzmann equation is derived and shown as Eq. (2.1)

4 2 ' ' 1 1 0 ( ) ( ) ( ) ( ) i i c i i i nf nf nf f u F n f f ff g d dU t x u x π σ ∞ −∞ ∂ ∂ ∂ ∂ + + = = − Ω ∂ ∂ ∂ ∂

∫ ∫

(2.1)

Meaning of particle phase-space distribution function f is the number of particles with center of mass located within a small volume d r3 near the point r , and velocity within a ranged u3 , at time t . F_iis an external force per unit mass and t is the time and u_iis the molecular velocity. σ is the differential cross section and dΩ is an element of solid angle. The prime denotes the post-collision quantities and the subscript 1 denotes the collision partner. Meaning of each term in Eq. (2.1) is described in the following;

(28)

1. The first term on the left hand side of the equation represents the time variation of the distribution function of the particles (unsteady term).

2. The second term gives the spatial variation of the distribution function (flux term). 3. The third term describes the effect of a force on the particles (force term).

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

In general, it is difficult to solve the Boltzmann equation directly using numerical method because the difficulties of correctly modeling the integral collision term. Instead, the DSMC method was used to simulated problems involving rarefied gas dynamics, which is the main topic in the current thesis.

2. 2 General Description of the Standard DSMC

In order to the expected rarefaction caused by the rarefied gas flows, the direct simulation Monte Carlo (DSMC) method which is a particle-based method developed by Bird during the 1960s and is widely used an efficient technique to simulate rarefied gas regime [Bird, 1976 and Bird, 1994]. In the DSMC method, a large number of particles are generated in the flow field to represent real physical molecules rather than a mathematical foundation and it has been proved that the DSMC method is equivalent to solving the Boltzmann equation [Nanbu, 1986 and Wagner, 1992]. The assumptions of molecular chaos and a dilute gas are required by both the Boltzmann formulation and the DSMC method [Bird, 1976 and Bird, 1994]. An important feature of DSMC is that the molecular motion and the intermolecular collisions are uncoupled over the time intervals that are much smaller than the mean collision time. Both the collision between molecules and the interaction between molecules and solid boundaries are computed on a probabilistic basis and, hence, this method makes extensive use of random numbers. In most practical applications, the number of

(29)

simulated molecules is extremely small compared with the number of real molecules. The general procedures of the DSMC method are described in the next section, and the consequences of the computational approximations can be found in Bird [Bird, 1976 and Bird, 1994].

In DSMC, there are three molecular collision models for real physical behavior and imitate the real particle collision, which are the Hard Sphere (HS), Variable Hard Sphere (VHS) and Variable Soft Sphere (VSS) molecular models, in the standard DSMC method [Bird, 1994]. The collision pairs then are chosen by the acceptance-rejection method. The no time counter (NTC) method is an efficient method for molecular collision. This method yield the exact collision rate in both simple gases and gas mixtures, and under either equilibrium or non-equilibrium conditions.

Figure 2.2 is a general flowchart of the DSMC method. Important steps of the DSMC method include setting up the initial conditions, moving all the simulated particles, indexing (or sorting) all the particles, colliding between particles and sampling the molecules within cells to determine the macroscopic quantities. The details of each step will be described in the following;

Initialization

The first step to use the DSMC method in simulating flows is to set up the geometry and flow conditions. A physical space is discredited into a network of cells and the domain boundaries have to be assigned according to the flow conditions. An important feature has to be noted is the size of the computational cell should be smaller than the mean free path, and the distance of the molecular movement per time step should be smaller than the cell dimension. After the data of geometry and flow conditions have been read in the code, the numbers of each cell is calculated according to the free-stream number density and the current cell volume. The initial particle velocities are assigned to each particle based on the Maxwell-Boltzmann distribution according to the free-stream velocities and temperature, and

(30)

the positions of each particle are randomly allocated within the cells.

Molecular Movement

After initialization process, the molecules begin move one by one, and the molecules move in a straight line over the time step if it did not collide with solid surface. For the standard DSMC code by Bird [Bird, 1976 and Bird, 1994], the particles are moved in a structured mesh. There are two possible conditions of the particle movement. First is the particle movement without interacting with solid wall. The particle location can be easy located according to the velocity and initial locations of the particle. Second is the case that the particle collides with solid boundary. The velocity of the particle is determined by the boundary type. Then, the particle continues its journey from the intersection point on the cell surface with its new absolute velocity until it stops. Although it is easier to implement by using structured mesh, it is difficult for those flows with complex geometry.

Indexing

The location of the particle after movement with respect to the cell is important information for particle collisions. The relations between particles and cells are reordered according to the order of the number of particles and cells. Before the collision process, the collision partner will be chosen by a random method in the current cell. And the number of the collision partner can be easy determined according to this numbering system.

Gas-Phase Collisions

The other most important phase of the DSMC method is gas phase collision. The current DSMC method uses the no time counter (NTC) method to determine the correct collision rate in the collision cells. The number of collision pairs within a cell of volume V_C over a time interval ∆ is calculated by the following equation; t

c r T N c t V F N N ( ) / 2 1 max∆ σ (2.2)

(31)

is the particle weight, which is the number of real particles that a simulated particle represents.

T

σ and c_r are the cross section and the relative speed, respectively. The collision for each pair is computed with probability

max

) /( )

(σ_Tc_r σ_Tc_r (2.3)

The collision is accepted if the above value for the pair is greater than a random fraction. Each cell is treated independently and the collision partners for interactions are chosen at random, regardless of their positions within the cells. The collision process is described sequentially as follows:

1. The number of collision pairs is calculated according to the NTC method, Eq. (2.2), for each cell.

2. The first particle is chosen randomly from the list of particles within a collision cell.

3. The other collision partner is also chosen at random within the same cell.

4. The collision is accepted if the computed probability, Eq (2.3), is greater than a random number.

5. If the collision pair is accepted then the post-collision velocities are calculated using the mechanics of elastic collision. If the collision pair is not to collide, continue choosing the next collision pair.

6. If the collision pair is polyatomic gas, the translational and internal energy can be redistributed by the Larsen-Borgnakke model [Borgnakke and Larsen, 1975], which assumes in equilibrium.

(32)

and then progress to the next step.

Sampling

After the particle movement and collision process finish, the particle has updated positions and velocities. The macroscopic flow properties in each cell are assumed to be constant over the cell volume and are sampled from the microscopic properties of each particle within the cell. The macroscopic properties, including density, velocities and temperatures, are calculated in the following equations [Bird, 1976 and Bird, 1994];

nm = ρ _(2.4a) ' c c c c_o = = _o + _(2.4b) ) ' ' ' ( 2 1 2 3 2 2 2 w v u m kT_tr = + + (2.4c) ) ( 2 r rot rot k T = ε ζ (2.4d) ) ( 2 v v v k T = ε ζ (2.4e) ) 3 ( ) 3

( _tr _rot _rot _v _v _rot _v

tot T T T

T = +ζ +ζ +ζ +ζ (2.4f)

n, m are the number density and molecule mass, receptively. c, co, and c’ are the total

velocity, mean velocity, and random velocity, respectively. In addition, Ttr, Trot, Tv and Ttot are

translational, rotational, vibration and total temperature, respectively. ε_rotandε_vare the rotational and vibration energy, respectively. ζ_rot and ζ_v are the number of degree of freedom of rotation and vibration, respectively. If the simulated particle is monatomic gas, the translational temperature is regarded simply as total temperature. Vibration effect can be

neglect if the temperature of the flow is low enough. The flow will be monitored if steady state is reached. If the flow is under unsteady

(33)

As a rule of thumb, the sampling of particles starts when the number of molecules in the calculation domain becomes approximately constant.

2. 3 General Description of the PDSC

Although the large number of molecules in a real gas is replaced with a reduced number of model particles, there are still a large number of particles must be simulated, leading to tremendous computer power requirements and needing to cost a lot of computational time. As a result, parallel DSMC method is developed to solve the problem. Figure 2.3 illustrates a simplified flow chart of the parallel DSMC method used in the current study. The DSMC algorithm is readily parallelized through physical domain decomposition. The cells of the computational grid are distributed among the processors. Each processor executes the DSMC algorithm in serial for all particles and cells in its domain. Data communication occurs when particles cross the domain (processor) boundaries and are then transferred between processors.

Parallel DSMC Code (PDSC) is the main solver used in this thesis, which utilizes unstructured tetrahedral mesh. Figure 2.4 is the features of PDSC and brief introduction is listed in the following paragraphs.

1. 2D/2D-axisymmetric/3-D unstructured-grid topology: PDSC can accept either 2D/2D-axisymmetric (triangular, quadrilateral or hybrid triangular-quadrilateral) or 3D (tetrahedral, hexahedral or hybrid tetrahedral-hexahedral) mesh [Wu et al.’s JCP paper, 2006]. Computational cost of particle tracking for the unstructured mesh is generally higher than that for the structured mesh. However, the use of the unstructured mesh, which provides excellent flexibility of handling boundary conditions with complicated geometry and of parallel computing using dynamic domain decomposition based on load balancing, is highly justified.

2. Parallel computing using dynamic domain decomposition: Load balancing of PDSC is achieved by repeatedly repartitioning the computational domain using a multi-level graph-partitioning tool, PMETIS [Wu and Tseng, 2005] by taking

(34)

advantage of the unstructured mesh topology employed in the code. A decision policy for repartition with a concept of Stop-At-Rise (SAR) [Wu and Tseng, 2005] or constant period of time (fixed number of time steps) can be used to decide when to repartition the domain. Capability of repartitioning of the domain at constant or variable time interval is also provided in PDSC. Resulting parallel performance is excellent if the problem size is comparably large. Details can be found in Wu and Tseng [Wu and Tseng, 2005].

3. Spatial variable time-step scheme: PDSC employs a spatial variable time-step scheme (or equivalently a variable cell-weighting scheme), based on particle flux (mass, momentum, energy) conservation when particles pass interface between cells. This strategy can greatly reduce both the number of iterations towards the steady state, and the required number of simulated particles for an acceptable statistical uncertainty. Past experience shows this scheme is very effective when coupled with an adaptive mesh refinement technique [Wu et al.’s CPC paper, 2004].

4. Unsteady flow simulation: An unsteady sampling routine is implemented in PDSC, allowing the simulation of time-dependent flow problems in the near continuum range [JCP paper submitted in June 2007]. A post-processing procedure called DSMC Rapid Ensemble Averaging Method (DREAM) is developed to improve the statistical scatter in the results while minimizing both memory and simulation time. In addition, a temporal variable time-step (TVTS) scheme is also developed to speed up the unsteady flow simulation using PDSC. More details can be found in [JCP paper submitted in June 2007]. Details of the idea and implementation are described next.

5. Transient Sub-cells: Recently, transient sub-cells are implemented in PDSC directly on the unstructured grid, in which the nearest-neighbor collision can be enforced, whilst maintaining minimal computational overhead [JFM paper under preparation, 2007].

2. 4 General Description of Unsteady Sampling Method in DSMC [JCP

paper in March 2008]

(35)

therefore some modification is need for unsteady sampling. The unsteady sampling method has been described in detail in the paper [Cave, et al., 2008].

There are two methods for unsteady sampling, the differences illustrated in Figure 2.5 and the details will be described in the following;

7. The “Ensemble-Average” method, require multiple simulation runs. The flow flied is sampled at the suitable sampling times during the run. The sampling simulation outputs from each run are averaged over the runs. There the results are vary precise, but the method is very computational expensive. Because a large of runs is required to reduce the statistical scatter to smooth data and a large amount of memory is needed to record the sampling data for each simulation.

8. The “Time-Average” method, require one simulation run. It averages a number of time steps over an interval before the sampling time. However it suffers a potential disadvantage in that the results will be “smeared” over the time over which samples are taken. Hence the sample time must be sufficiently short to minimizes time “smearing” and yet long enough to obtain a good statistical sample. This method of time averaging has been used previously by Auld to model shock tube flow [Auld, 1992]

In PDSC, the method of time-averaging was implemented [Cave, et al., 2008]. Here a technique called the temporal variable time step (TVTS) method was used to reduce the simulation time by increasing the time step between sampling. The code has an option for the user to choose specific output flow times or for output at regular intervals. Figure 2.6 shows the flow chart of the PDSC method with the unsteady sampling procedures implemented. Here M is the output matrix for sampling interval M. Most parts of the procedure are the same as the steady simulation except the sampling data must be reset after completing each

(36)

simulation interval.

2. 5 DSMC Rapid Ensemble Averaging Method (DREAM) [JCP paper

March 2008]

Because reducing the statistical scatter greatly, in time-average data requires a very large number of simulation particles with consequent large computational times. In the thesis, we have adapted DREAM code which has been described in detail in this paper by [Cave, et al., 2008]. The illustrated in Figure 2.7

First, we select a raw data set X-n produced by PDSC n sampling intervals prior to the sampling interval of interest X. New particle data is generated from the macroscopic properties in data set X-n by assuming a Maxwellian distribution of velocities. The standard PDSC algorithm is then used to simulate forward in time until the sampling period of interest

X is reached. The time steps close to the sampling point are time-averaged in the same way as in PDSC and this process is repeated a number of times, thus building up a combination of ensemble-averaged and time-averaged data without having to simulate from zero flow time for each run. This process reduces the statistical scatter in the results by adding to the number of particles in the sample, rather than by some artificial smoothing process. Because only a short period of the flow is processed in this way, the scheme has significant memory and computational advantages over both ensemble-averaging and using a greater number of particles in the time-averaging scheme.

(37)

Chapter 3 Results and Discussion

A large number of experimental and numerical studies have been reported on the vortex-shedding flows in the continuum limit, while there have been very few studies focusing on similar flows in the rarefied gas regime. Major obstacle of the investigation in rarefied regime mostly came from the difficulties of experiments and also numerical simulations for unsteady flows in this regime. A general-purpose Parallel Direct Simulation Monte Carlo Code, named PDSC, is used to simulate the subsonic flow pasts a 2D vertical plate for studying the vortex-shedding phenomena.

3. 1 Problem Description and Test Conditions

To demonstrate the capability of the unsteady function in PDSC, a number of simulations of the flow past a vertical plate were conducted. A systematic study, including the effects of TVTS (Unsteady time average with temporal variable time step) factor, particle per cell, number of temporal node, domain size and Reynolds number, was first undertaken. Only one parameter was varied at a time for a set of simulations to determine its effect on the flow solution while the other simulation parameters remained unchanged from the control case.

3.1 1 Test Flow past a Vertical Flat Plate with Different TVTS Factors

Figure 3.1 (b) shows the computational domains for the developing vertical flat plat flow. Flow and simulation conditions are summarized in Table 1. The control parametric study is about the number of simulation TVTS (Unsteady time average with temporal variable time step) factor. Five cases, such as 100, 150, 198, 220 and 300 simulated TVTS factor are proposed here. The calculation employed the hard sphere molecular model for which Reynolds number is related Mach number and Knudsen number by Equation give this all case. At time t=0, Mach number =0.77 of flow past a vertical flat plate with a height of 50

以蒙地卡羅法模擬次音速流通過垂直平板的現象

國 立 交 通 大 學

機

械

工

程

機

械

工

程

機

械

工

程

機

械

工

程

學

學

學

學

系

系

系

系

DSMC Simulation of the Subsonic Flow Past a Vertical

Plate

以蒙地卡羅法模擬次音速流通過垂直平板的現象

研

研

研

研 究

究

究 生

究

生

生

生：

：

：

：林士傑

林士傑

林士傑

林士傑

指導教授

指導教授

指導教授

指導教授：

：

：吳宗信

：

吳宗信

吳宗信 教授

吳宗信

教授

教授

教授

中華民國九十

中華民國九十

中華民國九十

中華民國九十七

七

七

七年

年

年

年七

七

七月

七

月

月

月

以蒙地卡羅法模擬次音速流通過垂直平板的現象

DSMC Simulation of the Subsonic Flow Past a Vertical Plate

研 究 生：林士傑

Student: Shin-Chieh Lin

指導教授：吳宗信 博士 Advisor: Dr. Jong-Shinn Wu

國立交通大學

國立交通大學

研究

究生

吳宗信教授

研究生：林士傑

指導教授：吳宗信博士 Advisor: Dr. Jong-Shinn Wu

碩士論文