中 華 大 學

中 華 大 學 碩 士 論 文

微分法於熱傳問題之分析與應用

Analysis and Applications of Heat Transfer Problems Based on the Differential Scheme

系 所 別：機械工程學系碩士班 學號姓名：M10008018 陳冠霖 指導教授：楊 一 龍 博 士

中華民國 102 年 8 月

ii

中文摘要

本研究利用數值方法探討二維自然熱傳之現象，透過文獻結果之比對以驗證程式 之開發，再進行參數之設定以了解流場結構與熱傳增益之變化。其中稜形夾角之改變 對孔穴內流場結構之影響上，以 90^o夾角能提供最強之對流效果，另外瑞利數越大夾 角改變所產生之渦流擠壓現象也越複雜，但當夾角縮小至 15^o或是放大至 165^o時，其 壁面摩擦之影響會使對流減弱成為單一漩渦之結構。將造型改為同心圓時，其瑞利數 增加其平均紐塞數會增加但伴隨流場結構之變化，其圓柱左右迴旋數目與迴旋內部之 渦流數目皆會增多，此一現象可呈現穩態流場發展至不穩定流動之歷程。對偏心圓柱 之問題上採用偏下放置，其熱對流渦流不受阻礙情況下會有熱傳增益。上下週期運動 圓柱上，透過局部疊代數目與物理時間之設定可以控制時間精度之需求。最後針對開 放圓柱之熱傳分析上，主要是探討上下游範圍對圓柱表面紐塞數之影響。

關鍵字：瑞利數、紐塞數、稜形夾角、同心圓柱、偏心圓柱、動態圓柱、開放圓柱

iii

Abstract

Two-dimensional natural heat transfer problems were studied in the work. Numerical validation with the previous results was performed. The flow structure and heat transfer enhancement were analyzed through a series of parametric constant. In the rhombic angle study of the cavity, the slant angle at 90^o gives the highest heat transfer rate. The higher Rayleigh number, the more complex is the flow structure within the cavity. At slant angle at 15^o or 165^o, the smearing effect of the walls slows down the vortex motion within the cavity, and the structure of the flow becomes a single vortex. For the concentric cylinders problem, the higher th Rayleigh number the larger the mean Nusselt number. At the same time, the flow structure in the annulus diverse to more and more rolls and also more vortices within each roll. The phenomena explain the possible trail of the steady laminar flow to unsteady turbulent flow. By putting the hot cylinder lower gives the buoyant jet more space to travel and transfers more heat to the exterior for the eccentric cylinders problem. Dual time stepping was used to solve the unsteady problem, the local iteration numbers and physical time step can be used to control the temporal resolution of the results.

For the heat transfer of cylinder in the open space, the distance of the far field is examined on the accurate of the local Nusselt number.

Keyword: Rayleigh number, Nusselt number, rhombic angle, concentric cylinders, eccentric cylinders, moving cylinder, flow over a cylinder

iv

誌謝

首先誠摯的感謝指導教授 楊一龍博士，老師的悉心指導使得我得以在碩班時期 間了解人生理念與實務領域的深奧，使我在碩士兩年期間受益匪淺，老師對做人處事 與學問上的嚴謹更是我學習的典範。

本論文承蒙傅武雄博士、鄭藏勝博士，對本文不厭其煩的細心指出我研究中的缺 失，且為我解惑，因為有個委員的建議與指正，使得本論文更加完善，學生在此深表 謝意。

在我求學的日子裡受到機械系上許多老師的照顧，在此特別感謝，陳春福老師、

蔡永培老師、許隆結老師、陳俊宏老師、黃國饒老師、葉明勳老師、邱奕契老師、牛 仰堯老師、鄭藏勝老師與在這段時間一直給我諄諄教誨的蔡博章老師，因為有你們使 的學生在學業上與做人處事上都以各位老師為表率，沒有各位老師也就沒有今天的我，

在此特別感謝。

在這碩班兩年期間特別感謝我的實驗室同學祐達，學弟柏昶、易達，室友芳瑜，

及碩班的福恩、佑政、穩哲、國哲、威昇、小狗、孝杰、證勳、正君、仕達、健聖、

阿凱，還有遠在台科的翔洲、元智的志安，以及碩士班的學長宏澧、冠儒、經勝、富 弟、志凡、大帥在研究的過程中的協助與勉勵，並且陪我走過在實驗室裡的共同生活 點滴。

最後在此感謝我的父親陳忠勇先生、母親黃麗樺女士、兄陳冠穎以及妹陳貞如，

有你們的默默支持與鼓勵，才能讓我在無後顧之憂的完成學業，謹以此文獻給我最摯 愛的家人。

陳冠霖 2013/8

v

目錄

目錄

中文摘要 ... ii

Abstract ... iii

誌謝 ... iv

目錄 ...v

表目錄 ... vii

圖目錄 ... viii

符號定義 ... ix

第一章 緒論... 1

1-1 前言 ...1

1-2 文獻回顧...1

1-3 採用方法...5

1-4 文章安排...6

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

2.1 不同傾斜角度下自然對流流場 ...7

2.2 靜止與動態圓柱之自然對流流場 ...8

2.3 開放空間下圓柱強制對流流場 ...9

2.4 紐塞數 ... 10

2.4 基本假設 ... 11

第三章 數值方法 ... 12

3.1 Navier-Stokes 方程式... 12

3.2 時間積分 ... 13

vi

3.3 流通量計算 ... 14

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

4.1 非正交網格下自然對流分析驗證 ... 16

4.2 菱形孔穴傾斜角度對自然對流之影響 ... 17

4.3 同心圓柱之自然對流流場... 19

4.4 動態圓柱之自然對流流場... 20

4.5 不同上下游邊界設定對流場與熱傳之影響 ... 22

第五章 結論 ... 23

參考文獻 ... 24

vii

表目錄

表一、壁面傾角 90^O下之平均紐塞爾數 ...28

表二、 不同壁面傾角下之平均紐塞爾數 ...28

表三、 不同瑞利數下內外圓柱之平均紐塞爾數 ...28

表四、 在 = 20、不同上下游距離之阻力係數與文獻比較圖 ...29

viii

圖目錄

圖 2.1 傾斜孔穴物理模型示意圖 ...7

圖 2.2 圓柱物理模型示意圖 ...8

圖 2.3 開方空間下流體流過圓柱之物理模型示意圖 ...9

圖 4-1 傾斜角度30 、 = 10 、 = 1000與文獻[16]流線溫度之比較圖 ...30

圖 4-2 在傾斜角度30 、 = 10 、 = 1000與文獻[16]局部紐賽爾數比較圖 ...30

圖 4-3 在 = 10 下流線分布圖 ...31

圖 4-4 在 = 10 下流線分布圖 ...32

圖 4-5 在 = 10 下流線分布圖 ...33

圖 4-6 在 = 10 下流線分布圖 ...34

圖 4-7 在 = 10 下溫度分布圖 ...35

圖 4-8 在 = 10 之溫度分布圖 ...36

圖 4-9 在 = 10 之溫度分布圖 ...37

圖 4-10 在 = 10 之溫度分布圖 ...38

圖 4-11 不同傾斜角度與瑞利數下平均紐塞數 ...39

圖 4-12 在 = 5 × 10 同心圓柱 120^O 處溫度與文獻[4]比較圖 ...40

圖 4-13 在 = 5 × 10 局部紐塞數與文獻[4]比較圖 ...40

圖 4-14 在 = 1 × 10 下流線(左)與溫度(右)分布圖 ...41

圖 4-15 在 = 5 × 10 下流線(左)與溫度(右)分布圖 ...41

圖 4-16 在 = 1 × 10 下流線(左)與溫度(右)分布圖 ...42

圖 4-17 在 = 5 × 10 下流線(左)與溫度(右)分布圖 ...42

