中文摘要 本研究以有限體積求解穩態可壓縮

中文摘要

本研究以有限體積求解穩態可壓縮 Navier-Stokes 方程式應用於正方形孔穴內自 然與混合對流之求解。其中 Navier-Stoke 方程式之流通量採用二階中央差分方式計 算 ， 速 度 與 溫 度 梯 度 則 採 用 四 階 等 向 元 素 差 分 方 式 ， 加 總 之 淨 流 通 量 以 四 階 Runge-Kutta 方式配合預置矩陣進行時間積分，以收斂至穩態流動。文章首先利用文 獻所提供之標準解，驗證剪切流場、自然對流等流場之解答，接著進行不同馬赫數與 不同溫差下正方形孔穴內密度變化造成相較不可壓縮流場之差異，也針對不同傾斜角 (30

^o

-150

^o

)探討重力方向與波形壁波數(0-3)對自然對流之影響。最後則進行剪切流與 浮力作用下混合對流之分析，以探討波形壁波數於混合對流之影響。

關鍵字：雷諾數、瑞利數、葛拉斯赫夫數、馬赫數、波形壁、有限體積法。

Abstract

A numerical study of natural and mixed convection inside a cavity based on a finite volume scheme for the steady compressible Navier-Stokes equations was carried out. The numerical fluxes are obtained based a second order central scheme while the gradients of velocity and temperature are calculated from a fourth order isoparametric element. The overall residual is integrated using the fourth-order Runge-Kutta scheme with a preconditioning matrix. The numerical code is validated by comparing our results with previously published results of shear and thermal driven flows. The compressible effect of Mach number of lid and the temperature difference of the differential heated vertical wall problems was examined and compared with the solutions based on incompressible flow.

For the natural convection inside an inclined cavity with a hot wavy wall problem, the results are presented for different angle of inclination (30

^o

-150

^o

) for four different undulations (0-3). Finally, the current scheme is conducted to analyze mixed convection heat transfer in lid-driven cavity with a sinusoidal wavy bottom surface.

Keywords: Reynolds number、Rayleigh number、Grashof number、Mach number、

Wavy-wall、Finite volume method.

致謝

本論文得以順利完成，首要感謝指導教授楊一龍博士，在兩年半的論文指導期間 不辭辛苦的用心教導，給予學生學習與磨練的機會，對於學術研究及待人處事方面均 獲益良多，此外承蒙口試委員鄭藏勝博士，洪哲文博士及傅武雄博士的指導，並提出 許多寶貴的建議，使得本論文內容更加充實完備，在此線上無窮的敬意與感謝。

在兩年半的研究生活當中，感謝學長中祺，及同學博志、俊良、家豪、佑承、廷 揆、文祥、李安、文建、世勛、咸佑...等人以及學弟經勝、志凡、宇志、宬岳、偉廷、

鵬洋、逸咸、思銳、冠儒、富地...等人對我的勉勵與扶持，因為有你們的陪伴使我的 人生增添許多色彩，實驗室裡共同的生活點滴，佔了日常生活的大部分，感謝你們。

最後，感謝我親愛的父母親及我最愛的家人，感謝你們所給予的提攜教誨及支持 鼓勵，還有經濟與精神上的支持，讓我能專心於課業研究中，使我可以順利完成學業，

願將本文獻給關心我的每一個人，共同來分享我的成果與喜悅。

目錄

中文摘要 ... 1

Abstract ... 2

致謝 ... 3

目錄 ... 4

表目錄 ... 6

圖目錄 ... 7

符號說明 ... 10

第一章 緒論 ... 12

1.1 前言 ... 12

1.2 文獻回顧 ... 12

1.3 採用方法 ... 16

1.4 文章安排 ... 16

第二章 物理問題描述 ... 17

2.1 剪切流 ... 17

2.2 自然對流 ... 17

2.3 混合對流 ... 18

2.4 網格選用 ... 18

第三章 數值方法 ... 20

3-1 統御方程式（Navier-Stokes 方程式） ... 20

3-2 時間積分 ... 20

3-3 流通量計算 ... 21

3-4 穩定性分析 ... 21

3-5 人工黏滯 ... 25

3-6 固體壁面 ... 26

3-7 預置矩陣 ... 26

第四章 結果與討論 ... 27

4.1 剪切流場驗證 ... 27

4.2 不同馬赫數於剪切流之影響 ... 28

4.3 自然對流比較 ... 28

4.4 不同溫差於自然對流之影響 ... 29

4.5 傾斜角度之影響 ... 30

4.6 混合對流比較 ... 32

第五章 結論與展望 ... 32

5.1 結論 ... 32

5.2 未來展望 ... 33

參考文獻 ... 35

附錄 A ... 37

附錄 B-1 ... 38

附錄 B-2 ... 40

表目錄

表一、瑞利數 10

³

不同網格分析結果比較。 ... 42

表二、瑞利數 10

⁴

不同網格分析結果比較。 ... 43

表三、瑞利數 10

⁵

不同網格分析結果比較。 ... 44

表四、瑞利數 10

⁶

不同網格分析結果比較。 ... 45

表五、瑞利數 10

³

不同溫差之結果比較。 ... 46

表六、瑞利數 10

⁴

不同溫差之結果比較。 ... 46

表七、瑞利數 10

⁵

不同溫差之結果比較。 ... 47

表八、不同傾斜角度下平均紐塞數變化與文獻[16]比較（A=0.05，ΔT=1

^O

C）。 ... 47

表九、平均紐塞數與文獻結果之比較(Gr=10

⁴

、Re=1000)。 ... 48

表十、平均紐塞數與文獻結果之比較(Gr=10

⁴

、Re=100)。 ... 48

圖目錄

圖 2-1 剪切流之正方形孔穴示意圖。 ... 49

圖 2-2 自然對流之正方形孔穴計算格點與傾斜角度示意圖（81 x 81）。 ... 49

圖 2-3 混合對流之正方形孔穴示意圖。 ... 49

圖 4-1（a）網格 41x41， =0.035

β

。（b）為網格

81x81

， =0.025

β

。 ... 50

圖 4-2 不同雷諾數，水平速度分佈圖，實線(81x81)與虛線(41x41)為實驗結果，圓點 為參考文獻結果。 ... 51

圖 4-3 不同雷諾數，垂直速度分佈圖，實線(81x81)與虛線(41x41)為實驗結果，圓點 為參考文獻結果。 ... 52

圖 4-4 不同雷諾數下，流線分佈圖，(a) Re=100，(b) Re=400 ，(c) Re=1000 ，(d)

Re=3200。 ... 53

圖 4-5 Re=3200，不同馬赫數之流線分佈圖。(a) M=0.001，(b) M=0.5 ，(c) M=0.8。 ... 54

圖 4-6 Re=3200，不同馬赫數下之密度分佈，(a) M=0.5 (b) M=0.8。 ... 55

圖 4-7 Re=3200，不同馬赫數下之馬赫數分佈，(a) M=0.5 (b) M=0.8。 ... 56

圖 4-8 Re=3200，不同馬赫數各面摩擦係數分佈。 (a)左面 (b)右面 (c)下面 (d)上面。 ... 57

圖 4-9 Re=3200，不同馬赫數下之水平速度(x=0.5)與垂直速度(y=0.5)比較圖。... 57

圖 4-10 不同瑞利數下溫度分佈圖，(a)Ra=

10

³(b) Ra=

10

⁴(c) Ra=

10

⁵(d) Ra=

10

⁶。 . 58 圖 4-11 不同瑞利數下水平速度分佈比較圖。 ... 59

圖 4-12 不同瑞利數下垂直速度分佈比較圖。 ... 59

圖 4-13 不同瑞利數下流線圖，(a) Ra=

10

³(b) Ra=

10

⁴(c) Ra=