圖 4-18 在 = 1 × 10 下流線(左)與溫度(右)分布圖 ...43

圖 4-19 在不同瑞利數下外圓局部紐塞數分布圖...44

圖 4-20 瑞利數 = 5 × 10 、 = 1、 = 1，內外圓柱平均紐塞數發展歷程 ...45

圖 4-21 = 0 時之溫度(上)、流線(下)分布與文獻(左)之比較圖 ...46

圖 4-22 = 14 時之溫度(上)、流線(下)分布與文獻(左)之比較圖 ...47

圖 4-23 = 12 時之溫度(上)、流線(下)分布與文獻(左)之比較圖 ...48

圖 4-24 = 34 時之溫度(上)、流線(下)分布與文獻(左)之比較圖 ...49

圖 4-25 不同週期時局部紐賽數分布圖 ...50

圖 4-26 在 = 20、不同上下游大小下，溫度與VORTICITY和文獻[?](A)比較圖 ...51

圖 4-27 在 = 20、不同上下游距離設定下圓柱流線圖 ...52

圖 4-28 在 = 20、不同上下游距離設定下溫度分布 ...53 圖 4-29 在 = 20、不同上下游大小下，局部紐塞數分布與文獻[29]、[30]比較圖 .53

ix

符號定義

局部音速

, 流速流通量

, 壓力流通量

內插函數 總焓

局部熱傳導係數

平均溫度下熱傳導係數

, 特徵長度

馬赫數, M = ⁄

平均紐塞數, ( ) = ∫

局部紐塞數, Nu =

∆

耦合之點數 壓力

普朗特數, Pr = ⁄ 密度

瑞利數(Rayleigh number),Ra = g∆TL T ν α⁄

雷諾數(Reynold number),Re = ⁄

x

溫度 平均溫度 高溫溫度 低溫溫度

高溫與低溫溫差 x 方向之速度 y 方向之速度

平均溫度下動黏滯係數 平均溫度下熱擴散係數 黏滯係數

法線方向 流體剪應力 預置矩陣 時間步階 物理時間步階 夾角角度

= ⁄ ⁄

1

第一章 緒論

1-1 前言

由於近年來數值模擬之蓬勃發展，對於計算能力的方面越來越要求，對於計算元 件之幾何造型越來越複雜，且隨著工程上的要求對於模擬的精度也越來越要求，但傳 統的數值計算方法大多以二階精度來處理且多為不可壓縮流，對於元件分析上可靠度 不足且無法觀察到溫差與流動馬赫數之影響，本研究則希望透過高扭曲的幾何造型，

以可壓縮流場之統御方程式模擬元件之熱傳現象，並且從高度扭曲的幾何造型下探討 程式之解析能力，接著再以圓形模型探討其流場與熱傳特性，並探討振動圓柱下之動 態流場與熱傳特性，建構一動態分析程式以供工程應用，最後再以開放空間下流過一 圓柱流場來進行不同的 Fair field 大小對於熱傳分析結果所造成之影響作探討。

1-2 文獻回顧

本文採用微分法於熱傳問題之分析與應用，首先進行菱形孔穴內自然對流流場與 熱傳現象之探討，再來以同心圓柱與偏心圓柱之自然對流進行程式的驗證，並進行不 同參數之比較，進而再模擬圓柱內圓上下移動之問題，最後進行開放圓柱的上下游邊 界來進行探討，本節將與本研究有關之文獻依章節之先後順序整理如下：

1964 年 Noh[1]提出了一種結合兩種座標系統特色的方法，是以有限差分法來處 理後續學者 Hirt 等人[2]於 1974 年將此稱為 ALE(Arbitrary Lagrangian -Eulerian)法，

並且利用有限差分法探討其穩定性與精度等問題。但由於有限差分法在不規則曲面時 較不易處理，所以 Thomas 等人[3]於 1981 年就改善了有限差分處理不規則曲面不易 的問題，改用容易處理不規則曲面的有限元素法來處理。

2

1976 年 Kuehn and Goldstein[4]進行了實驗與模擬同心圓柱自然對流流場，實驗 部分固定內外徑半徑大小比例為0.8，不同瑞利數範圍(Ra = 2.11 × 10 ~9.76 × 10 )，

並利用有限差分法進行不同瑞利數(Ra = 10 ~10 )與普朗特數(Pr = 0.7~100)之模 擬，提供了不同角度之溫度分布，與局部紐塞數分布圖之數據，供後人研究。

1983 年 De Vahl Davis [5]分析左右垂直壁面有溫差而頂部與底部為絕熱之正方形 自然對流問題，使用二階中央差分法處理不可壓縮 Navier-Stokes 方程式，初始網格 數為 21x21，並漸增加至 41x41 及 81x81，再利用網格外差的方法去推算下一個網格 更高精度的結果，以此達到誤差值小於 1%以內的預測值。提供瑞利數Ra = 10 ~10 之流場結構與熱傳特性。

1997 年 Cheng and Hong[6]利用非慣性座標系(Non-inertial reference fram)進行了 開放空間下圓柱上下振動之熱傳研究，其研究不同範圍之雷諾數(Re = 0~300)、振 幅(0~0.7)、振動頻率(0~0.3)對熱傳所造成之影響，結果顯示雷諾數(Re)越大熱傳越 好，且隨振幅增加而增加，當頻率與圓柱之自然頻率相同時其熱傳此時為最佳，但是 採用不可壓縮流，所以並未考慮摩擦生熱所帶來之效應，此一效應或許在高頻振動下 以可壓縮流場分析能產生一些差異。

1999 年 Yoo[7]分析不同瑞利數(Ra = 10 ~8 × 10 )與不同普朗特數(Pr = 0.3~1) 下同心圓柱多重解發生之範圍，並找出了熱不穩定的臨界瑞利數，其研究指出多重解 發生可由於普朗特數(Pr)和初始條件所造成，且在普郎特數 Pr=0.3~0.5 區間內可由零 的初始條件而得到多重解，但在 Pr=0.6~1 之間則不行，且隨普朗特數(Pr)的增加，使 得臨界瑞利數(Ra)上升讓多重解存在之機會降低。

2000 年 Cheng and Hung [8]利用有限體積法加上 two-stage pressure correction method 求解二維可壓縮 Navier-Stokes 方程式，模擬了水平往復汽缸運動，並且提供 了在不同轉速(RPM=300,600,1800)與有無考慮重力下週期運動的暫態過程，其結果發

3

現當轉速越高時流場要達到週期運動所需之週期數增加，而當沒有考慮重力項時流場 會呈現對稱的情形，且對稱於水平中心軸。

2001 年 Fu and Tong[9]於利用有限元素法求解不可壓縮流 Navier-Stokes 方程式，

探 討 一 圓 柱 於 空 氣 流 經 水 平 渠 道 內 上 下 週 期 運 動 之 熱 傳 問 題 ， 其 雷 諾 數 為 Re = 100,200,500，頻率Ω = 0.1,0.2,0.4，與振幅0.2,0.4,0.8，其結果顯示了當 Vortex shedding 的頻率與擺動頻率相近時，流場進入 Lock-in regime 速度場與溫度場會成週 期狀態，且紐塞數會明顯增加，然而紐塞數也會隨著雷諾數與擺動速度的增加而增 加。

2008 年 Alawadhi[10]利用 ALE 法求解一圓形空穴內圓柱上下振動之自然對流熱 傳分析，其研究不同瑞利數(Ra = 1 × 10 ~1 × 10 )、振幅(ε = 0.25~0.75)與頻率 (Ω = 0.5~2)，結果顯示當瑞利數(Ra)降低時，內圓紐塞數變化也會減弱，但是對外圓 則是相反，反而會有很大的影響；當提高振幅時，內圓圓柱之平均紐塞數波動會明顯 提高，但對外圓確是有減弱的趨勢；最後當頻率提高時對內圓的平均紐塞數影響增加，

但對外圓來說會隨著頻率提高而影響減少。

2008 年 Barozzi 等人[11]以商業軟體 Fluent 對一方形空穴內部一高溫圓柱之模型 下，進行不同特徵長度定義下的長寬比(d = 0.2,0.4,0.6,0.8)與瑞利數(Ra = 10 ~5 × 10 )下之流場特性探討，其流場主要可以分成穩態、穩態對稱、穩態非對稱、週期性 非穩態、非週期非穩態，並建立圖表供後人驗證與應用。

2010 年 Yang 與 Lin [12]採用微分法處理不同傾斜角度、不同波數以及不同深度 之波型壁之自然對流影響，其分析之結果與實驗結果吻合，亦即在120 傾斜角度下熱 傳最佳，至於波型壁之波數在傾斜角度小時越多越好，其熱傳以傳導為主，但在90 傾 斜角度下時則以光滑壁最佳，但波深度大於 0.3 後，波型壁仍會優於光滑壁，至於傾 斜角度超過120 採用 2 個正弦波之造型其對流效果最佳，其左右兩個角落之迴流效應

4

很弱，浮力與對流作用能加成傳遞至上壁，雖然振輻越大越好，但熱傳增益僅有幾個 百分比，反觀180 傾斜角度其浮力效應與渦流對衝，造成孔穴內有三個角落之明顯 反轉漩渦，形成局部絕熱之現象，因此雖然浮力效應最大但無法配合對流則熱傳遞仍 是較少。

2010 年 Saleha 等人[13]進行了添加奈米粒子在不同壁面傾斜角度(0 ~23.96 )與 不同 Grashofnumber (10 ~10 )及不同體積分率(ϕ = 0~0.05)下的左熱右冷梯形自然 對流流場分析，其結果顯示體積分率越高與傾斜角與 Grashof number 越高熱傳就越高，

並且也提供了經驗公式供工程上使用。

2010 年 Fu 等人[14] 採用 Roe 解法透過預置矩陣與雙時間步階積分，分析三維管 道內以一週期移動塊對熱傳增益之影響，由於移動塊對壁面之擠壓以及捲繞出之渦流 等效應，皆會造成下游壁面紐塞數之增加，分析之實例最大之熱傳增益可達 50%。

2011 年 Liao[15]採用微分法之方法探討具波形壁之二維渠道於可壓縮流場下強 制熱傳分析，其入口邊界採用完全發展流之設定，出口則以流場內部之解答外插出口 處之u、v 與 s(熵)，至於四個面邊界之壓力皆利用法線方向之動量方程式外插求出，

最後再平移整體之壓力分佈以符合入口中央處之壓力設定值(P = 10 Pa)，本文之開 放問題則是採用此邊界條件之設定再加以改良。

2012 年 Anandalakshmi and Basak[16]採用有限元素法，進行左右冷壁面，下壁為 熱 壁 ( 均 溫 ， 非 均 溫 ) 之 菱 形 空 穴 ， 在 不 同 瑞 利 數 (Ra = 10 ~10 ) 、 普 朗 特 數 (Pr = 0.015~1000)、壁面傾斜角度(30 , 45 , 75 ) 下自然對流流場分析，由於傾斜角 與壁溫設定所造成之溫度不連續導致在四個角落之平均紐塞數均為無窮大。當 Ra = 10 時流場主要是以傳導為主，在不同普朗特數、不同壁面傾斜角下其流線與溫 度分布圖都維持一個主迴旋；在高瑞利數(Ra = 10 )、普朗特數(Pr = 0.015、0.7)、

所有壁面傾斜角下都會產生多個渦旋，但在Pr = 7.2、1000和壁面傾斜角度為75 時

5

會產生一兩個渦漩對稱的情況，且在所有普朗特數與不同壁面傾斜角下，其壁面平均 紐塞數最大值皆出現在左面，主要原因是由於左側冷壁靠近底面熱壁所造成，也因為 要達到熱平衡的情況下最小壁面平均紐塞數皆出現於右側冷壁面，且在傾斜角度30 時為最低，且總體熱傳量則是在高瑞數(Ra = 10 )與高普朗特數(Pr = 1000)以及傾 斜角度大(75 )時為最佳。

2012 年 Lin and Liao[17]利用沉浸邊界法求解不可壓縮流，橢圓形圓柱在方形孔 穴內之熱傳特性之研究，其分別進行了不同瑞利數下(Ra = 10 ~10 )、不同長短軸 ( ⁄ = 0.1~1；b 為短軸，a 為長軸)比例、不同內部橢圓傾斜角度(0 ~90 )、與不同 的里查森數(Ri = ⁄ , = 2 ⁄ )與不同迴轉速度(Rpm)對熱傳增益之影ν 響。在相同瑞利數，不同長短軸比例下，平均紐塞數會隨長短軸比例增加而增加，而 局部紐賽數最大值與最小值的差異會隨長短軸比例越小而越大，且會隨著瑞利數增加 而趨近明顯；在相同瑞利數與長短軸比例下，不同傾斜角度下，平均紐塞數以傾斜角 度為90 (即垂直擺放)的情況下為最佳，傾斜角度為0 (即水平擺放)最差，且紐塞數會 隨長短軸比例增加而減少和隨瑞利數增加而增加；最後在內部橢圓形圓柱旋轉的情況 下以 Ri=1.7、瑞利數為10 、長短軸比例為 0.5、Ω < 240RPM時為最佳，可達 16%的 熱傳增益。

1-3 採用方法

本研究延續先前等向性元素微分法之數值方法[12]，在速度與溫度微分採用九個 點的微分，由先前研究可知採用九個點有較好的數值穩定度，且能提供較精確之數值 解答，因此本研究在擴散項上將採用此方法來計算，對流項方面採用三、五、七、九 方式之中央配置方式，當在中央處維持九個點，在壁面則以三個點來進行計算，壓力 計算皆以法線方向之動量方程式外插求出，在封閉情況下以疊代滿足內部質量為一定 值之方式來調整絕對壓力，使得內部質量滿足初始質量達到質量守恆，開放方面維持

6

先前學長[15]之於開放上游入口處之壓力為一大氣壓，至於上游之流速設定上採用 Potential flow 之假設，而出口邊界則是以流線方向進行外插。另外，本研究之結果需 與文獻不可壓縮流數據比對，因此採用低馬赫數設定，因此需運用預置矩陣之方式，

以縮小可壓流場之特徵值差異過大以避免收斂困難之發生，另外也加入時間項，並與 預置矩陣整合以符合暫態流場時間精度之需求。

1-4 文章安排

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

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

第二章：介紹本研究中所探討之物理現象與無單位係數定義。

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

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

第五章：為結論歸納出研究中所得之成果。

7

第二章 物理問題描述

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

2.1 不同傾斜角度下自然對流流場

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

Ra = ^∆ (2.1)

其中 g 為重力加速度，β為熱膨脹係數，T 為平均溫度，H為特徵長度，ν 為平 均溫度下動黏滯係數，α 為平均溫度下熱擴散係數。模型為一左右垂直壁面固定溫度 下，上下壁面為絕熱之方形孔穴，探討不同傾斜角度(15 ~165 )，瑞利數(10 ~10 )之 穩態流場結構與熱傳特性，左側為高溫壁面且右側為冷溫壁面，兩側固定溫差為 ΔT = 1 ，相關物理性質均以300K來計算以配合傳統不可壓縮流流場，模型示意圖 如下：

圖 2.1 傾斜孔穴物理模型示意圖

8

2.2 靜止與動態圓柱之自然對流流場

為內圓半徑為 0.625，外圓半徑為 1.625，外圓為高溫壁面，內圓為冷溫壁面，兩側固 定溫差為ΔT = 1 ，ε為偏心距離，探討不同瑞利數(10 ~10 )和動態下之流場結構 與熱傳特性，相關物理性質均以300K來計算，其物理模型如下：