10

⁵(d) Ra=

10

⁶。 ... 60

圖 4-14 Ra=

10

³不同溫差下溫度分佈圖，(a)溫差 1K (b)10K (c)100K (d)300K。 ... 61

圖 4-15 Ra=

10

³不同溫差下溫度分佈(Y=0.5)。 ... 62

圖 4-16 Ra=

10

³不同溫差下左壁面局部紐塞數分佈。 ... 62

圖 4-17 Ra=10 不同溫差下水平速度分佈圖(X=0.5)。 ... 63 ³ 圖 4-18 Ra=10 不同溫差下垂直速度分佈圖(Y=0.5)。 ... 63 ³ 圖 4-19 Ra=10 不同溫差下流線分佈圖，(a)溫差 1K (b)10K (c)100K(d)300K。 ... 64 ³ 圖 4-20 Ra=10 不同溫差下溫度分佈圖，(a)溫差 1K (b)10K (c)100K(d)300K。 ... 65 ⁴ 圖 4-21 Ra=10 不同溫差下溫度分佈(Y=0.5)。 ... 66 ⁴ 圖 4-22 Ra=10 不同溫差下左壁面局部紐塞數分佈。 ... 66 ⁴ 圖 4-23 Ra=10 不同溫差下水平速度分佈圖(X=0.5)。 ... 67 ⁴ 圖 4-24 Ra=10 不同溫差下垂直速度分佈圖(Y=0.5)。 ... 67 ⁴ 圖 4-25 Ra=10 不同溫差下流線分佈圖，(a)溫差 1K (b)10K (c)100K(d)300K。 ... 68 ⁴ 圖 4-26 Ra=10 不同溫差下溫度分佈圖，(a)溫差 1K (b)10K (c)100K(d)300K。 ... 69 ⁵ 圖 4-27 Ra=10 不同溫差下溫度分佈(Y=0.5)。 ... 70 ⁵ 圖 4-28 Ra=10 不同溫差下左壁面局部紐塞數的分佈。 ... 70 ⁵ 圖 4-29 Ra=10 不同溫差下水平速度分佈圖(X=0.5)。 ... 71 ⁵ 圖 4-30 Ra=10 不同溫差下垂直速度分佈圖(Y=0.5)。 ... 71 ⁵ 圖 4-31 Ra=10 不同溫差下流線分佈圖，(a)溫差 1K (b)10K (c)100K(d)300K。 ... 72 ⁵ 圖 4-32 不同傾斜角度下溫度分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=0)。

... 73 圖 4-33 不同傾斜角度下溫度分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=1)。

... 74 圖 4-34 不同傾斜角度下溫度分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=2)。

... 75 圖 4-35 不同傾斜角度下溫度分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=3)。

... 76 圖 4-36 不同傾斜角度下流線分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=0)。

... 77

圖 4-37 不同傾斜角度下流線分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=1)。

... 78

圖 4-38 不同傾斜角度下流線分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=2)。 ... 79

圖 4-39 不同傾斜角度下流線分佈圖:(a) 30° (b) 60° (c)90° (d)120° (e)150° (N=3)。 ... 80

圖 4-40 傾角30 ，不同波數下局部紐塞數分佈。 ... 81 ^

圖 4-41 傾角60 ，不同波數下局部紐塞數分佈。 ... 81 ^

圖 4-42 傾角90 ，不同波數下局部紐塞數分佈。 ... 82 ^

圖 4-43 傾角120 ，不同波數下局部紐塞數分佈。 ... 82 ^

圖 4-44 傾角150 ，不同波數下局部紐塞數分佈。 ... 83 ^

圖 4-45 Re=1000 不同波數下溫度曲線圖與流線圖。 ... 85

圖 4-46 Re=100 不同波數下溫度曲線圖與流線圖。 ... 86

圖 4-47 (a)雷諾數 1000，不同波數下 NU 變化。(b)雷諾數 100，不同波數下 NU 變化。 ... 88

符號說明

M：馬赫數（Mach number），

RT

top

M U

= γ

^。

g：重力加速度。

Gr：葛拉斯赫夫數（Grashof number），

Gr

=

g

β

H

³

( T

_h −

T

_c

)

ν² 。 H：腔體高度。

L：腔體寬度。

Nu：紐塞爾數（Nusselt number），

x T TK Nu HK

∂

∂

= ∆

0

。

Pr：普朗特數（Prandtl number），Pr=ν α 。

Ra：瑞理數（Rayleigh number），

Ra

=

g

β

H

³

( T

_h −

T

_c

)

αν 。 Re：雷諾數（Reynolds number），Re=

UH ν

。

Ri：理查森數（Richardson number），

Ri

=

Gr

Re²。 t：時間。

T：溫度。

U：上蓋速度。

u：x 軸方向之速度。

U：x 軸方向之相對速度。

v：y 軸方向之速度。

V：y 軸方向之相對速度。

x：水平座標。

y：垂直座標。

希臘符號

α：熱傳導係數。

β：熱膨脹係數。

γ：傾斜角度。

θ：無因次溫度，

( T

−

T

c

) ( T

h −

T

c

)

。 ν：動黏滯係數。

τ：流層間的剪應力。

下標符號

c：低溫壁面。

h：高溫壁面。

第一章 緒論 1.1 前言

傳統不可壓縮流場在熱傳之分析已發展相當完備，應用在空氣上在溫差與壓降低 之情況下也是相當理想，本研究則希望透過傳統的自然對流或混合對流，以可壓縮流 場之統御方程式模擬空氣的熱傳現象。透過這些分析能掌握溫差與壓降等參數之影 響，以協助軟體工程應用能力之提昇。在本研究中以溫差之設定，波型壁之效應來探 討可壓縮流場中，密度受溫度與壓力作用下之影響，浮力效應之增減。至於程式之開 始則先以剪切流場對流速分佈做一驗證，並以不同馬赫數觀察密度梯度之影響，接著 以自然對流結果了解熱傳分析能力，並以波型壁之波數、傾斜角度等參數探討熱傳增 益，最後結合剪切流場與溫差之混合對流做一分析，以了解波數與雷諾數對熱傳之影 響。

1.2 文獻回顧

Ghia 等人 [1] 於 1982 年提供正方形剪切流於雷諾數由 100 至 10000 的流場解 答，利用渦度與流函數處理二維不可壓縮 Navier- Stokes 方程式，其結果顯示當雷諾 數越高則需要越多之網格數，且格點間的疏密分佈要越均勻越好，顯示解析能力取決 於網格數目。該物理模型雖然很單純，但由於剪力效應的影響，計算結果會有一個主 漩渦搭配角落數個漩渦，此文獻提供各雷諾數下，水平與垂直速度斷面、頂面渦度分 佈以及主漩渦、反轉漩渦之座標位置等，方便後續學者比對。

De Vahl Davis [2] 於 1983 年分析左右垂直壁面有溫差而頂部與底部為絕熱之正 方形自然對流問題，使用二階中央差分法處理不可壓縮 Navier-Stokes 方程式，初始 網格數為

21x21

，並漸增加至

41x41

及

81x81

，再利用網格外差的方法去推算下一個 網 格 更 高 精 度 的 結 果 ， 以 此 達 到 誤 差 值 小 於 1% 以 內 的 預 測 值 。 提 供 瑞 利 數

3 6

10 ≤ Ra ≤ 10

之流場結構與熱傳特性。

Adjlout 等人 [3] 於 2002 年在一個具波形壁的二維層流自然對流傾斜方腔體內，