圖 2.2 圓柱物理模型示意圖 至於動態偏心的距離之描述為一三角函數，可表示為：

ε = ε sin(2 Ωt) (2.2) 其中ε 為最大偏心距離，即振幅。t 為物理時間，為了配合傳統不可壓縮流之 解答這裡的Ω為無因次化頻率，其定義如下：

= ⁄ (2.3)

其中α 為平均溫度下熱擴散係數，ω為有單位之頻率、L 為(R − R )。

T_H



T_L

g X

Y

=_maxSin(2t)

9

2.3 開放空間下圓柱強制對流流場

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

Re = (2.4)

其中U 為入口速度，D為圓柱直徑大小， 為平均溫度下動黏滯係數。至於入口 速度大小是採用入口馬赫數值來控制，其定義如下：

M = (2.5) 其中T 為平均溫度。邊界條件為外界高溫，圓柱表面為冷溫壁面，入口端採用 Potential flow之設定其定義如下：

U = U cos 1 − V = −U sin 1 − (2.6) 其中R為徑向方向之距離，ϕ為角度。本節選定Re = 20，在不同Far field大小 (5D, 10D, 15D, 20D)下探討對流場所造成的影響，以供未來進行分析時參考。，其物 理模型如下：

圖 2.3 開方空間下流體流過圓柱之物理模型示意圖

D u=U

v=V T=T_L N-mom.

(u)_S=0 (v)_S=0 (Entropy)_SS=0 N-mom.

u=0,v=0,T_H Stationary cylinder

n-D

X Y

10

至於圓柱表面所受到之阻力，皆採用無因次之阻係數(Drag Coefficient)表示，其 定義如下：

C = (2.7)

2.4 紐塞數

在分析壁面熱傳導量之多寡上，皆採用無因次之紐塞數( Nusselt number)表示，

其定義如下：

Nu =∆ (2.8)

其中

n

Nu(x) = ∫ (2.9)

本論文中為了與文獻進行比較在此定義局部等效傳導率k (Local Equivalent Conductivity)，其定義如下：

k = R ln( ⁄ ) Nu at R = R

k = R ln( ⁄ ) Nu at R = R (2.10) 其中R 與R 分別代表圓柱之內圓半徑與外圓半徑。

11

2.4 基本假設

本研究採用二維 Navier-Stokes 方程式描述自然對流對流流場之分析，並做了以 下之假設：

1. 流場為層流。

2. 假設空氣為理想氣體。

3. 忽略重力效應(開放空間問題)

4. 黏滯係數(μ)、熱傳導係數(k)以平均溫度T = 300 計算。

5. 比熱的比(γ)為一常數。

12

第三章 數值方法

3.1 Navier-Stokes 方程式

本研究採用 ALE (Arbitrary Lagrangian-Eulerian)的方法並將其加入二維非穩態 Navier-Stokes 方程式，其守恆方程式可表示如下:

+ + + + = (3-1)

其中U為守恆變數[ , u, v, ] ， , 為壓力流通量， , 為流速流通量，S為源 項，其符號定義為

E = [0, , 0,0]

F = [0,0, , 0]

E = ρU, ρuU − τ , ρvU − τ , ρUH − uτ − vτ − ∂T⁄

F = ρV, ρuV − τ , ρvV − τ , ρVH − uτ − vτ − ∂T⁄

S = [0,0, −ρg, −ρgv − ẋ ⁄ − ̇ ⁄ ]

其中ρ為密度，u為水平速度，v為垂直速度，ẋ為網格水平速度，ẏ為網格垂直速 度，U = u − ẋ為相對水平速度，V = v − ẏ為相對垂直速度，P為壓力，T為溫度，E為 總內能， H為總焓，τ 流體之剪應力。其中流速流通量與黏性流通量無論高階與低階 微分皆採用中央差分進行計算，至於壓力流通量在近壁面採用外插方式計算，最後在 疊代以滿足初始質量來滿足質量守恆。

13

為了改善在低馬赫數情況下收斂困難，因此將方程式(3-1)轉換成以下形式：

Γ + + + + + = (3-2)

其中U 為原始變數[ , u, v, ] ，Γ為預置矩陣。Choi 與 Merkle [18] 最早提出使 用預置矩陣以修改對流項之特徵值以改善收斂狀況。後續 Weiss 與 Smith [19]之研 究。建議

Γ =

⎣⎢

⎢⎢

⎡ θ + 1 RT⁄ 0 0 − ρ T⁄ u[θ + 1 RT⁄ ] ρ 0 − ρu T⁄ v[θ + 1 RT⁄ ] 0 ρ − ρv T⁄ H[θ + 1 RT⁄ ] − 1 ρu ρv ρ C − (H T⁄ ) ⎦⎥⎥⎥⎤

(3-3)

其中θ = 1⁄ − 1⁄ ， C為音速。

並參考 Edwards and Liou [11]所提出的參考速度V ，其定義如下

V = [ , ( )] (3-4) 其中V 為流道中最大流速而K為一常數，並設為 1。

3.2 時間積分

本研究中在物理時間項採用二階 backward 差分，並配合預置矩陣，進行雙時間 步階之積分，方程式(3-2)可以展開為下列式

Γ + + + + + = (3-5)

其中 可與預置矩陣項合併成為新的預置矩陣Γ = Γ + ，最後虛擬時間項則

以三階 Runge-Kutta 方式進行整體時間積分以獲得非穩態流場解答

14

U = + Δ (Γ ) ( ) U =3

4 +1

4 +1

4Δ (Γ ) ( )

U = + + Δ (Γ ) ( ) (3-6)

其中  為虛擬時間，而

R = − − + − + + (3-7)

3.3 流通量計算

Irons [20] 首先以等向性元素來進行任意網格之微分與積分。此任意四邊形之元 素，由於方便描述任意曲線之邊緣或是網格的疏密控制，因此廣泛應用於固體力學之 結構分析上，其等向性元素座標之處理在於物理座標( , )，與數值計算之自然座標 (ξ, η)以多項式方式進行耦合(mapping)，以二維幾何可表示為

x( , ) = ∑ ( , ) ( , ) = ∑ ( , ) (3-8)

其中H 為內插函數(Interpolation Function)，N為耦合之點數，本研究在壁面採取N = 3，

而中央區域最多加至N = 9處理。至於流場變數與通量計算也採用相同之方式處理，

以ϕ( , )做一代表，其流場之物理量分佈與計算網格ξ與η之耦合，亦可表示為

ϕ( , ) = ∑ ( , ) (3-9)

而此一物理量之計算網格二維微分可透過下式計算之

15

= +

= + = [ℐ] (3-10)

其中[ ]為 Jacobin 矩陣，透過反矩陣可將(3-8)之結果轉換成物理座標下之微分值，

亦即

= [ ] (3-11) 而方程式(3-9)中[ ] 矩陣數值會隨幾何造型之大小、方向與採用之耦合之點數而變 化。

16

第四章 結果與討論

本論文所採用之數值方法可透過文獻之結果進行驗證，再以此程式探討複雜流 場之熱傳現象，首先以稜形孔穴之自然對流熱傳問題，進行網格扭曲對數值分析結 果之影響，由文獻之流場與熱傳結果可驗證非正交網格下之解析能力，接著再以此 程式探討稜形孔穴固定垂直壁面溫差下，隨稜形夾角下之自然對流問題。接著再以 同心圓形之造型探討其流場與熱傳特性，並以上下運動之圓柱進行動態流場與熱傳 之分析，建構一動態分析程式以供工程應用，最後以開放空間下流過一熱圓柱之流 場來進行不同的上下游計算距離對於熱傳所造成之結果進行分析。

4.1 非正交網格下自然對流分析驗證

首先以方型孔穴內之自然對流在正交網格下的平均紐賽爾數來驗證程式，如表 一所列，可以看出與文獻的計算結果相當吻合，接著再以 Anandalakshmi and Basak [16]所提供之自然對流流場解答，其物理模型為上壁絕熱，左右壁面為冷面，下壁 為熱壁面，溫差為 1^oC，來進行非正交網格下的自然對流分析驗證，圖 4-1 為瑞理 數 10⁵、傾斜角度 30^o、Pr = 1000之流線溫度分布圖，圖中可以看出流線與溫度分 布圖上與文獻相當吻合，且由於壁面較靠近底面熱壁，浮力導致溫度梯度較集中於 左邊冷壁面，相對的右邊冷壁之溫度梯度較平緩，從圖 4-2 局部紐賽爾數也可以觀 察出目前結果與文獻也相當吻合，左冷壁面(a)的數值比右冷壁面(b)來的高，證實 了冷壁靠近熱壁時熱傳確實會提高，與文獻比較可以發現左右角落的局部紐賽爾數 差異較大，主要是由於下壁角落溫度不連續，而上壁角落點以均溫方式處理，所以 對數值處理時之計算方法會影響其結果之呈現。

17

4.2 菱形孔穴傾斜角度對自然對流之影響

透過前節之數值驗證，此一技術可用於不同傾斜壁面之夾角改變，對自然對流 之影響，本研究採用15 、30 、45 、60 、75 、90 、105 、120 、135 、150 ^o與165 等十一個角度進行比較，所採用之普朗特數皆為0.71，圖 4-3 至圖 4-6 分別繪製瑞 利數在10 、10 、10 與10 時不同傾斜角度之流線分布，從這四組圖可以看出當 瑞利數越大時對流現象越強因此浮力所造成之流動越貼近左右傾斜壁面，而孔穴中 央處則產生不同數目之迴漩渦流，其中瑞利數10 時不論菱形傾斜角度皆只有單一 個迴漩渦流，至於瑞利數10 時傾斜角介於60 至135 會有兩個迴漩渦流，但傾斜 角小於60 或大於135 則受壁面擠壓之效應減弱孔穴內對流，因此中央處也是只有 一個迴漩渦流，在瑞利數10 時傾斜角介於60 至150 其左右垂直壁面之對流更 強，因此有兩個迴漩渦流會更貼近左右壁面外，於中央處也產生第三個迴漩渦流，

且於傾斜角105 與120 時因左熱壁浮力效應與右冷壁下沉增強下並配合菱形壁之 擠壓，左上角及右下角之又各自再產生出另一小渦漩，形成內部五個渦流之結構，

至於在傾斜角介於30 與45 受菱形面擠壓，其左右迴漩向中央靠近，原先中央處 之旋渦消失，形成只有兩個迴漩渦流之結構，至於15 與165 之夾角下孔穴內之對 流也只剩下一個漩渦，最後在瑞利數為10 時因對流效應強烈許多，不僅渦流更貼 近壁面，且內部皆有多個渦旋產生而流場結構複雜，在60 與150 時，左右兩邊之 渦漩因對流效應強烈而拉扯出來另一小渦旋，且在此壁面傾斜角大於90 之夾角情 況下，因配合浮力效應其小渦旋皆比傾斜角小於90 之情況下來的大，而在傾斜角 介於30 與45 受菱形面擠壓，流場內部渦旋漸漸往中央靠近，導致中央處渦旋消 失，至於15 與165 之夾角下孔穴內，因瑞利數為10 ，對流效應較強烈，因此即 使左右壁之壓制下，其內部流場渦旋數目仍會超過一個以上的現象。

圖 4-7 至圖 4-10 分別顯示瑞利數在10 、10 、10 與10 下不同傾斜角度之溫 度分布，從這四組圖可以看出當瑞利數越大時對流現象越強因此造成中央處溫度水

18

平分層之現象越來越明顯，壁面邊界層也越來越薄，至於壁面傾斜角度大或小時壁 面摩擦對於流場之抑制造成溫度場的分布會較均勻之傳導結構。從瑞利數10 的溫 度分布圖來看，當傾斜角於15 時溫度梯度主要集中在左上角，隨著傾斜角度趨近 於90 時越向左下角集中，到了傾斜角大於90 時溫度梯度皆集中於左下角，至於水 平分層現象傾斜角大於90 的水平分層現象皆較傾斜角小於90 得來的明顯，在瑞利 數10 時與瑞利數10 時的溫度分布現象大致上一致，但當瑞利數提高時對於流場熱 傳之進行會有所助益。從圖 4-8 中不同壁面傾斜角圖中，可以看出瑞利數10 於傾 斜角度45 與150 時水平分層的現象依然很明顯，反觀圖 4-7 瑞利數10 ，傾斜角度 60 與150 時水平分層的現象已經很不明顯，由圖 4-9 之瑞利數等於10 之溫度分 布，其對流效應強以至於夾角由 30^o至165 下溫度水平分層現象皆很明顯，其溫度 邊界層較先前低瑞利數的圖形要薄。圖 4-10 為瑞利數等於10 時之溫度分布，因對 流效應更為強烈，導致在壁面傾斜角90 溫度斷面在中央處幾乎完全水平分層，在 壁面傾斜角15 與165 之溫度分布上，中央處仍有局部水平分層之現象，顯示在此 瑞利數下，溫度受對流效應仍有局部效果。

圖 4-11 為在不同瑞利數下平均紐賽爾數分布圖，其中可以看出在低瑞利數10 時傾斜角的大的比小的平均紐賽爾數皆只有略高一點(大約 5%，詳細如表二所列)，

但隨著瑞利數增加到10 時傾斜角大高於傾斜角小的趨勢就越來越明顯(增加到 21%)，雖然菱形壁面傾斜角在大夾角與小夾角下會對於流場有壓制的效應，但傾斜 角大於90 的配置相較於小於90 而言浮力效應在熱下冷上之有限空間其熱傳仍會 較佳。

19

4.3 同心圓柱之自然對流流場

本節探討同心圓柱在不同瑞利數下之自然對流流場，首先透過文獻所提供之溫 度斷面與平均紐塞數(K )來進行程式之驗證，之後再針對瑞利數10 ~10 ，的流 場穩態歷程來觀察瑞利數對流場結構與熱傳增益之影響，圖 4-12 為在Ra = 5 × 10 同心圓柱由0 至180 角度每格30 繪出溫度隨半徑之分布圖，與文獻之結果比較，

可以看出其分布很接近，在角度0 與30 時之溫度分布略低於文獻[4]之模擬結果，

但在角度120 、150 與180 時之模擬結果則略高於文獻之結果，圖 4-13 則為內圓 柱表面與外圓柱表面之局部紐塞數分布圖，與文獻之結果比較也很近似，唯在內圓 柱底部( = 150 ~180 )其值略低於文獻之結果。接著再進行不同瑞利數下對流場 結構與熱傳增益之探討，圖 4-14 為瑞利數為10 時之流線與溫度圖分布圖，由於10 時對流效應很弱導致流場內只有一渦旋且渦漩中心靠近水平中央處，且在內外圓柱 上之溫度邊界層較厚，溫度分布許多區域仍維持線性傳導之特徴，並未開始捲繞。

圖 4-15 為瑞利數在Ra = 5 × 10 時之流線、溫度分布圖，隨著瑞利數的提高，流場 浮力效應變強，渦旋由中央位置慢慢向上移動，流場內雖只有一渦旋之結構，但溫 度分布上已有明顯捲繞發生，靠近外圓上半段已有浮力效應導致之蕈狀溫度結構，