進 行了 左 右壁面熱 差異 的 數值研 究 ，並針 對 不同的傾 斜角度 (0 ~^ 180 ) 、 振 幅^ (A=0.05、0.06、0.075、0.08) 、波數(N=1、3) 和瑞利數(10⁴ ≤

Ra

≤10⁶)，而普朗特數 保持不變。此文章說明不同傾角、振幅和波數明顯影響熱壁面的流場結構與熱傳增 減，其中在傾角影響為0 ，此時只有浮力的效應，流場幾乎趨近於靜止狀態，紐塞數^ 趨近 1。隨著傾角慢慢增至120 平均紐塞數為最高值，因此熱邊界層受到對流之影響^ 而增厚或變薄，此外，熱壁面上的波數與振幅深淺更是影響了局部紐塞數的分佈。

Al-Amiri 和 Khalil Khanafer [4] 於 2007 年在頂部為一流動 U 強制對流與一正弦 曲線的底部表面，進行分析二維混合對流的正弦波表面結構對流場和傳熱特性的數值 模擬研究，針對Gr=10 且垂直壁面為絕熱邊界條件。波數(N=0、1、2、3)改變使得⁴ 底部表面積增加，探討熱傳增益。波形壁振福(A=0、0.025、0.05、0.075)變深，熱壁 內部產升了更多的渦旋。理查森數(Ri=10、1、0.01)的降低(即雷諾數增加)使得強制對 流的流速變快，增加了強制對流的效應。結果表示平均紐塞數顯著的受到這些影響，

並發現在所有參數中以 N=2、Ri=0.01 時，平均紐塞數是最佳的。

Rathish Kumar 和 Shalini [5] 於 2003 年此文獻在描述一個物理模型上下壁面絕 熱，左側壁面為一波形壁的方腔體且左右壁面溫差∆ = 進行了熱傳分析。流場使用

T

1 有限元素法求解，而主要探討的參數:低瑞利數的流場(

Ra ≤ 500

)，波形壁的個數 (N=2、4、8、16)、角度(φ =0 ~ 300^{ ��})與振幅，其中葛拉斯赫夫數定義

*

*

2

K g tK

Gr

β

ν

= ∆ ，

* 12

K =C_f

K

。結果顯示壁面波的個數影響紐塞數很小，主要影響在於振幅大小，發現 波形壁在角度 300^�紐塞數會是最佳的，至於波形壁造型之啟始角度φ =60^�、

a

≤0.1則 時紐塞數會更佳。

Moallemi 和 K. S. Jang[6] 於 1992 年在一個二維混合對流的方腔體進行了數值模

擬，邊界左右側為絕熱，上壁為低溫面且有一流速

U ，下壁為高溫面。頂部的流速

_L

增強了強制對流的效應，此文主要在配合不同的 Gr(葛拉斯赫夫數)、Pr(普朗特數)、

Ri(理查森數)和 Re(雷諾數)等參數，針對雷諾數 500 以下之所有參數互相比對並完整 的整理出數據，由這些數據就可反應出其溫度分佈、流場結構與平均及局部紐塞數的 變化，主要控制參數以雷諾數與葛拉斯赫夫數為主。

Ercan 等人 [7] 於 2007 年進行二維穩態不可壓縮流使用 Navier - Stokes 曲線坐標 方程式求解驅動的傾斜腔流動與非正交網格，即使在夾角很小之造型下，此數值方法 仍能有效且穩定的收斂。驅動傾斜腔流動的網格精度採用 512x512，選定了雷諾數 100 和 1000 及孔穴夾角(15 ~165^{ })數值，流函數(ψ)和渦度(ω)的方程解也得到相當高的精 準。此文章也列出詳細結果，可提供未來學者參考。

Ozoe 等 人 [8] 於 1975 年 測 定 在 長 方 形 通 道 的 層 流 自 然 對 流 介 於 瑞 利 數

3 5

3x10 ~10 熱傳特性，物理模型是從下面加熱且頂部為冷面，另外左右兩側絕緣。在 不同的高寬比( H / W )的橫截面下進行比值=1、2、3、4.2、8.4 和 15.5，以及旋轉傾 斜角度(ψ= 0 ~180^{ })和瑞利數從3x10 ~10 增長，分析各種設定下的熱傳特性以及影^{3 5} 響。此文章表明其瑞利數提高將會提高平均紐塞數，增強熱傳的效應，而且在比較高 寬比和傾斜角之結果影響發現ψ=

⁶⁰

^和

H W

=1 時，平均紐塞數是最佳的。

本論文之數值方法參考 Arnone 等人 [9] 於 1995 年以中央差分法搭配 Jameson 人工黏滯與 Runge Kutta 時間積分處理非穩態流場，透過虛擬時間步階方式，非穩態 之求解可用局部時間步階大小求解，而不受限於邊界層內過小之時間步階所限制，在 圓柱自然週期渦流排放(Karman vortex shedding)預測上相當吻合，以及穿音速氣流通 過葉片之暫態壓力分佈預測上相當一致。

Guoliang He 與 Yinnian He [10] 於 2007 年採用上述之有限體積法與有限元素法 於二維 Navier- Stokes 方程式進行了研究。過程中 FVM 與 FEM 在解 Navier- Stokes 方程具有良好穩定性。此外，對於四邊形(16 個點)和三角形(12 個點)這兩種方案的誤 差估計解進行了介紹，提供了一個數據表證實了有限體積法的效率。

Core 等人 [11] 於 2003 年專門研究非均勻網格或非結構化網格下流通量的計算 方法，主要探討的是在低馬赫數的流動流體在任意接口下如何進行有效率的計算與修 正。此文章使用一個多層次的局部加密的投影方法，進行接口流量的校正，結果顯示 了此方法有效率的、精準的計算非均勻網格或非結構化網格的流量。

Hartnett 與 Minkowycz [12]於 1999 年在方腔體內上下面為絕熱，左側有一流速

U 從下向上移動，此外左右兩面

0 T 、_H

T 可互相對調進行流場分析。主要參數

_C

Ri

( ₂

Re

Ri = Gr

)從 0.01 增加至 100，此文獻結果水平速度會非常劇烈的變化，並且歸納

出三種流場特性:

1

Re

²

>>

Gr

為強制對流(極小的自然對流效應)，

10

1 Re .

0 ≤ Gr

₂

≤

為

混合對流（強制與自然對流互相影響），

1

Re

²

<<

Gr

為自然對流（強制對流可以忽略

不計）。

1.3 採用方法

本研究延續先前有限體積法之數值方法，將二階之速度與溫度微分更改為四階計 算，由先前學長[17]之研究可知擴散項若採高階能免除人工黏滯之需求，能提供較平 滑之數值解答，或許能改善數值之精度，因此本論文初步嘗試將二階之有限體積流場 解答，搭配四階等向性元素之微分法在擴張項之計算，如果數值上能提高其精度與收 斂效率，則可以提供後續程式開發上，採用更高階算則之基礎。由於本研究之題目設 定在自然對流場與低馬赫數之剪切流與混合對流，所以也將有限體積法之積分項，運 用預置矩陣之方式，以縮小特徵值差異過大收斂困難之情形。

1.4 文章安排

本研究論文分為五個章節，架構簡介如下：

第一章：為緒論部份，說明本論文先前有關文獻與採用方法。

第二章：本章在介紹物理問題描述。

第三章：本章介紹統御方程式與數值方法並對時間積分、有限體積法、人工黏滯、邊 界條件的設定並加以說明。

第四章：為模擬結果驗證，將數值模擬答案與參考文獻之結果進行比 對，並探討不同參數變化之影響。

第五章：為結論及未來展望，歸納研究中所得之成果並說明未來發展 之方向。

第二章 物理問題描述

本論文所探討之問題共有三個，各別程式所需輸入之無單位係數定義、邊界條件 設定、初始條件與探討參數之範圍將分述於下，另外分析時所建議之網格數與幾何造 型也將做一說明。

2.1 剪切流

其雷諾數(Reynolds number)定義如下：

Re UH