表示溫度傳遞已受對流效應主導之趨勢。圖 4-16 為瑞利數提高到10 時之流線、溫 度分布圖，流場內部因為對流效應的作用之下渦漩開始貼近外圓壁面，且內部之渦 漩受到內圓熱壁上升與外圓冷壁下沉增強之影響拉扯產生出另一小渦流之結構，溫 度上可以看出蕈狀結構變得更加明顯，且深入到圓柱下半段，在內圓柱底部區域之 溫度邊界層開始變得較薄。圖 4-17 為瑞利數在5 × 10 下之流線、溫度分布圖，因 對流效應更加強烈，導致速度邊界層貼向外圓壁面，上方之渦漩上移變得更加扁 平，並下側分裂成另一個漩渦結構，在此主渦流之下方，於圓柱底部開始產生局部 兩個小渦流結構，形成內部有四個渦流之流場結構，瑞利數在5 × 10 相較瑞利數 1 × 10 時之溫度邊界層更薄，尤其是在外圓壁面頂部附近，至於內圓柱上方之蕈 狀結構更加細長，顯示溫度衝擊於外圓上壁熱面越來越大。最後圖 4-18 為瑞利數

20

1 × 10 時之流線、溫度分布圖，與瑞利數在5 × 10 之流線、溫度圖做一比較其流 場結構與溫度分布相當類似，雖然在左上之渦旋結構變得更加扁平細長，其下方拉 扯出的三個小渦流結構更加貼近外圓壁面，且略為變大外，此三個渦流附近之流線 已演變得更加扭曲的現象，因此欲計算更大的瑞利數將會需要更多的網格以處理此 一高度扭曲之流場，至於在溫度分布上，因內圓熱壁上升與外圓冷壁下沉增強導致 溫度貼近兩圓壁面，相較5 × 10 時之溫度分布，溫度結構變得更加扁平細長之薄 熱邊界層外，其中間分布上又多一些扭動的特徴。

圖 4-19(a)、(b)分別為不同瑞利數下內圓與外圓局部紐塞數分布圖，圖中紐塞 數分布與流線與溫度分布之現象一致，隨著瑞利數之增加其溫度邊界層變薄之發展 而增加其內外圓柱表面之紐塞數，其中在內圓柱對流體加熱之區域，主要來至下半 側，亦即 = 90 ~180 ，反觀外圓柱之散熱僅集中在頂面附近亦即 = 0 ~30 ， 其噴流衝擊點 = 0 會最大，由於內圓柱底側之小渦流之低速流動的效果，形成底 部溫度分層之現象，造成外圓壁面在 > 140 其散熱量明顯變得很微小。

4.4 動態圓柱之自然對流流場

針對暫態熱傳之分析探討可以透過前一問題之設定，並將內圓柱改成上下週期 之運動，並利用文獻[10]之數值結果，比對目前程式在暫態分析之能力，所採用之 圓柱擺動之振幅(ε = 0.5)與頻率(Ω = 1)與文獻中所述相同，初始條件流速為 零，溫度為線性分布，圖 4-20 為目前計算之前三個週期之內外圓柱平均紐塞數發 展歷程，相較文獻之解答其趨勢接近，但目前分析之內圓最大熱傳發生在 0.64T 與 文獻 0.56T 有較延後之差異，而外圓表面散熱最大發生之瞬間目前發生在 0.98T，

而文獻中圖檔估計也是在 1T 附近，至於目前分析之內圓最小熱傳發生在 0.23T 與 文獻 0.16T 有較延後之差異，而外圓表面散熱最小發生之瞬間目前發生在 0.35T，

而文獻中圖檔估計在 0.28T 附近，由這些比較顯示目前之內外圓柱之加熱散熱的峰 值略有不同外且較延後一些。圖 4-21 至 4-24 則繪製最後一個週期中每四分之一週

21

期瞬間之溫度與流線分布，由此四張系列發展圖之比較，與文獻結果是相當雷同 的，其內部渦流受內圓柱之擠壓與浮力對流之交互影響，所以漩渦中心位置與強弱 會隨時間改變，由瞬間溫度分布可觀察出 0.25T 時溫度邊界層較厚且無明顯蕈狀結 構，與先前加熱與散熱最低瞬間(0.23T/0.35T)相符合，至於內圓柱在 0.5T 時以向下 運動，因此下側之溫度邊界層會變較薄，再辦隨內外圓柱間距之減小，因此內圓柱 之最大熱傳會發生在 0.64T 瞬間，至於外圓柱之溫度邊界層顯示在 0T 瞬間頂面最 薄，乃因浮力對流與圓柱上移之加總效果，也證明目前外圓柱最大熱傳發生於 0.98T 瞬間。最後圖 4-25 為一個週期中四個瞬間之內外圓柱表面紐塞數分布，當內圓柱 位於中央開始向上移動時(0T)，其內圓柱最大熱傳量發生於135 ，主要是因為圓柱 上移對內圓正下方產生真空作用所以熱傳量會變弱一些，配合浮力所產生之對流作 用其邊界層最薄之位置會向下游移動，至於 0.25T 瞬間之紐塞數如同前述之流場結 構，相較其他瞬間其量為最小，在 0.5T 瞬間之紐塞數因圓柱向下運動因此整體值 皆比 0T 要高峰值也是在135 附近，最後在 0.75T 瞬間之紐塞數因此時圓柱瞬間靜 止在下死點，此時空間窄產生之擠壓加浮力對流，氣流較強因此邊界層較薄處會後 移至100 處。對於外圓柱局部紐塞數分布圖在此四個瞬間趨勢都雷同，主要集中於 0 頂面處，因圓柱向上運動瞬間，運動方向與浮力方向一致，有助於外圓散熱，所 以 0T 瞬間會最高，反觀 0.25T 瞬間當圓柱靜止於上死點，空間太窄造成對流渦流 被擠開外移，因此0 頂面處之紐塞數最低，但是在 30 與 50 區間會較其他瞬間要 強一些，至於在 140 ~180 外圓柱表面其穩態紐塞數皆很低，對於動態分析上，0.75T 時因空間之變窄會增加此處之熱傳量，透過上述之局部紐塞數之觀察可以明確整合 流場與熱傳之發展，提供暫態熱流分析之可能。

22

4.5 不同上下游邊界設定對流場與熱傳之影響

對於不同上下游邊界設定對外流場之影響上，採用先前學長[15]所使用之邊界 條件，上游入口處之壓力為一大氣壓，於開放入口邊界上採用 Potential flow 之設定 其流速，而出口邊界則是以流線方向進行外插，並進行不同上下游5D, 10D, 15D, 20D 四個不同的距離進行對流場與熱傳之影響，圖 4-24 為Re = 20參考文獻[28]所提供 之溫度與渦度分布圖與目前四個不同上下游距離設定下之溫度與 vorticity 比較圖，

由這些圖形之比較可發現，雖然上下游之邊界不同對圓柱附近之流場與溫度分布影 響不大，所有圖形皆相當雷同。圖 4-25 為不同上下游距離設定時之尾流長度比較，

其中我們可以發現流線分布上，不論是分離位置上與尾流長度上也幾乎沒有太大差 異。若將溫度分布以更大範圍繪出，並重疊四個上下游距離之結果，圖 4-26 為其 整體溫度之比較，在此圖中則能較明顯看出溫度上略有不同，主要原因是在於上下 游距離太短導致垂直方向之限制，所以5D之設定其溫度擴散之張角是略小的，反 觀當上下游距離為10D, 15D, 20D可以觀察出其最外側溫度擴散之張角是較為一致 的，而在下游處其溫度尾流擴散之長度也是略有不同的，採用5D之上下游是最短 的，以20D之上下游是是最長但。對於阻力係數之預測上不同上下游設定也是略有 影響的，表三為在不同上下游距離下之阻力係數與文獻比較表，表中在上下游距離 5D時之阻力係數為2.35，與各文獻阻力係數平均值相比約有 8%誤差，當上下游距 離達到20D左右時其阻力係數變動值就小於1%，圖 4-27 為Re = 20在不同上下游設 定下，局部紐塞數之分布比較圖，在上下游距離5D時，其局部紐塞數相較於文獻 與其它不同上下游距離的差異較大，主要原因是在於上下游距離太短導致出口邊界 影響到內部溫度分布所造成，但在其餘不同上下游距離與文獻比較上並沒有太大差 異(< 2%)。由前述之結果顯示，本論文提出以流線方向進行外差的方式處理邊界 條件雖然可行，但仍要以10D以上之距離對分析結果的影響會較小。

23

第五章 結論

本研究透過數值程式之分析進行了非正交網格之驗證得以進行稜形夾角改變 對孔穴內自然對流之影響，接著以同心圓柱熱傳問題觀察瑞利數增加下內部流場結 構之變化與熱傳增益之情形，進而處理內部圓柱於上下週期運動下之動態熱傳分 析，最後進行開放圓柱之熱傳分析。由前面之結果可得下述結論：

1. 當稜形夾角在 30^o下自然熱傳問題中，本研究之結果與文獻所提供之紐塞數、

溫度分布與流線分布皆相當一致。

2. 在相同瑞利數下改變稜形之夾角，其最大熱傳量發生在 90^o，不受左右壁面擠 壓之造型下，若相同夾角則瑞利數越大平均紐塞數越高，另外流場內部之漩渦 數目會隨對流效應之增強而變多，因此當夾角縮小至 30^o或是放大至 165^o時其 左右壁面之摩擦增強，其內部只剩一個迴旋。

3. 針對同心圓之自然熱傳問題也利用文獻之數據進行紐塞數、溫度分布與流線分 布之比對，並進行不同瑞利數之探討，當瑞利數小於10 其流場結構左右各一 大迴旋，當瑞利數於5 × 10 則左右各一大迴旋內部開始分裂成為兩個漩渦，

若瑞利數增至1 × 10 則同心圓下方產生另一迴流且內部有兩個漩渦，並將原 先一個大迴旋內部兩個漩渦之結構向上擠壓，形成更複雜之流場現象，此一渦 流與迴旋數目增多之現象可以說明層流開始轉換至不穩定流動之發展歷程。

4. 當內部圓柱處於上下週期運動時，目前分析程式以雙時間步階方式進行非穩態 之計算，透過文獻之結果與目前分析結果在瞬間流線與溫度分布很相近，但在 瞬間最大與最小之熱傳發展上有較延後之差異。

5. 最後在開放圓柱之熱傳分析上，針對上下游之邊界設定上至少需有 10 直徑對 局部紐塞數之影響將不會超過 2%。

24

參考文獻

[1]W.F. Noh, “A Time-dependent Two-space-dimensional Coupled Eulerian- Lagrangian Code”, Methods in Computational Physics, vol.3, p.117, B. Alder, S.

Fernbach and M. Rotenberg, Academic Press, New York, 1964.

[2]C.W. Hirt, A.A. Amsden, and H.K. Cooks, “An Arbitrary Lagrangian -Eulerian Computing Method for all Flow Speeds,” Journal of Computational Physics, vol.14, pp.227-253, 1974.

[3]Thomas J.R. Hughes, Wing Kam Liu, Thomas K. Zimmermann“Lagrangian-Eulerian finite element formulation for incompressible Viscous flows”Computer Methods in Applied Mechanics and Engineering, Volume 29, Issue 3, December 1981, PP.329-349 [4]T.H. Kuehn, R.J. Goldstein, “An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74 (1976) 695–719.

[5]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.

[6]C-H CHENG, J-L HONG, “Numerical prediction of lock-on effect on convective heat transfer from a transversely oscillating circular cylinder,”International Journal of Heat and Mass Transfer 40 (1997) 1825–1834.

[7]J.S Yoo, “Prandtl number effect on bifurcation and dual solution in natural

convection in a horizontal annulus,”International Journal of Heat and Mass Transfer 42 (1999) 3279–3290.

25

[8]C.H. Cheng and K.S. Hung, “ Numerical Predictions of Flow and Thermal Fields in a Reciprocating Piston-Cylinder Assembly, ”Numerical Heat Transfer, Part A, Vol. 38 (2000) pp. 397-421.

[9]W.S. Fu, B.H Tong, “Numerical investigation of heat transfer from a heated

oscillating cylinder in a cross flow,” International Journal of Heat and Mass Transfer 45 (2002) 3033–3043.

[10]Esam M. Alawadhi, “Natural convection flow in a horizontal annulus with an oscillating inner cylinder using Lagrangian-Eulerian kinematics”Computers and fluids 37 (2008) 1253–1261.

[11]D. Angeli, P. Levoni, GS. Barozzi, “Numerical Predictions for stable buoyant regimes within a square cavity containing a heated horizontal cylinder,” International Journal of Heat and Mass Transfer 51 (2008) pp.553–565.

[12]Y.Y. Lung, Y.C. Lin., “Numerical investigation of natural convection in an inclined wavy cavity using a high-order differential scheme,” Journal of Aeronautics,

Astronautics and Aviation, Series A, Vol. 43, 2011, No.1, pp. 045-056.

[13]H. Saleha, R. Roslanb, I. Hashima.“Natural convection heat transfer in a

nanofluid-filled trapezoidal enclosure.”International Journal of Heat and Mass Transfer Volume 54, Issues 1–3, 15 January 2011, PP. 194–201.

[14]W.S Fu, J.C Huang, C.G Li., “Enhancement of forced convection heat transfer in a three-dimensional laminar channel flow with insertion of a moving block,” International Journal of Heat and Mass Transfer, Vol. 53, 2010, pp.3887-3897.

[15]廖經勝., “具波形壁之二維渠道於可壓縮流場下強制熱傳分析,” 中華大學碩士

26

論文, 2011-7 月.

[16]R. Anandalakshmi, Tanmay Basak., “Heat flow visualization for natural convection in rhombic enclosures due to isothermal and non-isothermal heating at the bottom wall,”

International Journal of Heat and Mass Transfer, Vol. 55, 2012, pp.1325-1342.

[17]C.C. Liao, C.A. Lin.,“Influences of a confined elliptic cylinder at different aspect ratios and inlinations on the laminar natural and mixed convection flows,” International Journal of Heat and Mass Transfer, Vol. 55, 2012, pp. 6638-6650.

[18]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.

[19]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.

[20]Irons, B. M., “Engineering applications of numerical integration in stiffness methods,”AIAA Journal, Vol. 4, No. 11, 1966, pp. 2035-2037.

[21]M. Hortmann, M. Peric, and G. Scheuerer,“Volume Multigrid Prediction of Laminar Natural Convection; Bench-Mark Solutions”, Int. J. Numerical Methods in Fluids, 11 (1990), 187-207.

[22]A.Nag, A. Sarkar, and V. M.K. Sastri, “Natural Convection in a Differentially Heated Square Cavity with a Horizontal Partition Plate on the Hot Wall”, Comput.

Methods Appl. Mech. Eng., 110(1993), 143-156.

[23]X. Shi and J. M. Khodadadi, “Laminar Natural Convection Heat Transfer in a

Differentially Heated Square Cavity Due to a Thin Fin on the Hot Wall”, Journal of Heat Transfer, 125 (2003),624-634.

[24]E. Bilgen, “Natural convection in cavities with a thin fin on the hot wall”, Int. J.

Heat and Mass Transfer, 48 (2005), 3493-3505.

27

[25] Z. Wang, J. Fan, K. Luo, K. Cen, Immersed boundary method for the simulation of flows with heat transfer, Int. J. Heat Mass Transfer. 52 (2009) 4510–4518.

[26] D.J. Triton, Experiments on the flow past a circular cylinder at low Reynolds number, J. Fluid Mech. 6 (1959) 547–567.

[27] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian gird method for viscous imcompressible flows with complex immersed boundaries, J. Comput. Phys.

156 (1999) 209–240.