= ν

（2-1）

其中 U 為孔穴頂部平板向右移動之速度，H 為正方形孔穴之寬度，

ν

為孔穴內空氣 之動黏滯係數。至於頂面速度大小之輸入採用馬赫數值來控制，其定義如下：

RT

top

M U

= γ

（2-2）

其中 T

top

為頂面中央之溫度其值固定在 300K，溫度雖然取 300K 作為初始條件，而四 個壁面以絕熱條件方式處理，故最終其溫度仍會因渦流帶動而有一溫度場。執行之流 場雷諾數有 100、400、1000、3200 四組，最後以雷諾數 3200 進行馬赫數 0.001、0.5、

0.8 等變化以了解可壓縮對流場結構的影響與剪力的變化進行了探討。

2.2 自然對流

瑞利數(Reynolds number)定義如下：

g TL 3

Ra β να

= ∆

（2-3）

其中 g 為重力加速度為熱膨脹係數，T 為溫度，L 為特徵長度，

ν

動黏滯係數，

α

熱 擴散率。一個左右垂直壁面固定溫度下，而上下面為絕熱之二維自然對流，進行瑞利

數10³ ~10⁶之穩態分析。左側為高溫壁面且右側為冷溫壁面，兩垂直壁面固定溫差共 計處理 T∆ =1、10、100 與 300 四個不同溫差，例如溫差為 10K，則左側垂直面將維 持在 305K 而右側垂直面將固定在 295K。另外，選定 Ra=10 探討傾角⁵

30°

、

60°

、

90°

、

120°

、

150°

以及波狀壁面不同波數 N=0~3 於自然對流場結構與熱傳特性的影響。至 於波型壁之描述為一三角函數，可表示為:

 

  −

= y

H A N

x 2 π

cos

1

(2-4)

2.3 混合對流

左右垂直壁面為絕熱，頂部為低溫壁面且有一 U 速度在頂部平板向右移動，底 部為波狀高溫壁面的二維混合對流。其理查森數(Richardson)定義如下:

Re 2

Ri = Gr

（2-5）

其中

Gr

為葛拉斯赫夫數，

Re

為雷諾數。目前 Gr=10 此為固定值，處理 Re=100、⁴ 1000，以及波狀壁面 N=0~3 為主要參數改變，分析其流場結構與熱傳特性。在分析 壁面熱傳導量，皆採用無因次之紐塞數(Nusselt number)表示，其定義如下：

x T TK Nu HK

∂

∂

= ∆

0

（2-6）

其中 K 為局部熱傳導係數，而K 亦為 300K 下之熱傳導係數。 ₀

2.4 網格選用

本文以二維可壓縮 Navier-Stokes 方程式來描述正方形孔穴內空氣之流動現 象，正方形孔穴如圖 2-1 所示。物理模型為高寬比 AR=1，討論各種參數對二維流場 的影響，主要的參數包括瑞利數、雷諾數、理查森數、馬赫數以及溫差影響、傾角影

時間雖短但容易造成解析度不足且無法正確計算出迴流區之淨熱傳量，然而使用網格

數

81x81

的結果與文獻結果相較之下相當的接近，其誤差值也都在合理接受範圍之

內，因此

81x81

（變化率 1.025）為本研究目前採用的網格數做所有數值模擬計算，

此選擇能縮短計算時間且大部分的解析度已能符合文獻所公佈之數據值。

第三章 數值方法 3-1 統御方程式（Navier-Stokes 方程式）

採用 Reynolds-averaged 之二維穩態流場，其守恆型式之統御方程式表示為:

c c v v

E F E F

U

t x y x y

∂ ∂ ∂ ∂

∂ + + = +

∂ ∂ ∂ ∂ ∂

（3-1）

其中

U

、

E _c

、

F _c

、

E _v

及

F _v

如附錄 A 所示。

3-2 時間積分

本論文採用 Arnone [9] 等人所提出之計算概念並搭配 4 階 Runge-Kutta 進行時間積分 運算

( ) ( ) ( ) ( ) ³

1

2 3

1 2

1

2 1 3 1 4 1

U tR U

U

U tR U

U

U tR U

U

U tR U

U

n n

n n

n n

∆ +

=

∆ +

=

∆ +

=

∆ +

=

+

（3-2）

其中，

∆ t

為物理時間。

∑

−

=

∂ + + ∂

∂ + ∂

∂

− ∂

∂

− ∂

= Flux

R A 1

y S F x

E y

F x

E _{C C V V}

（3-3）

其中，

A

為計算網格之面積。

3-3 流通量計算

二維流場其任一截面之流通量有

E

_,j

12

i+ 及

12 , ij+

F

，將採用平均值計算，則如下:

( )

( i, j i , 1 )

1 2 , i

, 1 i j , j i

2 , i 1

2 1 2 1

+ + + +

+

=

+

=

j j

j

F F F

E E E

（3-4）

因此，該截面積上之總流通量為

2 2 , 1 1 i 2 ,

i 1 A F A

E

Flux = + j + j +

_（3-5）

其中，

A ₁

與

A ₂

為該流通量之投影面積。

若以四邊形之網格而言，該計算體積之總流通量為

2 2 , 1 1 i 1 2 , 2 i 2 , i 1 , 1

1 2

i A F A E A F A

E

Flux = + j + + j + j + + j +

∑

^（3-6）

3-4 穩定性分析

在在使用顯式法時，未避免發散採用 von Neumann 法來作分析。考慮統御方程式 如下：

c c v v

E F E F

U

t x y x y

∂ ∂ ∂ ∂

∂ + + = +

∂ ∂ ∂ ∂ ∂

在 x-y 座標下，所採用的計算格點為非均勻網格，為便於分析則將物理 x-y 座標系統 轉換為計算均勻網格點

ξη

座標系統，其結果如下：

c c v v

E F E F

J U

t ξ η ξ η

′ ′ ′ ′

∂ ∂ ∂ ∂

∂ + + = +

∂ ∂ ∂ ∂ ∂

（3-7）

其中

U,F ,G ,F ,G , J _c ′ ′ ′ ′ _{c v v}

_{如附錄 B-1。}

(3-7)式為非線性方程式，將其推導為線性方程式：

 

 





∂

∂ + ∂

∂ + ∂

∂ + ∂

 

 





∂ + ∂

∂

− ∂

∂ =

∂

η ξ η

ξ η

ξ

T Q S Q

S Q A Q

A Q t

Q ²

2 1 2 2 2

2 1 2

1

（3-8）

而

Q

A F Q ,

A ₁ E ^c ₂ ^c

∂

∂ ′

∂ =

∂ ′

=

其中

Q, A ₁ , A ₂ , S ₁ , S ₂ , T ₁

如附錄 B-2。將中央差分法應用於

上式可得：

( ) ( )

 





 





∆

− −

∆

− + ∆

∆ + + −

∆ + + −

 





 





∆ + −

∆

− −

∆ =

−

−

− +

−

− + +

+

− +

− +

− +

+ +

η η

ξ

η ξ

η ξ

2 Q Q

2 Q Q

2 T

Q 2Q S Q

Q 2Q S Q

2 Q A Q

2 Q A Q

t Q Q

n 1 j 1, i n

1 j 1, i n

1 j 1, i n

1 j 1, 1 i

2 n

1 j i, n

j i, n

1 j i, 2

2 n

j 1, i n

j i, n

j 1, i 1

n 1 j i, n

1 j i, 2 n

j 1, - i n

j 1, i 1 n

j i, 1 n

j

J i,

（3-9）

使用 von Neumann 法：令

Q ⁿ _{i, j} = e ^at e ^ik

^x

^ξ e ^ik

^y

^η

代入(3-9)式且等號兩邊各除以Q ，則可

ⁿ _{i, j}

得

( ) ( )

 





 



 −