[28]C.C. Liao, C.A. Lin.,“Simulations of natural and forced convection flows with moving embedded object using immersed boundary method ,” Comput. Methods Appl.

Mech. Engrg. 213–216 (2012) 58–70.

[29]N. Zhang, Z.C. Zheng, S. Eckels, Study of heat-transfer on the surface of a circular cylinder in flow using an immersed-boundary method, Int. J. HeatFluid Flow 29 (2008) 1558–1566.

[30] R.P. Bharti, R.P. Chhabra, V. Eswaran, A numerical study of the steady forced convection heat transfer from an unconfied circular cylinder, Heat Mass Transfer. 43 (2007) 639–648

28

表一、壁面傾角 90^o下之平均紐塞爾數 Reference Average Nusselt number, Nu

Ra = 10 Ra = 10 Ra = 10 Ra = 10 Davish[5] 2.234 4.51 8.789 - Hortmann[21] 2.4468 5.5231 8.8359 -

Nag[22] 2.24 4.51 8.82 -

Shi[23] 2.247 4.532 8.839 16.935 Bilgen[24] 2.245 4.521 8.8 16.629

Present 2.25 4.52 8.83 16.58

表二、 不同壁面傾角下之平均紐塞爾數 Angle(α) Average Nusselt number,

Ra=10⁴ Ra=10⁵ Ra=10⁶ Ra=10⁷

15⁰ 1.13 1.31 2.87 7.48

30⁰ 1.18 2.43 5.75 12.16

45⁰ 1.53 3.45 7.32 14.44

60⁰ 1.90 4.06 8.21 15.71

75⁰ 2.14 4.40 8.67 16.38

90⁰ 2.25 4.52 8.83 16.58

105⁰ 2.23 4.47 8.72 16.39

120⁰ 2.09 4.23 8.34 15.77

135⁰ 1.78 3.78 7.63 14.61

150⁰ 1.31 2.99 6.41 12.71

165⁰ 1.19 1.69 4.10 9.10

表三、 不同瑞利數下內外圓柱之平均紐塞爾數

Ra = 10 Ra = 5 × 10 Ra = 10 Ra = 5 × 10 Ra = 10 Inner 3.31 4.88 5.81 8.44 9.88

Outer 1.27 1.88 2.23 3.23 3.79

29

表四、 在Re = 20、不同上下游距離之阻力係數與文獻比較圖

Present Liao and Lin[28] Triton [26] Ye et al. [27] Wang et al. [25]

5D 2.35

2.21 2.22 2.03 2.18

10D 2.12 20D 2.07 30D 2.05

30

(a)流線圖(Anandalakshmi and Basak)

(b)流線圖(Present)

(c)溫度圖(Anandalakshmi and Basak) (d)溫度圖(Present)

圖 4-1 傾斜角度30 、Ra = 10 、Pr = 1000與文獻[16]流線溫度之比較圖

(a)Left Wall (b)Right Wall

(c)Bottom wall

圖 4-2 在傾斜角度30 、Ra = 10 、Pr = 1000與文獻[16]局部紐塞數比較圖

Distance

LocalNusseltNumber

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10 12 14 16 18 20 22 24

Present

Anandalakshmi and Basak

Distance

LocalNusseltNumber

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10 12 14 16 18 20 22 24

Present

Anandalakshmi and Basak

Distance

LocalNusseltNumber

0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10 12 14 16 18 20 22 24

Present

Anandalakshmi and Basak

31

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰

(F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-3 在Ra = 10 下流線分布圖

32

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-4 在Ra = 10 下流線分布圖

33

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-5 在Ra = 10 下流線分布圖

34

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-6 在Ra = 10 下流線分布圖

35

(A)15⁰ (B)30⁰

(C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-7 在Ra = 10 下溫度分布圖

0.9 0.8

0.7 0.6 0.5

0.4 0.3

0.2 0.1

36

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-8 在Ra = 10 之溫度分布圖

0.9 0.8

0.7 0.6

0.5 0.4 0.3

0.2 0.1

37

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-9 在Ra = 10 之溫度分布圖

0.9

0.8

0.7 0.6 0.5 0.4 0.3 0.2

0.1

38

(A)15⁰ (B)30⁰ (C)45⁰

(D)60⁰ (E)75⁰ (F)90⁰

(G)105⁰ (H)120⁰ (I)135⁰

(J)150⁰ (K)165⁰

圖 4-10 在Ra = 10 之溫度分布圖

0.9

0.8

0.7 0.6 0.5 0.4 0.3 0.2

0.1

39

圖 4-11 不同傾斜角度與瑞利數下平均紐塞數

 N u

0 15 30 45 60 75 90 105 120 135 150 165 180 0

2 4 6 8 10 12 14 16 18 20

Ra=10⁷ Ra=10⁶ Ra=10⁵ Ra=10⁴

40

圖 4-12 在Ra = 5 × 10 同心圓柱不同角度處溫度與文獻[4]比較圖

圖 4-13 在Ra = 5 × 10 局部紐塞數與文獻[4]比較圖

R

T

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1















 keq

0 30 60 90 120 150 180

0 1 2 3 4 5 6 7 8 9 10 11

Kuehn and Goldstein(Outer) Kuehn and Goldstein (Inner) Present(Inner)

Present(outer)

41

圖 4-14 在Ra = 1 × 10 下流線(左)與溫度(右)分布圖

圖 4-15 在Ra = 5 × 10 下流線(左)與溫度(右)分布圖

42

圖 4-16 在Ra = 1 × 10 下流線(左)與溫度(右)分布圖

圖 4-17 在Ra = 5 × 10 下流線(左)與溫度(右)分布圖

43

圖 4-18 在Ra = 1 × 10 下流線(左)與溫度(右)分布圖

44

(a)Inner

(b)Outer

圖 4-19 在不同瑞利數下外圓局部紐塞數分布圖



LocalNusseltNumber

0 30 60 90 120 150 180

0 2 4 6 8 10 12 14

16 Ra=1X10⁴

Ra=5X10⁴ Ra=1X10⁵ Ra=5X10⁵ Ra=1X10⁶



LocalNusseltNumber

0 30 60 90 120 150 180

0 5 10 15 20 25

Ra=1X10⁴ Ra=5X10⁴ Ra=1X10⁵ Ra=5X10⁵ Ra=1X10⁶

45

(a)Inner

(b)Outer

圖 4-20 瑞利數Ra = 5 × 10 、 = 0.5、Ω = 1，內外圓柱平均紐塞數發展歷

Dimensionless Time AvreageNusseltNumber(keq)

0.5 1 1.5 2 2.5 3

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Present Alawadhi.

Dimensionless Time AverageNusseltNumber(keq)

0.5 1 1.5 2 2.5 3

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Present Alawadhi.

46

程

圖 4-21 t = 0T時之溫度(上)、流線(下)分布與文獻(左)之比較圖

47

圖 4-22 t = T時之溫度(上)、流線(下)分布與文獻(左)之比較圖

48

圖 4-23 t = T時之溫度(上)、流線(下)分布與文獻(左)之比較圖

49

圖 4-24 t = T時之溫度(上)、流線(下)分布與文獻(左)之比較圖

50

(a)Inner

(a)Outer

圖 4-25 不同週期時局部紐賽數分布圖

 LocalNusseltNumber(keq)

0 30 60 90 120 150 180

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

0T T/4 T/2 3T/4

 LocalNusseltNumber(keq)

0 30 60 90 120 150 180

2 4 6 8 10

0T T/4 T/2 3T/4

51

(a)

(b)5D (c)10D

(d)15D (e)20D

圖 4-26 在Re = 20、不同上下游大小下，溫度與 vorticity 和文獻[28](a)比較圖

52

(a)5D

(b)10D

(c)15D

(d)20D

圖 4-27 在Re = 20、不同上下游距離設定下圓柱流線圖