− − +

+ + −

+ + −

 

 



 − + −

−

− =

−

−

−

− −

− −

∆

η η

ξ

η ξ

2 Δ 2 Δ

2 Δ T

Δη S 2

Δξ S 2

2 Δ 2 Δ A

Δt A

Δ Δ Δ

Δ Δ Δ

1

2 Δ Δ

2 2 Δ Δ

1

Δ Δ

2 Δ Δ

1 t

η η ik

ξ ik ik η η ik

ξ ik ik

η η ik

ξ ik ξ ik

ik

η η ik

ξ ik ξ ik

a ik

y y

x y

y x

y x y

x

y x y

x

e e e

e e e

e e

e e

e e

e e

I J e

其中 I 為單位矩陣，

k _x , k _y

為誤差相角。

再應用三角函數關係

2 e sin e

2 , e cos e

i i

i

i θ θ θ θ

θ

θ = + ⁻ = ^{− −}

_可得

( ) [ ] _{( )} [ ]

) sin

) Δ sin

T

1 ) 2S cos

1 ) 2S cos

) A sin

) A sin

Δt

1

2 y 2 2

1

2 1

t

η η ξ

ξ

η η ξ ξ

η η ξ ξ

∆

∆ ∆

−

−

∆ ∆ +

−

∆ ∆ +

 

 



 ∆

+ ∆

∆ ∆

−

− =

∆

y x

x

y x

a

(k (k

(k (k

(k i (k

I i J e

化簡上式可得：

[ ]

( ) ( )

[ ]

) sin ) sin T

1 cos

2S 1

cos t 2S

) sin A ) sin t A

I

y 1

y 2

1

y 2

1 t

(k (k

k J k

i (k J (k

e

x x

x a

−

− +

∆ − +

∆ +

−

∆ =

（3-10）

上式

e ^{a ∆ t}

為擴大矩陣 t Z) P(

G

G,

J

= ∆ 將(3-10)式對流項與擴散項兩部份探討，以求

得最小時間步階，設定Z= Z_I +Z_R則：

[ ]

( ) ( )

[ ^{2S cos 1 2S cos 1 T sin ) sin )} ]

Z

) sin A ) sin A Z

y 1

y 2

1 R

y 2

1 I

k ( k

( k

k

k ( k

(

x x

x

−

− +

−

=

+

=

為符合穩定性條件

λ

_G ≤1則：

1. 考慮

t Z ) 1 ( _I  ≤



 



 ∆

J

^，亦即

t ( ) Z ) 1

( _{I max}  ≤



 



 ∆ λ

J

^。

^{Z I}

如附錄 B-2 所示，則

Z _I

的四個 特徵值為：

2

2 b

a c U , U ,

U ± +

當

k _x , k _y

_為

π 2

時，可得最大之特徵值為：

( ) Z I max = MAX [ U ± c a ² + b ² ]

λ

因此對流項時間步階

∆ t _c

為：

^{∆ t c = λ 2} ( ) ^{Z 2} _I ^× _max ^J

（3-11）

其中2 2為四階 Runge-Kutta 之容許值。

2.

t Z ) 1 ( _R  ≤



 



 ∆

J

^，則

( ) ( )

[ ]

( ) ( )

[ ^{2S cos 1 2S cos 1 T sin ) sin )} ]

t

) sin ) sin T 1 cos 2S 1 cos t 2S

t Z

y 1

y 2

1

y 1

y 2

1 R

k ( k ( k

J k

k ( k ( k

J k J

x x

x x

−

− +

∆ −

≤

−

− +

∆ −

∆ =

當

k

_x

, k

_y為

π

時，

cos k _x − 1 , cos k _y − 1

有最大值 2，而當

k _x , k _y

為

π

2 _時，

) sin ),

sin ( k _x ( k _y

有最大值 1，因此

[ 1 2 1 ]

R t 4 S 4 S T

t Z ∆ + +

∆ ≤

J J

由於T 的特徵值非常複雜，因此假設格點是近乎直交的，則₁ T 的四個特徵值如下： ₁

( ) ( ) ( ) ( ) ( )

ξ

ρ η

µ λ µ

λ λ

λ

x 3 y

T

0 T T

T

t l 1 4

1 3 1 2 1 1

J

= +

=

=

=

因此

( )

ξ

ρ η

µ

µ

y x

T₁ 3^{l t}

J

= +

S 1

的四個特徵值為

( ) ( ) ( ) ( ) 1 ( ) ( ^{2 2} )

4

1 3 1 2 1 1

x 1 y

S

0 S S

S

η

ρ η

λ

λ λ

λ

− +

=

=

=

=

RJ

γ

K

因此

( ) [ ( ^{2 2} ) ( ^{2 2} ) ]

2

1 4 1 y x y x

S 4 S

4 _{η η ξ ξ}

ρ ^{+ + +}

= −

+ RJ

γ K

所以

( ) ( ) [ ( ) ( ) ]

 



 

 + + − + + +

≤ ∆

∆ l t 2 2 2 2

R 4 1 y x y x

x 3 y

Z t t

ξ ξ η η ξ

η ρ

ρ µ µ

RJ γ K J

J J

則擴散項的時間步階為：

R

v Z

t = J

∆

（3-12）

因此時間步階取兩者中最小之值。

∆ t = MIN ( ∆ t _c , ∆ t _v )

（3-13）

3-5 人工黏滯

在數值計算時，若加上人工黏滯可避免數值震盪。本論文使用 Jameson 的人工黏 滯[9]。考慮一維模擬之方程式

= 0 + _x

t cu

u

_（3-14）

透過人工黏滯之導入

xxxx xx

x

t cu v xu v x u

u + = ₂ ∆ − ₄ ∆ ³

_（3-15）

而(3-14)式之解為

u = e ^{− ikx − i ω t}

，將其帶入(3-15)式，則可得

( )

( ^{v 2 xk 2 v 4 x 3 k 4} )

kct x

i e e

e ^{− ω} = ^{− − ∆ + ∆}

_（3-16）

此時，當

v ₂

與

v ₄

之值皆大於零時，則可使數值震盪現象呈現指數衰退。如果在 方程式中加入

v ₃ ∆ x ² u _xxx

_{，則將產生}

e ^{− i ∆ x}

²

^k

³

^v

³

^t

的項，此項只會修改解答相位的改 變，對於避免數值的震盪並無實質上幫助。由於人工黏滯

u _xx

_{項的加入，產生一階的}

誤差，而

u _xxxx

項則造成三階的誤差，因此

v ₂

不可取的過大，以免影響數值結果的 精確度，本研究只使用四階微分以避免數值震盪的現象，其值皆固定在 0.001。最後，

加上人工黏滯的方程式為

( ) _{x , y}

y , x y ,

x Flux AV

dt

V dU + ∑ =

^（3-17）

3-6 固體壁面

考慮黏滯性流時，由於在固體壁面上之無滑動條件，所以流通量之計算可表示成：

( ) ( )

( )

0

0

wall xx xy

xy wall yy

P dy

E F dS

P dx

τ τ

τ τ

 

 

 − − 

   

+ =  

 − + − 

   

 

 

∫ ∫

∫

（3-18）

其中

P

_wall 為壁面上之靜壓力，

τ _ij

_{為壁面上之剪應力。}

3-7 預置矩陣

Choi 與 Merkle [13] 最早提出使用此矩陣以修改對流項之特徵值以改善收斂狀 況。後續 Weiss 與 Smith [14]之研究。建議

[ ]

[ ]

[ ] ( )

1 0 0

1 0

1 0

1

_p

RT T

u RT u T

v RT v T

H RT u v C H T

ρ

ρ ρ

ρ ρ

ρ ρ ρ

Θ + −

 

 Θ + − 

 

Γ =  Θ + − 

 

 

Θ + −

   

 

（3-19）

其中

Θ = ¹ V

ref²

− ¹ c

² ，

C

為聲速，且參考 Edwards and Liou [15]所提出的參考速度

V

ref ，其定義如下

^{V ref 2 = min } _ ^{c 2 , max} ( ^{V 2 , KV max 2} ) ^ _

（3-20）

其中

V _max

為正方形孔穴內之最大流速，

K

為一常數項於本研究中皆固定在 10。

第四章 結果與討論

本論文探討不同驅動腔上壁馬赫數於剪切流場之影響、不同溫度差下正方形孔穴 之自然對流現象、不同傾斜角度下對正方形孔穴內自然對流之影響，最後再以剪切流 與溫差兩種效應下之混合對流做一分析，透過文獻之比對來驗證程式之開發，並透過 網格數之多寡來確認分析之精度，並利用流線及溫度分佈來了解流場的熱傳機制，提 供熱交換裝置設計開發之參考。

4.1 剪切流場驗證

本研究利用文獻[1]中上方驅動之正方形孔穴內之剪切流場之數值解做為本研究 之程式開發驗證，採用雷諾數 100、400、1000 與 3200 之結果進行比較，由於採用可 壓縮流分析，

ν

之溫度取

300 K

作為初始條件，而四個壁面以絕熱條件方式處理。疊 代後皆將頂面中點之溫度歸零成

300 K

，孔穴正中央之壓力皆歸零為

10 ⁵ Pa

。 本研究以 41×41 與 81x81 兩種網格處理剪切流場之分析，其計算網格如圖 4-1 所 示，其疏密度在四個壁面角落略為密。利用文獻整理之正中央水平速度分佈與計算結 果之比較繪製在圖 4-2，採用 41x41 格點數在壁面附近差異較大，但以格點數 81x81 則非常吻合，顯示以目前四階等向性元素微分之速度與溫度計算仍需以較多之網格數 方能改善其之數值精度，另外半腰垂直速度分佈之比較繪製在圖 4-3，其結果也是需 要格點數 81x81 方能吻合，孔穴內流線之分佈繪製在圖 4-4，與文獻[1]網格數為 129x129 比較，其圖形很相似且在角落部分的迴流區域大小也吻合，由角落處迴流區 之解析能力，可知雷諾數 100 時網格之疏密差異需增加以處理較小之反轉渦流，但對 雷諾數 3200 時，則需增加網格數以改善多層之反正轉漩流的解析。

4.2 不同馬赫數於剪切流之影響

本研究上方驅動之馬赫數大小探討密度變化之可壓縮流場與不可壓縮分析之差 異，採用 0.001、0.5 與 0.8 等三種不同馬赫數進行比較，圖 4-5 為孔穴內流線之分佈，

由此三張圖之比較可以發現流場結構差異很小，只有在右下角與左上角迴區之大小，

因馬赫數較大時動量影響使其區域略為縮減，以及右上角之流線更貼近角落，至於孔 穴中央區域則變化不大，由於頂面馬赫數 0.001 在密度變化之影響只有 10

^-6

之效應所 以幾乎是不可壓縮結果，反觀馬赫數 0.5 與 0.8 應該在密度影響上較明顯，如圖 4-6 所示為其密度分佈，主要影響之範圍只集中在右上角，針對此二馬赫數分別有 19%

與 51%之最大增加量，至於中央區域之密度影響分別只有減少-2%與-8%，因此對整 體流線重新分配之影響有限，則馬赫數分佈如圖 4-7 所示;至於四個牆壁上摩擦係數之 分佈如圖 4-8 所示，有可以透過此圖上摩擦係數為零之點判斷迴流區域之大小隨馬赫 數變化之情形，其結果如同上述之流線觀察結論，至於摩擦係數之變化，左壁面由馬 赫數 0.001 加大至 0.8 其最大摩擦係數增加 11%，至於右、底與頂面則分別增加 8.6%、

27%與 6.9%，圖 4-9 為不同馬赫數正中央水平速度與半腰垂直速度分佈之比較，採用 馬赫數 0.8 其速度邊界層有薄一點，但在外流場之速度大小會減小約 4%。

4.3 自然對流比較

對於熱傳分析之能力則利用文獻[2]之自然對流結果進行驗證，採用之瑞利數有 10

³

、10

⁴

、10

⁵

與 10

⁶

等四個，左右垂直壁面之溫差取 1℃，圖 4-10 為溫度分佈圖與文 獻之結果很類似，不同瑞利數下之正中央水平速度與半腰垂直速度分佈則分別繪於圖 4-11 與 4-12，此一分佈說明隨瑞利數之增加其速度邊界層將貼近壁面，中央區域形 成低速靜止狀態，且溫度梯度方也會漸漸旋轉至垂直方向，亦即對流作用主導之熱傳，

至於流線分佈圖 4-13 在各個瑞利數上也與文獻非常吻合，唯目前 81x81 網格數在瑞 利數 10

⁶

分析上，於中央迴流區之解析稍嫌不足，另外利用文獻之中央最大平速度、

之誤差分析可知，在最大流速、最大與最小紐塞數之預測上誤差只有幾個百分比以 內，只有在最大垂直速度發生之位置很靠近壁面因此誤差會較高，需考慮將目前格點 之疏密度做一調整。

4.4 不同溫差於自然對流之影響

目前採用統御方程式為可壓縮流場來求解，流場中有溫度、壓力、水平方向 速度與垂直方向速度等四個未知數，透過左右壁面溫差的大小孔穴內之流速也隨之增 加，由於密度可隨著溫度與壓力的改變，可以透過此一參數觀察可壓縮流場之影響，

採用之溫差有 1K、10K、100K 與 300K 等四組以及瑞利數 10

³

、10

⁴

與 10

⁵

等三種狀 況，圖 4-14 為瑞里數 10

³

於四種溫差下之溫度分佈圖，在溫差 1K 時溫度分佈幾乎對 稱，但當溫差增加時左側溫度梯度下降而右側溫度梯度增加，透過無單位溫度分佈圖 4-15 可以更明顯觀察出此一現象，顯示左半側之密度較低因此流線寬度略增，相反，

右半側密度略大流線會變窄。圖 4-16 則繪製不同溫差下左側壁面之局部紐塞數分佈，

由此圖可觀察到溫差 1K 與 10K 之結果幾乎相同，但隨著溫差之加大，最大局部紐塞 數會減小而最小局部紐塞數會增加之變化，但是積分此面積之平均紐塞數則幾乎不受 溫差之影響(見表 5)，顯示左側垂直壁面之溫度梯度減小而熱傳導係數增加之結果使 其局部紐塞數趨於均勻，不同溫差下正中央水平速度與則繪於圖 4-17，由此圖可明顯 觀察出上方高溫區之流線向下移，造成向左迴流之低溫流體流線縮減，同樣類似之情 形發生在圖 4-18 半腰垂直速度分佈上，左側高溫區之流線向右移，造成向下流動之 低溫流體流線縮減，並略為增加其垂直方向之流速，圖 4-19 之流線分佈隨不同溫差 變化之結果也描述低溫差下同心圓之結構，隨溫差增加演變至扭曲不同疏密分佈的情 形。

圖 4-20 為瑞利數 10

⁴

於四種溫差下之溫度分佈圖，在溫差 1K 時左下側與右上側 之溫度梯度幾乎對稱，但當溫差增加時左下側溫度梯度下降而右上側溫度梯度增加，

透過無單位溫度分佈圖 4-21 可以更明顯觀察出此一現象，顯示左半側之密度較低因 此流線寬度略增，相反，右半側密度略大流線會變窄。圖 4-22 則繪製不同溫差下左 側壁面之局部紐塞數分佈，由此圖可觀察到溫差 1K 與 10K 之結果幾乎相同，但隨著 溫差之加大，最大局部紐塞數會增加而最小局部紐塞數也會增加，唯中段之局部紐塞

數會減小，而積分此面積之平均紐塞數由溫差 1K 增至 300K 會略為增加 1%(見表 6)，

不同溫差下正中央水平速度與則繪於圖 4-23，由此圖可明顯觀察出上方高溫區之流線 向下移，造成向左迴流之低溫流體流線縮減，同樣類似之情形發生在圖 4-24 半腰垂 直速度分佈上，左側高溫區之流線向右移，造成向下流動之低溫流體流線縮減，並略 為增加其垂直方向之流速，圖 4-25 之流線分佈隨不同溫差變化之結果也描述低溫差 下水平橢圓對稱之結構，隨溫差增加演變至扭曲左右不對稱有疏密分佈的情形。

圖 4-26 為瑞利數 10

⁵

於四種溫差下之溫度分佈圖，在溫差 1K 時左下側與右上側 之溫度梯度幾乎對稱，但當溫差增加時左下側溫度梯度下降而右上側溫度梯度增加，

透過無單位溫度分佈圖 4-27 可以更明顯觀察出此一現象，顯示左半側之密度較低因 此流線寬度略增，相反，右半側密度略大流線會變窄。圖 4-28 則繪製不同溫差下左 側壁面之局部紐塞數分佈，由此圖可觀察到溫差 1K 與 10K 之結果幾乎相同，但隨著 溫差之加大，最大局部紐塞數會增加而最小局部紐塞數也會增加，唯中段之局部紐塞 數會減小，而積分此面積之平均紐塞數由溫差 1K 增至 300K 會略為增加 0.8%(見表 7)，不同溫差下正中央水平速度與則繪於圖 4-29，由此圖可明顯觀察出上方高溫區之 流線向下移，造成向左迴流低溫高密度之流體流線縮減，並明顯減小向左流動之無因 次速度，同樣類似之情形發生在圖 4-30 半腰垂直速度分佈上，左側高溫區之流線向 右移，造成向下流動之低溫流體流線縮減，但隨近牆處密度之增大其垂直方向之無因 次流速會減少，圖 4-31 之流線分佈隨不同溫差變化之結果也描述低溫差下左右兩個 對稱之漩渦結構，隨溫差增加演變至左漩渦增大而右側漩渦扭曲並縮擠至右下角，形 成不對稱有疏密分佈的情形。

表 5、表 6 與表 7 則分別整理出以上圖形中之最大水平速度、最大垂直速度、左 壁最大與最小局部紐塞數以及上述物理量發生之位置，如同前面分析之結果，在高溫 差下因左右壁面密度增減，造成中央渦流向右下側移動，造成高溫側溫度梯度減小，

低溫側溫度梯度增加，但在高溫側熱傳導係數會隨之增大，因此平均紐塞數幾乎維持 不可壓縮流之大小，在瑞利數 10

⁴

與 10

⁵

下僅有 1%之增益，但對局部流場之影響上，

高瑞利數之邊界層更靠近壁面，因此左右溫差加大所產生之密度效應會更明顯，所以 表 7 之結果相較表 5 與表 6 有比較大之增減影響。

4.5 傾斜角度之影響

至於傾斜角度對自然對流之影響，分別採用 30

^o

、60

^o

、90

^o

、120

^o

與 150

^o

五個 角度進行，而加熱壁面之波數則以 N=0、1、2、3 等四種造型分別加以討論，圖 4-32 至圖 4-35 分別為不同波數在五個傾斜角之溫度分佈圖，當傾斜角為 30

^o

至 60

^o

為上熱

下冷，因此浮力作用不強對流較弱，因此溫度邊界層皆較厚，90

^o

以上則對流明顯，

因此溫度邊界層在渦流衝擊處皆較薄，且孔穴中央形成等溫區，但不同波數對中央區 域之影響很小，只有在波形壁面附近有一些出入。圖 4-36 至圖 4-39 則分別為不同波 數在五個傾斜角之流線分佈圖，由這些流線繪製可以發現，當波數增加至 3 個其流線 對中央區域之漩渦解析能力開始變弱，因此越複雜之造型需適度增加格點數以維持其 精度，以瑞利數 10

⁵

而言其中央處有兩個同轉向之漩渦，但隨者傾斜角度增加會越來 越靠近，且會縮小，當傾斜角度在 150

^o

時，此中央兩個漩渦完全消失，形成單一個 迴流。至於波數之影響可以觀察到光滑熱壁之前段邊界層最薄，但隨著邊界層成長而 變厚，若採用一個波則熱壁之前段下陷因此熱傳量會減弱，但主迴流與波峰處之衝擊 會提高熱傳，但後段下陷之波谷又會減少後段之傳熱。若採用兩個波則邊界層需通過 三個谷與兩個峰，會較困難，但主迴流配合到正中央的第二個谷，因此熱傳量因兩個 波峰與中央波谷之協助，會改善熱傳現象。至於採用波數三個之造型，則有四個波谷，

其邊界層在第一個與最後一個波谷明顯增厚，因此需更高之瑞利數才會有優點。

圖 4-40 至圖 4-44 則為不同傾斜角下局部紐塞數分佈圖，圖中比較光滑與不同波 數之差異，與溫度與流線分佈之現象一致，亦即在接近第一波峰前都有一個局部紐塞 數最高點，其衝擊位置越前段其紐塞數峰值越高，因此波數為三之造型有最大之局部 紐塞數，但其波谷之低點也是所有造型最弱之紐塞數，採用光滑之造型其局部紐塞數 之最大值與最小值之差異會最少。至於最大局部紐塞數發生之位置在這五個傾斜角下 可以發現以波數 3 而言不論傾斜角度皆發生在 0.16 處，波數 2 而言不論傾斜角度皆 發生在 0.24 處，但對波數 1 而言則隨傾斜角度之增加由上游 0.1 處向下游移動至 0.37 處，代表高傾斜角度時其漩渦變強，衝擊之位置會更靠近波峰，光滑壁也有類似之情 形，則隨傾斜角度之增加由上游 0.04 處向下游移動至 0.27 處。

表八比對不同波型壁面上平均紐塞數隨傾斜角度之變化，表中之數值與文獻[16]

之結果做一比較，發現誤差皆在 1%以內，顯示不同分析網格與方法其結果仍相當一

致。分析此表可以發現在 30

^o

與 60

^o

傾斜角度，熱壁面隨著增加波型數，傳導之平均 間隔會縮短，因此可以增加平均紐塞數值，但增加非常小的百分比，正方型孔穴流內 傾斜角度為 90

^o

光滑壁可以得到較高的熱傳遞速率；但傾斜角度大於 90

^o

時，兩個波 型結構中所得到的平均紐塞數會稍為高一些，而傾斜角為 120

^o

發生最大平均紐塞數 值可達 4.63。

4.6 混合對流比較

本研究結合前述之剪切流與自然對流之處理能力進行混合對流之探討，採用文獻 [4]所提供之兩筆數據進行；頂面馬赫數M = 0.001，葛拉斯赫夫數10

⁴

但正方形孔穴之 雷諾數Re = 1000與Re=100兩種條件進行分析，所採用的網格數為81×81，其底部加熱 壁面之波數以N = 0，1，2，3等四種造型分別進行討論，雷諾數1000之溫度與流線分 佈如圖4-45所示，而圖4-46則為雷諾數100，比較這些分佈與文獻公佈之結果很接近，

唯在底面之反轉渦流之解析能力仍嫌不足，比較圖4-45與4-46可知高雷諾數之速度與 溫度邊界層皆較薄，且雷諾數1000之孔穴中央幾乎維持均溫之狀態，至於底面局部紐 塞數與文獻[4]之結果同時繪於圖4-47，在這些圖中雖然趨勢接近，但以雷諾數1000 條件差異較小，而雷諾數100採用光滑與波數1之分析差異最大，其最大局部紐塞數之 差異達11%(Pr=1)，當波數增加後，以波數3為例，與文獻之差異則減少至4%，若雷 諾數1000時更可減少至2.7%的差異。最後將文獻[4]與本研究所得之平均紐塞數整理 於表9(Re=1000)與表10(Re=100)中，由此表9可知目前設定之參數其最佳之熱傳發生 於兩個波，雖然波谷對傳熱有所阻礙，但中央之波谷確符合目前剪切流所產生之渦 流，因此透過波峰之增益，整體有最好之熱傳遞，另外在表10中雖然也呈現波數2有 較佳之平均紐塞數，但增加量非常微小，顯示邊界層較厚時，對波數之變化較不敏感。

第五章 結論與展望 5.1 結論

本研究採用二階有限體積法搭配四階等向性元素進行速度與溫度梯度微分，求解 可壓縮 Navier-Stokes 統御方程式，針對剪切流、自然對流與混合對流等問題，探討 不同馬赫數、不同溫差、不同傾斜角度以及波形壁之波數等流場之影響，可得到以下 之結論：

（1）目前分析程式處理正方形孔穴在 Re=3200 以下之剪切流，需採用網格 81x81，

其速度與流線分佈會與文獻公佈之結果較接近，若採用網格 41x41 則壁面邊界 層與角落之迴流場解析會略差。

（2）正方形孔穴頂面馬赫數由 0.001 增加至 0.8 雖然有造成右上角之密度較大之梯度 變化，但整體流場結構沒有明顯之改變，唯局部摩擦係數之變化幅度會微幅加 大(~10%)，且最大無單位流速會略為減少(~5%)。

（3）針對正方形孔穴左右垂直壁面有溫差之自然對流問題，於瑞利數 10

³

~10

⁶

之模 擬與文獻結果皆很一致。在最大流速、最大與最小紐塞數之預測上誤差只有幾 個百分比以內，只有在最大垂直速度發生之位置很靠近壁面因此誤差會較高，

需考慮將目前格點之疏密度做一調整。

（4）在不同溫差之探討上，在高溫差下因左右壁面密度增減，造成中央渦流向右下 側移動，造成高溫側溫度梯度減小，低溫側溫度梯度增加，但在高溫側熱傳導

係數會隨之增大，因此平均紐塞數幾乎維持不可壓縮流之大小，在瑞利數 10

⁴

與 10

⁵

下僅有 1%之增益。

（5）針對不同傾斜角度之自然對流模擬上，平均紐塞數與文獻結果比較最大誤差皆 在 1%以下，其結果建議當傾斜角度在 60

^o

以下時，波數越多越好，可以縮減傳 導之距離，當傾斜角度在 90

^o

時，相較振幅 0.05 以光滑壁最佳，至於傾斜角度 超過 120

^o

時，以波數二之造型最好，因為兩個波峰與中央之波谷皆對熱傳有助 益，整體最大之熱傳量發生在 120

^o

傾斜角度與波數二時。

（6）上方驅動以及上下壁具溫差之混合對流模擬上，與文獻局部最大紐塞數之比較，

光滑加熱面最大誤差有 11%，但波數增加至三個最大誤差只有 4%，主要是差異 量沒有太大變化，但局部最大紐塞數隨波峰增加而造成比例減少。其分析之趨 勢與文獻結果一致，亦即以波數二之造型最好，因為此流場結構對兩個波峰與 中央之波谷皆能有效增加熱傳。

5.2 未來展望

目前分析上遇到一些現象可供未來繼續研究的方向之整理如下：

（1） 探討完全對稱之 0

^o

與 180

^o

傾斜角度下自然對流模擬上無法收斂之原因。

（2） 增加網格數或是更高階之速度與溫度梯度微分以處理高雷諾數下之剪切流、高 瑞利數之自然對流以及低理查森數之混合對流分析。

（3） 探討不同紊流模式之熱傳預測能力。

參考文獻

[1] U. Ghia, K. N. Ghia, C. T. Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a Multigrid Method. Journal of computational physics 48 (1982) 387-411.

[2] G. De Vahl Davis. Natural connection of air in a square cavity a bench mark numerical solution. International Journal for Numerical Methods in Fluids, vol. 3 (1983) 249-264.

[3] L. Adjlout, O. Imine, A. Azzi, M. Belkadi. Laminar natural convection in an inclined cavity with a wavy wall. International Journal of Heat and Mass Transfer 45 (2002) 2141–2152.

[4] Abdalla Al-Amiri, Khalil Khanafer, Joseph Bull, Ioan Pop. Effect of sinusoidal wavy bottom surface on mixed convection heat transfer in a lid-driven cavity. International Journal of Heat and Mass Transfer 50 (2007) 1771–1780.

[5] B.V. Rathish Kumar, Shalini. Free convection in a non-Darcian wavy porous enclosure.

International Journal of Engineering Science 41 (2003) 1827–1848.

[6] M. K. Moallemi and K. S. Jang. Prandtl number effects on laminar mixed convection heat transfer in a lid-driven cavity. HF. J. Hear Mass Transfer. Vol. 35, No. 8 (1992) 1881-1892.

[7] Ercan Erturk and Bahtiyar Dursun. Numerical solutions of 2-D steady incompressible flow in a driven skewed cavity. No. 5 (2007) 377-392.

[8] Hiroyuki Ozoe, Hayatoshi Sayama Stuart and W. Churchill.Natural convection in an inclined rectangular channel at various aspect ratios and angles-experimental measurements. Int. J. Heat Mass Transfer. Vol. 18 (1975) 1425-1431.

[9] Arnone, A., Liou, M. S. and Povinelli, L. A., “Integration of Navier-Stokes equations using dual time stepping and a multigrid method,” AIAA J., Vol. 33, No. 6 (1995)985-990.

[10]Guoliang He and Yinnian He. The finite volume method based on stabilized finite element for the stationary Navier–Stokes problem. Journal of Computational and Applied Mathematics 205 (2007) 651 – 665.

[11]Xavier Core, Philippe Angot, Jean-Claude Latche. A multilevel local mesh refinement projection method for low Mach number flows. Mathematics and Computers in Simulation 61 (2003) 477–488.

[12] J.P. Hartnett and W.J. Minkowycz. Aiding and opposing mechanisms of mixed convection in a shear-and buoyancy-driven cavity. Int. Comm.Heat Mass transfer, Vol.26, Vol. 7 (1999) 1019-1028.

[13] Choi, Y. H. and Merkle, C. L., “The Application of Preconditioning in Viscous Flows,” Journal of Computational Physics, Vol. 105, No. 2, 1993, pp. 207-223.

[14]Weiss, J. M. and Smith, W. A., “Preconditioning Applied to Variable and Constant Density Time-Accurate Flows on Unstructured Meshes,” AIAA Paper 94-2209, June 1994.

[15] Edwards, J. R. and Liou, M. S., “Low-Diffusion Flux-Splitting Methods for Flows at All Speeds,” AIAA Journal, Vol. 36, No. 9, 1998, pp. 1610-1617.

[16] Yang, Y. L., Lin, Y. C., “Numerical Investigation of Natural Convection in an Inclined Wavy Cavity Using a High-order Differential Scheme,” 2010 航太學會學術 研討會，論文編號 C2-2, 2010, 開南大學。

[17] Yang, Y. L., Lin, C. C., “A Mixed Order of Galerkin Discretization for Steady Compressible Navier-Stokes Equations in Lid-Driven Cavity Flows,” 中國機械工程 學會第二十五屆全國學術研討會論文集,2009,論文編號：A06-17，大葉大學。