水平同心圓柱間自然對流之分歧現象

國立交通大學

機械工程研究所

碩 士 論 文

水平同心圓柱間自然對流之分歧現象

Bifurcation Phenomenon of Natural Convection Between

Horizontal Concentric Cylindrical Annulus

研 究 生：林 紀 元

指導教授：楊 文 美 博士

中華民國九十三年七月

水平同心圓柱間自然對流之分歧現象

Bifurcation Phenomenon of Natural Convection Between

Horizontal Concentric Cylindrical Annulus

研 究 生：林 紀 元 Student：Ji-Yuan Lin 指導教授：楊 文 美 Advisor：Wen-Mei Yang 國 立 交 通 大 學 機 械 工 程 研 究 所 碩 士 論 文

A Thesis

Submitted to Instuitute of Mechanical Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of The Requirments

for the Degree of

Master of Science

in

Mechanical Engineering

June 2004

Hsinchu, Taiwan, Republic of China

摘 要

本文使用數值方法探討兩水平同心圓柱間因溫度差所引發的自 然對流現象，及圓柱間流場在不同參數下的變化。文中首先建立描 述圓柱流場及溫度場的統御方程式，並將其無因次化。在數值方法上 ，先以 Chebyshev polynomials 及 Fourier series 對空間座標作雙重級 數 展 開 ， 並 將 無 因 次 化 後 之 統 御 方 程 式 中 非 線 性 項 利 用 Adam- Bashforth Method 作離散化，而線性項則採用 Crank-Nicolson Method 處理。最後將雙重級數展開式代入離散方程式再以 Collocation Method 與 Galerkin Approximation 將方程式轉換為代數方程式。由數值解中可 探討受普郎特數與半徑比兩參數影響的兩水平同心圓柱間流場與溫 度場間之行為。 對於水平同心圓柱間自然對流流場產生很多不同流場型態與轉 換現象，均取決於普郎特數與半徑比的變化。 本文中以兩種不同半徑比，分別改變普郎特數來觀察穩定流場由 單一渦漩流（crescent-shaped eddy）轉變為在圓柱上方分離成兩個相 反方向旋轉的渦漩流（two counter-rotating eddies）：當圓柱半徑為 0.2 ，隨著普郎特數增加，液動不穩定性及溫度不穩定性皆會增大，流 場會發生分歧現象的臨界瑞里數愈大，所計算出圓柱間的平均紐賽 數也隨著瑞里數增加而變大；當圓柱半徑比為 0.5，在 0.3≦Pr≦0.5 ，自然對流流場中內熱圓柱會產生二次流的現象，隨著普郎特數增加 ，由於圓柱間寬距變小，液動不穩定性會減小，溫度不穩定性會增大 ，流場會發生分歧現象的臨界瑞里數略為減小，所計算出圓柱間的 平均紐賽數也隨著瑞里數增加而略為減小。 i

誌 謝

一轉眼，交大碩士班生涯即將結束，在此感謝指導教授楊文美博 士的悉心指導與撰寫論文時字字句句的斧正，並不厭其煩的改正初稿 不盡完善之處，俾使本論文得以順利完成，特別向楊老師致上最誠摯 的謝意。 同時也感謝論文口試委員游明輝教授與盧定昶教授，在百忙之中 撥空審閱本論文，並於口試期間，給予論文指正與寶貴建議，使得本 論文更趨完善與嚴謹。 感謝博士班學長致廣、又永，碩士班學長志宏、光偉、朝任、俊 良、文忠、耀鋒，同學志凱、彥翔、正華、志銘、至恩、鈞賀、明諭， 以及學弟盈立、勤暐、昱志、勝文、宏文、宏仁、英棋、友約於課業 上及生活上互相的提攜照顧與關懷，使我在碩士班生涯擁有相當多的 回憶。 更感謝父母親多年來含莘如苦的照顧我和弟弟─紀綱，讓我無後顧 之憂地研究專業領域上的知識，也感謝弟弟在課業上的鼓勵，以及一 直陪伴在我身邊支持鼓勵我的女友─曉玲，最後，將本論文成果獻給我 親愛的家人！ ii

目 錄

摘要………..……….i 誌謝………...………...……ii 目錄……….iii 圖目錄………..……….….…..v 符號說明………...…..…………...…...….vii 第一章 緒論………..1 1. 1 文獻回顧……….….1 1. 2 研究目的………...………...5 第二章 數學模式………..7 2. 1 統御方程式………..7 2. 2 無因次化………..9 第三章 數值方法………11 3. 1 求取數值解的方法………11 3. 1. 1 雙重級數展開………...…...……….…….11 3. 1. 2 離散方程式…………...………...………….……….12 3. 1. 3 代數方程式…………...……….14 3. 2 數值解的分析方法………18 3. 2. 1 流線圖………18 3. 2. 2 計算平均紐賽數(Nu)………18 iii

第四章 結果與討論………..………..……....21 4. 1 準確度分析………21 4. 2 臨界瑞里數對圓柱間自然對流流場之影響………22 4. 3 探討普郎特數對自然對流流場之影響………22 4. 3. 1 半徑比η=0.2 時，普郎特數對自然對流流場之影響… ………22 4. 3. 2 半徑比η=0.5 時，普郎特數對自然對流流場之影響… ………24 4. 4 探討圓柱半徑比對自然對流流場的影響………26 第五章 結論………28 參考文獻………....30 iv

圖 目 錄

圖一 基本物理模式示意圖…………..………..…..33 圖二 η=0.2，Pr=0.5，Ra=5000 時，測試各展開項數所計算出的平 均紐賽數，並與Yoo（1999）比較。………...….34 圖三 η=0.2，Pr=0.3，Ra=4000 時，測試各展開項數所計算出的平 均紐賽數，並與Yoo（1999）比較。………35 圖四 η=0.2，Pr=0.3 時，每個瑞里數所計算出的平均紐賽數，並與 Yoo（1999）比較。……….………36 圖五 η=0.2，Pr=1.0 時，每個瑞里數所計算出的平均紐賽數，並與 Yoo（1999）比較。……….37 圖六 η=0.5，Pr=0.7 時，測試各展開項數所測的臨界瑞里數，並與 Yoo（1996）比較。………...38 圖七 當η=0.2，Pr=0.3 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=2000（c） Ra=2500（d）Ra=2700（e）Ra=2800（f）Ra=3600（g） Ra=4000………39 圖八 當η=0.2，Pr=0.5 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=3000（c） Ra=4000（d）Ra=6000（e）Ra=6800（f）Ra=6900（g） Ra=8000……….42 圖九 當η=0.2，Pr=0.7 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=3000（c） Ra=4000（d）Ra=6000（e）Ra=8000（f）Ra=9000（g） Ra=9500………45 圖十 當η=0.5，Pr=0.3 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1500（b）Ra=1600（c） v

Ra=1900（d）Ra=2000（e）Ra=2900（f）Ra=3000（g） Ra=3100……….48 圖十一 當η=0.5，Pr=0.5 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=2000（c） Ra=2900（d）Ra=3000（e）Ra=3500（f）Ra=4000（g） Ra=6000……….51 圖十二 當η=0.5，Pr=0.7 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=2000（b）Ra=2500（c） Ra=2800（d）Ra=2900（e）Ra=3500（f）Ra=4000（g） Ra=5000……….54 圖十三 當 Pr=0.3 時，不同瑞里數所計算出的平均紐賽數。……….57 圖十四 當 Pr=0.5 時，不同瑞里數所計算出的平均紐賽數。……….58 圖十五 當 Pr=0.7 時，不同瑞里數所計算出的平均紐賽數。……….59 圖十六 當η=0.2 至η=0.5，臨界瑞里數與 Pr 之關係。……….60 vi

符 號 說 明

符 號 A d 展弦比(aspect ratio) 內外圓柱半徑差 g 重力加速度 Gr 格雷修夫數(Grashof number) L Nu Nu 寬距 紐賽數(Nusselt number)

平均紐賽數(mean Nusselt number)

p 壓力 P 無因次壓力 Pr 普郎特數(Prandtl number) r,θ 圓柱座標 R 圓柱半徑 Ra Rac 瑞里數(Rayleigh number)

臨界瑞里數(critical Rayleigh number)

T 溫度 Tn n 階 Chebyshev 多項式 t 時間 u,v 速度分量 U,V 無因次速度分量 X,θ 無因次座標 τ 無因次時間

Θ

無因次溫度 η 內外圓柱半徑比 ξ Chebyshev 多項式座標 ν 運動黏滯係數 vii

ρ 密度 k n φ 熱傳導係數 速度與溫度分量之基底函數 ψ α 流線函數(stream function) 熱擴散係數 β 熱膨脹係數 △ 2 2 2 2 2 X 1 X X 1 X ∂θ ∂ + ∂ ∂ + ∂ ∂ ▽2 2 2 2 2 2

r

1

r

r

1

r

∂

θ

∂

+

∂

∂

+

∂

∂

上 標 i 第 i 個時間間隔 平均值 下 標 1 2 內圓柱 外圓柱 c 臨界值 viii

第 一 章

緒 論

1.1 文獻回顧

多年來有關水平同心圓柱間的自然對流，研究方向由早期大部份 集中於熱交換器到現今在工業上的應用相當地廣泛，諸如：電子零件 的冷卻、材料製程的熱傳分析問題等。水平同心圓柱間的自然對流通 常發生在內圓柱為高溫，而外圓柱為低溫時，甚至在非常小的溫度差 異量下也會產生對流現象。 Beckmann (1931) 首先研究水平同心等溫圓柱間的自然對流， 在圓柱間研究空氣(air)、氫(hydrogen)和二氧化碳(carbon dioxide)並且求得 總熱傳係數(overall heat transfer coefficient)。 Voigt and Krischer (1932) 與Beckmann (1931) 也有相似的研究，如在普郎特數(Prandtl number) 很大時，圓柱間使用水(water)、transformer oil 和機油(machine oil)。

Crawford and Lemlich (1962) 第一個以數值方法研究水平同心等 溫圓柱間的自然對流，當時使用方法為高斯替代反覆逼近法 (Gauss Seidel iterative approach)。 Abbott (1964) 在非常薄的同心圓柱間中使 用反矩陣(Matrix Inversion)的技巧求解。 Powe et al. (1969) 在充滿空氣的圓柱間中藉由實驗研究探討流體 分類的型態與流體轉換的數範圍，發現瑞里數(Rayleigh number)和展 弦比(Aspect ratio) A 的影響可以劃分出不同流體形態的流場。而定 義展弦比 A 為內圓柱直徑與寬距(gap width) L 的比值。而寬距的定 義為外圓柱半徑與內圓柱半徑的差值。由實驗結果可將流體形態以 不同的範圍來做以下的分類： 1

(一) 0<A<2.8：在圓柱間的上部範圍的流場渦捲(cell)會從原來的單一 渦捲引發出其他多渦捲(multicellular)流場(也可解釋成穩定流體 轉換到不穩定)，並且可發現縱軸方向流體呈現振盪(oscillatory) 特性。 (二) 在 2.8<A<8.5：在瑞里數大於臨界瑞里數時(Critical Rayleigh number)，圓柱間的上部範圍流場渦捲會由穩定轉換到不穩 定振盪的特性，且研究螺旋形流場運動的特性。 (三) 在 A>8.5：在圓柱間的上部範圍流場渦捲不穩定結果形成時， 會有2-D 多渦捲組成的對流現象。 Powe et al. (1971) 以數值方法考慮普郎特數為 0.7 時，當開始形成 逆向旋轉渦流(counter-rotating eddy)時，研究流體穩定(steady)轉換到流 體不穩定(unsteady flow)現象並測定其瑞里數值。

Kuehn and Goldstein (1976) 以實驗與數值研究在考慮普郎特數為 多少時的流體會影響速度與溫度分佈和局部(local)熱轉換係數，例如在 空氣(Pr=0.7)和水(Pr=6.0)的熱對流現象。他們在實驗方面使用干涉 計(Mach-zehnder interferometer)來決定溫度分佈與局部熱轉換係數。另 外，在數值方面考慮水或空氣的圓柱間先假設寬矩與內圓柱直徑的 比率為0.8，且瑞里數介於 2.11× 104~9.76× 105的範圍下使用有限差分 法(Finite Different Method)可求得其結果。此後，使得各種水平同心圓 柱間自然對流的流場流動現象之實驗與數值研究開始有了廣泛的發 展。

Mojtabi and Caltagirone (1979) 研究同心圓柱間自然對流 2-D 流場 的穩定性，使用擾亂法(Perturbation Method)得到解，但其假設條件為

限制在要觀察的流場圓柱間的上層區域範圍內，並假設空氣在瑞里 數值很小時與外圓半徑與內圓半徑之半徑比率(R=R2/R1)很小時，使用

級數展開(power series expansion)可求得其解。

Tsui and Tremblay (1984) 研究以數值分析並使用波斯尼克近似法 (Boussinesq approximation)與渦度－流線函數逼近法(Vorticity-stream function approach)，來探討水平同心等溫圓柱間自然對流熱傳之暫 態(transient)解，而數值結果說明以直徑比或半徑比(R=R2/R1)當參數， 求得符合圓柱間自然對流熱傳計算之紐賽數(Nusselt number)－格雷 修夫數(Grashof number)曲線圖。 Rao et al. (1985) 以數值方法研究同心圓柱間自然對流之穩定性， 並計算出在普郎特數 Pr＝5000 時且展弦比 A＝4.7 時，得到瑞里數 為大於 2700 的穩定 3-D 螺旋形流場。與 Dyko (1999) 在實驗與數 值研究上，當瑞里數大於臨界瑞里數時會發生由浮力所驅動(buoyancy driven)的流體現象，而理論預測在展弦比 A 大於 2 時，瑞里數大於 臨界瑞里數時會產生二次流(secondary flow)的現象並轉換到 3-D 擾 動流。

Choi and Kim (1993) 在研究同心圓柱間自然對流之流體為空氣時

探討 2-D 流場，發現正在逐漸成長的渦捲會有穩定的現象。在數值方

面假設很小的普郎特數下研究從流體穩定轉換到不穩定的過渡流體， 並利用時間分列式法(Time Marching Method)解 3-D 擾動(disturbance) 下的線性方程式。而在實驗方面當同心圓柱間自然對流之流體為空氣 時(Pr=0.71)，假設展弦比 A 的值介於 2.1~10 之間並且瑞里數在 103≦ Ra≦105 的範圍下當瑞里數大於臨界瑞里數時，2-D 基態流體(base flow)會呈現不穩定狀態這與 3-D 擾動有關。

Cadiou et al. (1998) 研究充滿空氣狹窄的水平圓柱間自然對流 的不穩定，並在很小的外圓半徑與內圓半徑之比率下討論熱不穩定與 液動不穩定現象，在低瑞里數(Ra≦3000)時發現圓柱間上部範圍會 有熱不穩定發生，在高瑞里數(heigh Rayleigh number)時且外圓半徑

與內圓半徑之比率為1.14 時會有液動不穩定發生。這時當瑞里數繼續 升高時會發生多渦捲流轉換到單渦捲流的反向轉換現象。 而近年，Yoo (1998,1999a,b,c) 探討水平圓柱間自然對流現象，顯 示在內熱圓柱出現多種變化的流體轉換現象與普郎特數和圓柱間的形 狀大小有著強烈的關係，並且以有限差分法亦對其結果做了驗證。 Yoo (1999a) 在水平狹窄的圓柱間考慮自然對流問題，當展弦比 A=12 時， 探討普郎特數的影響。而熱傳導範圍中發現流體的穩定性可分成兩 種現象：1.當 Pr≦0.2，在垂直圓柱間穩定與振盪流體會有複雜的變 化像旋轉的渦捲。2.當 Pr≧0.3，在圓柱的頂端範圍，發現一個正在 逐漸成長的渦流(eddy)現象。而判斷是否呈現不穩定流體與逐漸形成 渦捲式的對流及分歧現象(Bifurcation phenomenon)都證明與普郎特 數值有關。 Yoo (1999b) 在寬的圓柱間時展弦比 A=2 探討普郎特數影響的範 圍(0.3≦Pr≦1)會發生分歧現象與雙重解(dual solutions)。且當瑞里數超 過臨界瑞里數時，其範圍為 0.3≦Pr≦0.5 時研究發現在內圓柱會有二 次流現象的流場。當 0.3≦Pr≦0.4 時研究發現流體會產生磁滯現象 (hysteresis phenomenon)，但是範圍在 0.5≦Pr≦1 則會產生分歧現象。 Yoo (1999c) 研究更寬的圓柱間時展弦比 A=0.5 探討普郎特數影 響：在普郎特數Pr≦0.2 時，多渦捲流型式會在高瑞里數時觀測到；普 4

郎特數 0.3≦Pr≦1，瑞里數大於臨界瑞里數時，流場穩定後單渦捲流 與雙渦捲流都會存在，而雙渦捲流的總紐賽數會小於單渦捲流的紐 賽數；當普郎特數增高，所得到的臨界瑞里數也會增高。

Mizushima (2001) 研究水平圓柱間自然對流的轉換並假設在 2-D 與 不 可 壓 縮 流 場 ， 當 瑞 里 數 大 與 臨 界 瑞 里 數 時 利 用 數 值 模 擬 出 dual stable steady solutions。

現今，水平圓柱間自然對流流體在低普郎特數(Pr≦0.3)時，暫 態對流和邊界條件為固定的熱通量下會產生多變化的現象，但大部 份集中探討有關在圓柱表面的熱轉換並討論不同半徑比的水平圓柱 間自然對流的問題。在穩定的熱傳導範圍下，來探討普郎特數的影響， 和流體型態從穩定轉換到不穩定時，流場會產生非常複雜的多渦捲 型態與分歧現象這些結果都證明與普郎特數有關。

1.2 研究目的

本研究重點主要探討兩水平同心圓柱間受內外溫度差異所產生自 然對流流場，在不同參數的改變下，觀察穩定流場之行為。主要可分 為下列數項： 1. 藉由數值模擬兩水平同心圓柱間自然對流之流場，在低普郎特數 (Pr≦1.0)的範圍之內與兩水平同心圓柱間不同半徑比，來探討流場 在自然對流下轉換之過程與分歧現象(Bifurcation phenomenon) 的行為。 5

2. 當自然對流流場產生二次流時轉換到分歧現象與二次解的發生，探 討半徑比與瑞里數值所影響流場行為的範圍。

第 二 章

數 學 模 式

本文所探討之物理模式如圖一所示，兩水平無限長的同心圓柱 ，其間充滿黏滯性流體，T1 和 T2 分別為內外圓柱之溫度，R1 和 R2 分別為內外圓柱之半徑，半徑比(radius ratio)定義為η=R1/R2。 本文主要探討兩水平同心圓柱在自然對流作用下，改變各參數對 自然對流流場所產生之影響。在本章中將先建立其數學模式，以作為 數值方法解題之基礎。

2.1 統御方程式

在開始分析問題前，吾人先作以下的假設： A. 兩同心圓柱為無限長，流場為 2-D 流場， B. 流體為牛頓流體， C. 流體的物理性質皆為定值， D. 忽略流體的黏滯性耗損， E. 根據波斯尼克近似法(Boussinesq approximation)，除了浮力項 ，忽略流體密度隨溫度的變化量。 根據上述假設，其圓柱座標 (r,θ) 的統御方程式如下： 連續方程式 (continuity equation) :

0

v

r

1

)

ru

(

r

r

1

₌

θ

∂

∂

+

∂

∂

( 2. 1 ) 7

動量方程式 (momentum equations) : θ − β +     θ ∂ ∂ − − ∇ ν + ∂ ∂ ρ − = − θ ∂ ∂ + ∂ ∂ + ∂ ∂ cos ) T T ( g v r 2 u ) r 1 ( r p 1 r v u r v r u u t u 2 2 2 2 2 ( 2. 2 a )

θ

−

β

−









θ

∂

∂

+

−

∇

ν

+

θ

∂

∂

ρ

−

=

+

θ

∂

∂

+

∂

∂

+

∂

∂

sin

)

T

T

(

g

u

r

2

v

)

r

1

(

p

r

1

r

uv

v

r

v

r

v

u

t

v

2 2 2 2 ( 2. 2 b ) 能量方程式 (energy equation) : T T r v r T u t T _{= α∇}2 θ ∂ ∂ + ∂ ∂ + ∂ ∂ ( 2. 3 ) 上式中 ₂ 2 2 2 2 2 θ r 1 r r 1 r ∂ ∂ + ∂ ∂ + ∂ ∂ = ∇ _， g 為重力加速度 (gravitational acceleration)， α為熱擴散係數 (thermal diffusivity)，

β為熱膨脹係數 (thermal expansion coefficient)， ρ為流體密度 (density)， ν為運動黏滯係數 (kinematic viscosity)。 邊界條件 (boundary conditions) : r = R1 ： u = v = 0 ， T = T1 ( 2. 4 a ) r = R2 ： u = v = 0 ， T = T2 ( 2. 4 b ) 8

2.2 無因次化

定義以下八個無因次參數可將上一節中四個統御方程式無因次 化： d r X = _， t d2 α = τ _， 2 1 2 T T T T − − = Θ _， α = ud U _， α = vd V _， 2 2 pd P ρα = _， α ν = Pr _， αν − β = g (T1 T2)d3 Ra 其中

d=R

2

－R

1

，

Pr 為普郎特數 (Prandtl number)， Ra 為瑞里數 (Rayleigh number)。 無因次化後之統御方程式表示如下： 連續方程式 (continuity equation) : 0 V X 1 ) XU ( X X 1 ₌ θ ∂ ∂ + ∂ ∂ ( 2. 5 ) 動量方程式 (momentum equations) : cosθ a R Pr Θ θ V X 2 )U X 1 (∆ Pr X P X V θ U X V X U U τ U 2 2 2 ⋅ ⋅ ⋅ +     ∂ ∂ − − + ∂ ∂ − = − ∂ ∂ + ∂ ∂ + ∂ ∂ ( 2. 6 a ) θ ⋅ ⋅ ⋅ Θ −     θ ∂ ∂ + − ∆ + θ ∂ ∂ − = + θ ∂ ∂ + ∂ ∂ + τ ∂ ∂ sin Ra Pr U X 2 V ) X 1 ( Pr P X 1 X UV V X V X V U V 2 2 ( 2. 6 b ) 9

能量方程式 (energy equation) : ∆Θ = θ ∂ Θ ∂ + ∂ Θ ∂ + τ ∂ Θ ∂ X V X U _{( 2. 7 )} 上式中 2 2 2 2 2 X 1 X X 1 X ∂θ ∂ + ∂ ∂ + ∂ ∂ = ∆ _。 邊界條件 (boundary conditions) : η − η = 1 X _：

_U

=

_V

=

₀

_，

Θ

=

1

_{( 2. 8 a )} η − = 1 1 X _：

_U

=

_V

=

₀

_，

Θ

=

0

_{( 2. 8 b )} 10

第 三 章

數 值 方 法

3.1 求取數值解的方法

在前一章方程式( 2. 5 )~( 2. 7 )中除了連續方程式為線性，其餘 三式皆為非線性，無法直接求得其解析解，所以使用數值方法來求 解此類問題。

3.1.1 雙重級數展開

在數值求解部份，由於流場在θ方向以週期性出現，因此以傅 立葉級數(Fourier series)展開，在徑向則選用具有高精確度及收斂快 速特性的契比希夫多項式(Chebyshev polynomials)，並以雙重級數展 開的形式，來模擬流場中的速度、壓力與溫度，而展開之係數則為 時間的函數。由於U 方向速度、壓力與溫度在水平圓柱間正視圖沿徑 向垂直線之對稱性，可將其表示為：

θ

ξ

φ

τ

=

∑

∑

+ = − =

m

cos

)

(

)

(

A

U

_{mn n} 1 N 2 n 1 M 0 m ( 3. 1 a )

θ

ξ

φ

τ

=

∑

∑

+ = =

m

sin

)

(

)

(

B

V

_{mn n} 1 N 2 n M 1 m ( 3. 1 b )

θ

ξ

τ

=

∑

∑

− = − =

m

cos

)

(

T

)

(

C

P

_{mn n} 1 N 0 n 1 M 0 m ( 3. 1 c )

2

1

m

cos

)

(

)

(

D

_{mn n} 1 N 2 n 1 M 0 m

ξ

−

+

θ

ξ

φ

τ

=

Θ

∑

∑

+ = − = ( 3. 1 d ) 11

上式中φn是滿足速度與溫度邊界條件的基底函數，表示如下：

( )

[

]

[

( )

]

2

T

1

1

2

T

1

1

T

n 1 n 0 n n

=

−

−

−

−

+

−

φ

_{， n= 2 , 3 , 4 , ……} Tn(ξ)為第一型 n 階契比希夫多項式(Chebyshev polynomials)， 其定義如下：

)

cos

n

cos(

)

(

T

1 n

ξ

=

⋅

−

ξ

， n= 0 , 1 , 2 , …… 其中ξ 之定義域為[-1,1]。為了方便處理，可將原方程式中Ｘ之定 義域經由 η − η + − = ξ 1 1 X 2 _{的關係式，由} η − ≤ ≤ η − η 1 1 X 1 轉換至

−

1

≤

ξ

≤

1

。

3.1.2 離散方程式

將方程式( 2. 5 )~( 2. 7 )式，加以離散化(discretization)，建立離散 方程式。為了大幅減少求得收斂解所需的計算時間，在連續方程式均 以隱式法(Implicit Method)來表示，動量方程式中，時間微分項 以ξ= − ₋+_ηη 1 1 X 2 _{表示。線性項部份，除壓力項採用Implicit Method 外，} 其它均使用Crank-Nicolson Method 處理，即 2 f f f i 1 i ₊ = + 。至於非線性 的部份，則採用Adam-Bashforth Method 處理，即 2 f f 3 f i 1 i ₋ = + 。則未 知數之最新一次的值，可由前兩次的結果求得。 12

經推導出的離散方程式，如下所示： 連續方程式(continuity equation)

0

V

U

)

X

X

1

(

i 1 i 1

₌

θ

∂

∂

+

∂

∂

+

+ + ( 3. 2 ) 動量方程式(momentum equations) θ ⋅ ⋅ ⋅ τ ⋅ Θ + − τ + − τ − θ ∂ ∂ ⋅ ⋅ τ + θ ∂ ∂ ⋅ ⋅ τ −     ₊ τ _{⋅ ⋅ ∆ −} = ∂ ∂ ⋅ τ +     ₋ τ _{⋅ ⋅ ∆ −} + − − − − + + cos Ra Pr d ) V V 3 ( X 2 d ) U R U R 3 ( 2 d V Pr X d V Pr X d 3 U ) X 1 ( Pr 2 d 1 X P d U ) X 1 ( Pr 2 d 1 1 i 1 i 2 i 2 1 i 1 i i i 1 i 2 i 2 i 2 1 i 1 i 2 ( 3. 3 a ) θ ⋅ ⋅ ⋅ τ ⋅ Θ − − τ − − τ − θ ∂ ∂ ⋅ ⋅ τ − θ ∂ ∂ ⋅ ⋅ τ +     ₊ τ _{⋅ ⋅ ∆ −} = θ ∂ ∂ ⋅ τ +     ₋ τ _{⋅ ⋅ ∆ −} + − − − − − + + sin Ra Pr d ) V U V U 3 ( X 2 d ) V R V R 3 ( 2 d U Pr X d U Pr X d 3 V ) X 1 ( Pr 2 d 1 P X d V ) X 1 ( Pr 2 d 1 1 i 1 i 1 i i i 1 i 1 i i i 1 i 2 i 2 i 2 1 i 1 i 2 ( 3. 3 b ) 能量方程式(energy equation) ) R R 3 ( 2 d ) 2 d 1 ( ) 2 d 1 ( ₋ τ_{∆ Θ}i+1 _{= +} τ_{∆ Θ}i ₋ τ i_Θi ₋ i−1_Θi−1 ( 3. 4 ) 上式中

=

_∂

∂

+

_∂

∂

_θ

X

V

X

U

R

i i i ，

=

_∂

∂

+

_∂

∂

_θ

− − −

X

V

X

U

R

1 i 1 i 1 i 。 13

3.1.3 代數方程式

將( 3. 1 )式模擬各速度分量、壓力與溫度之多項式代入上式的離 散方程式中，在徑向以定置法(Collocation methed)處理，而定置點取為 ) 1 N L cos( L ₊ π = ξ ，L= 1 , 2 , 3 , ….. N。在軸向則是以 Galerkin method處 理。如此可將統御方程式轉化為代數方程式，以矩陣表示如下： 1 i, i 1 i

_F

AX

+

=

− ( 3. 5 ) 其中，矩陣A 的大小為 ，矩陣F 的大小為 。 其組成如下： 2 ) N M 4 ( × × (4×M× N) 1 i, i 3 2 1 1 i mn mn mn mn 42 41 34 24 23 22 14 13 11

0

F

F

F

D

C

B

A

0

0

A

A

A

0

0

0

A

A

A

0

A

A

0

A

+ −

























=

















































其中，矩陣A 內各矩陣及矩陣 F 內各矩陣表示式如下，

)

(

)

X

1

X

m

D

X

1

D

(

Pr

2

d

1

2

1

A

_{2 2 n} 2 2 mn 11



⋅

φ

ξ











−

−

+

⋅

⋅

τ

−

=

)

(

DT

2

d

A

mn _n 13

⋅

ξ

τ

=

)

1

r

(

Ra

Pr

d

4

1

d

r

cos

]

)

m

1

cos(

)

m

1

[cos(

)

(

D

Ra

Pr

4

d

D

A

_n 2 n 1 i mn 0 m 1 i mn mn 14

=

⋅

⋅

τ

⋅

ξ

−

−

θ

θ

θ

+

+

θ

−

⋅

ξ

φ

⋅

⋅

⋅

⋅

π

τ

−

=

⋅

∫

∑

∑

π π − = + = +

)

(

)

X

1

X

m

D

X

1

D

(

Pr

2

d

1

2

1

A

_{2 2 n} 2 2 mn 22



⋅

φ

ξ











−

−

+

⋅

⋅

τ

−

=

14

)

(

T

m

X

2

d

A

mn _n 23

⋅

⋅

ξ

τ

−

=

)

1

r

(

Ra

Pr

d

4

1

d

r

sin

]

)

m

1

sin(

)

m

1

sin(

[

)

(

D

Ra

Pr

4

d

D

A

_n 2 n 1 i mn 0 m 1 i mn mn 24

=

⋅

⋅

τ

⋅

ξ

−

+

θ

θ

θ

+

+

θ

−

⋅

ξ

φ

⋅

⋅

⋅

⋅

π

τ

=

⋅

∫

∑

∑

π π − = + = +

)

(

)

X

m

D

X

1

D

(

2

d

1

2

1

A

_{2 n} 2 2 mn 34



⋅

φ

ξ











−

+

τ

−

=

)

(

)

XD

1

(

2

1

A

mn _n 41

=

+

⋅

φ

ξ

)

(

m

2

1

A

mn _n 42

=

⋅

⋅

φ

ξ

)

r

cos

,

V

V

3

(

X

2

d

)

r

cos

,

U

R

U

R

3

(

2

d

)

(

B

m

Pr

X

2

d

)

(

B

m

Pr

X

2

d

3

)

(

A

)

X

1

X

m

D

X

1

D

(

Pr

2

d

1

2

1

F

1 i 2 i 2 1 i 1 i i i n 2 n 1 i mn 1 m 2 n 2 n i mn 1 m 2 n 2 n i mn 0 m 2 2 2 2 mn 1

θ

−

τ

+

θ

−

τ

−

ξ

φ

⋅

⋅

⋅

⋅

τ

+

ξ

φ

⋅

⋅

⋅

⋅

τ

−

ξ

φ

⋅

⋅













−

−

+

⋅

⋅

τ

+

=

− − − = − = = = = =

∑

∑

∑

∑

∑

∑

)

r

sin

,

V

U

V

U

3

(

X

2

d

)

r

sin

,

V

R

V

R

3

(

2

d

)

(

A

m

Pr

X

2

d

)

(

A

m

Pr

X

2

d

3

)

(

B

)

X

1

X

m

D

X

1

D

(

Pr

2

d

1

2

1

F

1 i 1 i i i 1 i 1 i i i n 2 n 1 i mn 0 m 2 n 2 n i mn 0 m 2 n 2 n i mn 1 m 2 2 2 2 mn 2

θ

−

τ

−

θ

−

τ

−

ξ

φ

⋅

⋅

⋅

⋅

τ

+

ξ

φ

⋅

⋅

⋅

⋅

τ

−

ξ

φ

⋅

⋅













−

−

+

⋅

⋅

τ

+

=

− − − − = − = = = = =

∑

∑

∑

∑

∑

∑

15

)

r

cos

,

R

R

3

(

2

d

)

0

r

(

d

)

(

D

)

X

m

D

X

1

D

(

2

d

1

2

1

F

1 i 1 i i i n 2 n i mn 0 m 2 2 2 mn 3

θ

Θ

−

Θ

τ

−

=

Χ

τ

−

ξ

φ

⋅

⋅













−

+

τ

+

=

− − = =

∑

∑

上式中 n n n n n n

d

d

2

dX

d

D

ξ

=

=

_，

θ

θ

=

ξ

φ

sin

n

sin

n

2

)

(

D

_{n ( n= 2 , 4 , 6 , …. ) ，}

2

sin

n

sin

n

2

)

(

D

_n

−

θ

θ

=

ξ

φ

_{( n= 3 , 5 , 7 , …. ) ，}

)

sin

cos

n

sin

n

sin

n

cos

n

(

4

)

(

D

₃ 2 n 2

θ

θ

θ

+

θ

θ

−

=

ξ

φ

( n= 2 , 3 , 4 , …. ) ，其中

θ

=

cos

1

ξ

L ， − ) 1 N L cos( L ₊ π = ξ ( L = 1 , 2 , 3 , ….. N ) ， θ ⋅ π =

∫

π π − f gd 2 1 ) g , f ( ， θ ∂ ∂ + ∂ ∂ = i i i _V X U R ， θ ∂ ∂ + ∂ ∂ = − − −1 i 1 i 1 i _V X U R 。 F1、F２及F３分別為徑向、θ方向與溫度，前一次線性項及前兩次 非線性項的和。

利用高斯消去法(Gauss Elimination Method)解矩陣方程式 ( 3. 5 )式，以( 3. 6 )式表示之，先解出係數溫度 Dmn，再逐步地將係

數壓力 Cmn 、速度 Amn、速度 Bmn 解出，直到滿足收斂條件在相同

無因次時間(time step)，任一係數前後兩次的差值與前一次的絕對比 值小於 10-4 視為收斂。為了使程式收斂結果達到收斂條件，初始條 件經過測試為 10-5~10-7，初始值先以 10-5 代入程式中，隨著瑞里數 增加，初始值取愈小。 4 i i 1 i 10 − + < φ φ − φ 將達到收斂條件的係數( 3. 6 )式代回( 3. 1 )式之中，可求出在不同 普郎特數、瑞里數及半徑比下自然對流流場的速度、壓力與溫度。 以高斯消去法解矩陣方程式( 3. 5 )式的表示式如下：

{

i 1

}

mn 13 1 1 11 1 i mn

A

F

A

C

A

+

=

−

−

+ ( 3. 6 a )

{

i 1

}

mn 23 2 1 22 1 i mn

A

F

A

C

B

+

=

−

−

+ ( 3. 6 b )

{

}

1 23 1 22 42 13 1 11 41 43 1 i mn

A

A

A

A

A

A

A

C

+

=

−

−

−

− −

{

F

₄

−

A

₄₁

A

₁₁−1

F

₁

−

A

₄₂

A

₂₂−1

F

₂

}

( 3. 6 c ) 3 1 34 1 i mn

A

F

D

+

=

− ( 3. 6 d ) 因為自然對流流場中的壓力只與微分項有關，修正壓力項的 結果，所以( 3. 6 c )式中多了 A43與F4兩個小矩陣，才能求解出流場 中正確的壓力值，因此得到壓力值之後，接著解流場速度分量及最 後觀察整個自然對流流場的行為時，才不會產生錯誤的結果！ 17

3.2 數值解的分析方法

3.2.1 流線圖

由解出的速度( 3. 1 a )式或( 3. 1 b )式，根據流線函數的定義取 任一式代入( 3. 7 a )式或( 3. 7 b )式，經積分後轉換成流線函數(stream function)，得到( 3. 8 )式，進而畫出在圓柱間流場中的流線(streamline)， 即表示流場在每一瞬間，任一點的切向量平行於該點之速度向量， 並利用此圖來觀察流場的行為及判斷其結果與流場物理現象中的對 稱性是否互相符合。 流線函數的定義： θ ∂ ψ ∂ = X 1 U _{( 3. 7 a )} X V ∂ ψ ∂ − = ( 3. 7 b ) 兩同心圓柱間之流線函數： t tan cons m sin m 1 A X N 1 _{mn n} 2 n 1 M 0 m + θ ⋅ φ ⋅ = ψ

∑

∑

+ = − =

( 3. 8 )

3.2.2 計算平均紐賽數( Nu)

由於兩同心圓柱間流場同時具有熱傳導與熱對流的效應，吾人可 計算出內圓柱的紐賽數及外圓柱的紐賽數來分析收斂結果是否正確， 當結果收斂時，內外圓柱分別計算出的紐賽數誤差值應在 5%以下才算 準確，而內外圓柱的平均紐賽數若大於1.0，則顯示出此同心圓柱間流 18

場的熱對流效應較強烈。 紐賽數的定義： k X h Nu = ⋅ 流體運動純傳導無因次熱傳率： 1 2 conduction x x ln 1 Nu = 內外圓柱局部紐賽數：

1 1 X X 1 2 X ln _X X _X Nu = ∂ ⋅ − = X ∂Θ 2 X X X ln Nu = − ⋅ ∂Θ 2 2 X X 1 X X ∂ ₌ 內外圓柱平均紐賽數： θ Θ ∂ ⋅ − = 1

∫

π ln X X d Nu 2 2 1 ∂ π X X ₌ 2 1 X X 1 0 θ Θ ∂ ⋅ − = 1

∫

π ln X X d Nu 2 2 2 ∂ π X X ₌ 2 2 X X 1 0 19

推導之後的內外圓柱平均紐賽數：

( )

    − − ⋅ − ⋅ − = X ln X

∑

∑

2n 1 D

∑

∑

2D 1 Nu n 2 2 1   = = −1 = = X_{1 m 0 n 2} mn _{m 0 n 2} mn 1       _{⋅ − −} ⋅ − =

∑

∑

∑

∑

= = = = 1 D 2 D n 2 X X ln X Nu _mn 2 n 0 m mn 2 2 n 0 m 2 2 2 1 上式中中括號內第一項、第二項只取當m=0 時，第二項 n 只取奇數時 。 20

第 四 章

結 果 與 討 論

為了瞭解程式所得數據之精確性，在本章中將先對程式做準確度 分析，作為往後計算選擇級數展開項數的依據。再探討圓柱間自然對 流之分歧現象的穩定臨界值（Raｃ），分析普郎特數與半徑比二個參數 的改變對流場行為的影響。

4.1 準確度分析

對兩同心圓柱間自然對流流場的計算，前一章數值方法中軸向與 徑向的 Fouries series 及 Chebyshev polynomials 雙重級數展開項 數（Ｍ、N）將影響數值計算結果的準確度，因此在進行自然對流 流場計算之前，先測試計算時所需的項數，以作為往後計算時的依據。 理論上，展開項數取的愈多其結果之準確度愈高，但項數取太多將 耗費大量的時間在運算上，因此取最少的展開項數而能達到相當準 確度才是最佳的選擇。 考慮半徑比與普郎特數兩參數對兩同心圓柱間平均紐賽數的變化 作準確度測試，圖二及圖三為測試的結果並與 Yoo（1999）比較，當 M 及 N 均以 5 項展開所得到的平均紐賽數誤差為 3.44%，但 M 及 N 均取至 10 項以上所得到的平均紐賽數已減至 0.05%，收斂結果和 Yoo 以有限差分法解的結果十分接近。同樣的，圖四及圖五顯示在 低瑞里數時較少的展開項數(10×10)已有不錯的結果，但在較高的瑞 里數時則需要較多的展開項數(12×12)。 圖六測試穩定臨界值與展開項數之關係，結果也發現當M 及 N 均 21

以6、7、8、9、10 項展開時所得到的臨界瑞里數有些許誤差，一直到 M 及 N 均取至 11 項以上所得到臨界瑞里數誤差才縮減至 1.05%以內。 收斂結果和Yoo（1996）以有限差分法解的結果十分接近。 由以上的分析可知，在低瑞里數時M 及 N 各取 11 項展開準確度 已足夠，因此在本文中探討分歧現象之臨界瑞里數時M 及 N 各取 11 項計算，而對瑞里數較大時則增加展開項數以獲得準確的結果。

4.2 臨界瑞里數對圓柱間自然對流流場之影響

對於水平同心圓柱間自然對流流場產生很多不同流場型態與轉 換現象，均取決於普郎特數與半徑比的變化。 本文中以兩種不同半徑比，分別改變普郎特數來觀察穩定流場 由單一渦漩流（crescent-shaped eddy）轉變為在內圓柱上方分離成兩 個相反方向旋轉的渦漩流（two counter-rotating eddies），得到的結果有

：當圓柱半徑為 0.2，隨著增加普郎特數，液動不穩定性會增大，流 場會發生分歧現象的臨界瑞里數愈大；當圓柱半徑比為 0.5，在 0.3 ≦Pr≦0.5，自然對流流場中內熱圓柱會產生二次流的現象，隨著增 加普郎特數，由於圓柱間寬距變小，液動不穩定性會減小，流場會 發生分歧現象的臨界瑞里數愈小。分別就普郎特數及半徑比作以下 的探討。

4.3 探討普郎特數對自然對流流場之影響

4.3.1 半徑比η=0.2 時，普郎特數對自然對流流場之影響

22

考慮圓柱半徑比小於 0.5，圖七為η=0.2，Pr=0.3 時，流場由單 一渦漩流轉變至雙渦漩流情形：當 Ra=1000 時，流場單渦漩流受到 的浮力作用很小，流動方向為靠內圓柱向上流動，靠外圓柱向下流動 ；當 Ra=2000 時，流場單渦漩流形狀類似腎臟，逐漸向下沈；當 Ra= 2800 時，流場內熱圓柱上方發生逐漸向下沈的單渦漩流分離成兩個 逆時針旋轉的雙渦漩流，繼續增加瑞里數時，內熱圓柱上方分離出 的第二個渦漩流旋轉更劇烈。 隨著增加普郎特數，圖八為η=0.2，Pr=0.5 時，流場由單一渦 漩流轉變至雙渦漩流情形：當 Ra=1000 時，流場單渦漩流受到的浮力 作用比Pr=0.3 時來的大，流動方向為靠內圓柱向上流動，靠外圓柱向 下流動；當Ra=6000 時，流場單渦漩流形狀類似腎臟，才逐漸向下沈； 當 Ra=6900 時，流場內熱圓柱上方發生逐漸向下沈的單渦漩流分離成 兩個逆時針旋轉的雙渦漩流，繼續增加瑞里數時，內熱圓柱上方分離 出的第二個渦漩流旋轉會更加劇烈。 隨著增加普郎特數，圖九為η=0.2，Pr=0.7 時，流場由單一渦 漩流轉變至雙渦漩流情形：當 Ra=1000~9000 時，流場單渦漩流受 到的浮力作用比Pr=0.5 時來的大，逐漸向上浮，流動方向為靠內圓柱 向上流動，靠外圓柱向下流動；當 Ra=9500 時，流場內熱圓柱上方發 生逐漸向下沈的單渦漩流分離成兩個逆時針旋轉的雙渦漩流，繼續 增加瑞里數時，內熱圓柱上方分離出的第二個渦漩流旋轉會更加劇烈。 綜合以上觀察所知，當 Pr＝0.3 時，流場中分子動量的擴散率遠 小於分子熱量的擴散率，因此在流線圖上，由於單渦漩流受到向上 浮的作用沒有這麼強烈，所以液動不穩定性很小，流場單渦漩流分 23

離成雙渦漩流之臨界瑞里數較小。在等溫線圖上顯示出，當瑞里數 較小時，只有靠近內熱圓柱下方溫度較高的地方比較集中，內熱圓 柱上方溫度分佈較均勻；當瑞里數增加，除了內熱圓柱下方溫度呈 現較密集的分佈情形，內熱圓柱上方因受到分離出的另一個渦漩流 的影響，溫度呈現少許不穩定性的分佈情形，其餘部份就變得較均勻。 當 Pr=0.7 時，流場中分子動量擴散率略低於分子熱量擴散率， 因此在流線圖上，由於單渦漩流受到浮力作用較大，所以液動不穩 定性很大，流場單渦漩流分離成雙渦漩流之臨界瑞里數比 Pr=0.3 更大。在等溫線圖上顯示出當瑞里數較小時，只有靠近內熱圓柱下 方溫度較高的地方比較集中，內熱圓柱上方溫度受浮力影響，溫度 分佈較往上疏離；當瑞里數增大，除了內熱圓柱下方溫度呈現更密 集的分佈情形，內熱圓柱上方因受到分離出的另一個渦漩流的影響 ，溫度呈現極大不穩定性的分佈情形，其餘部份也變得較不穩定。

4.3.2 半徑比η=0.5 時，普郎特數對自然對流流場之影響

考慮圓柱半徑比為0.5，圖十為η=0.5，Pr=0.3 時，流場由單一渦 漩流轉變至雙渦漩流情形：當 Ra=1600 時，流場單渦漩流中上方產生 了小渦漩，這時還未分離成雙渦漩流，流動方向為靠內圓柱向上流動， 靠外圓柱向下流動；當 Ra=2000 時，流場內熱圓柱上方出現了一個小 渦漩，逐漸分離成兩個逆時針旋轉的雙渦漩流；當 Ra=2900 時，流場 雙渦漩流下方的渦漩流中，短暫產生了向下流的小渦漩，這與 Yoo （1999）的研究在 0.3≦Pr≦0.5，自然對流流場中內熱圓柱會產生二 次流的現象十分吻合；當 Ra=3000 時，流場內多渦漩流會組合成為 單渦漩流；當 Ra=3100 時，流場內熱圓柱上方出現了一個小渦漩， 又逐漸分離成兩個逆時針旋轉的雙渦漩流，繼續增加瑞里數時，內 24

熱 圓 柱 上 方 分 離 出 的 第 二 個 渦 漩 流 旋 轉 更 劇 烈 。 其 收 斂 結 果 與 Yoo（1999）以有限差分法解的結果十分接近。 隨著增加普郎特數，圖十一為η=0.5，Pr=0.5 時，流場由單一 渦漩流轉變至雙渦漩流情形：當 Ra=2900 時，流場單渦漩流中上方 產生了小渦漩，這時還未分離成雙渦漩流，流動方向為靠內圓柱向 上流動，靠外圓柱向下流動；當 Ra=3000 時，流場內熱圓柱上方出 現了一個小渦漩，逐漸分離成兩個逆時針旋轉的雙渦漩流；當 Ra= 4000 時，自然對流流場中內熱圓柱會產生二次流的現象，流場雙渦 漩流下方的渦漩流中，短暫產生了向下流的小渦漩，繼續增大瑞里 數時，內熱圓柱上方分離出的第二個渦漩流旋轉會更劇烈。 隨著增加普郎特數，圖十二為η=0.5，Pr=0.7 時，流場由單一 渦漩流轉變至雙渦漩流情形：當 Ra=2900 時，流場內熱圓柱上方出 現了一個小渦漩，逐漸分離成兩個逆時針旋轉的雙渦漩流，繼續增 加瑞里數時，內熱圓柱上方分離出的第二個渦漩流旋轉會更劇烈。 收斂結果和Yoo（1996）以有限差分法解的結果十分接近。 綜合以上觀察所知，當 Pr＝0.3 時，流場中分子動量的擴散率遠 小於分子熱量的擴散率，因此在流線圖上，由於單渦漩流受到向上 浮的作用沒有這麼強烈，但是受限於圓柱間寬距變小，所以液動不 穩定性很大，流場單渦漩流分離成雙渦漩流之臨界瑞里數較大。在等 溫線圖上顯示出，當瑞里數較小時，只有靠近內熱圓柱下方溫度較 高的地方比較集中，內熱圓柱上方溫度分佈較均勻；當瑞里數增加 ，除了內熱圓柱下方溫度呈現較密集的分佈情形，內熱圓柱上方因 受到分離出的另一個渦漩流的影響，溫度呈現少許不穩定性的分佈 情形，其餘部份就變得較均勻。 25

當 Pr=0.7 時，流場中分子動量擴散率略低於分子熱量擴散率， 因此在流線圖上，由於受限於圓柱間寬距變小，但是單渦漩流受到 浮力作用比 Pr=0.3 稍大，所以液動不穩定性減小，因此流場單渦漩 流分離成雙渦漩流之臨界瑞里數比 Pr=0.3 略為低些。在等溫線圖上 顯示出，當瑞里數較小時，只有靠近內熱圓柱下方溫度較高的地方 比較集中，內熱圓柱上方溫度分佈較均勻；當瑞里數增大，除了內 熱圓柱下方溫度呈現較密集的分佈情形，內熱圓柱上方因受到分離 出的另一個渦漩流的影響，溫度呈現極不穩定性的分佈情形，其餘 部份就變得較均勻。 圖十三、十四、十五為η=0.2 與η=0.5，普郎特數在三個情況 下所計算出圓柱間的平均紐賽數曲線圖，由圖上顯示出，當η=0.2 ，隨著普郎特數增加，臨界瑞里數會增大，所計算出圓柱間的平均 紐賽數也就愈大；當η=0.5，隨著普郎特數增加，臨界瑞里數會略 為減小，所計算出圓柱間的平均紐賽數也就略為減小。

4.4 探討圓柱半徑比對自然對流流場的影響

考慮圓柱半徑比為 0.2 時，由圖七、圖八、圖九三個不同普郎 特數作分析，發現普郎特數增加，液動不穩定性及溫度不穩定性皆 會增大，流場中單渦漩流轉變成雙渦漩流之臨界瑞里數愈大，所計 算出圓柱間的平均紐賽數也隨著瑞里數增加而變大。 另外，考慮圓柱半徑比為 0.5 時，由圖十、圖十一、圖十二三 個不同普郎特數作分析，發現普郎特數增加，液動不穩定性減小， 流場中單渦漩流轉變成雙渦漩流之臨界瑞里數略為減小，當 0.3≦ 26

Pr≦0.5，液動不穩定性很大，所以在內熱圓柱會產生二次流現象的 流場，但是溫度不穩定性會增大，所計算出圓柱間的平均紐賽數也隨 著瑞里數增加而略為減小。 由以上觀察得知，當固定普郎特數，圓柱半徑比為 0.5 以下時 ，如圖十六所示，隨著半徑比愈小，液動不穩定性會增大，流場由 單渦漩流轉變成雙渦漩流之臨界瑞里數也愈大。 當固定普郎特數，圓柱半徑比為 0.5 以上時，隨著半徑比愈大 ，圓柱間寬距變小，液動不穩定性減小，流體所受到的浮力作用愈小 ，流場由單渦漩流轉變成雙渦漩流的臨界瑞里數愈小。 27

第 五 章

結 論

本文以數值方法探討水平同心圓柱間自然對流流場之行為，以兩 種 不 同 半 徑 比 ， 分 別 改 變 普 郎 特 數 來 觀 察 穩 定 流 場 由 單 一 渦 漩 流（crescent-shaped eddy）轉變為在內熱圓柱上方分離成兩個相反方 向旋轉的渦漩流（two counter-rotating eddies），根據第四章的結果可 歸納以下幾點結論： （1） 當圓柱半徑為 0.2：隨著普郎特數增加，液動不穩定性及溫度 不穩定性皆會增大，流場發生分歧現象的臨界瑞里數愈大， 所計算出圓柱間平均紐賽數也隨著瑞里數增加而變大。 （2） 當圓柱半徑為 0.5：隨著普郎特數增加，由於圓柱間寬距變小 ，液動不穩定性會減小，流場發生分歧現象之臨界瑞里數愈小 ，但是溫度不穩定性會增大，所計算出圓柱間平均紐賽數也 隨著瑞里數增加而略為減小。另外，當 Pr=0.3 及 Pr=0.5 時， 內熱圓柱會有二次流現象的流場，這與 Yoo（1999）的研究 在 0.3≦Pr≦0.5，自然對流流場中內熱圓柱會產生二次流的現 象十分吻合；當 Pr=0.7，Ra=2900 時，流場內熱圓柱上方出 現了一個小渦漩，逐漸分離成雙渦漩流。其收斂結果與 Yoo （1999）、（1996）均以有限差分法解的結果十分接近。 （3） 圓柱半徑比為 0.5 以下時，隨著半徑比愈小，液動不穩定性 會增大，流場由單渦漩流轉變成雙渦漩流之臨界瑞里數也愈大 ；圓柱半徑比為 0.5 以上時，隨著半徑比愈大，圓柱間寬距 變小，液動不穩定性減小，流場由單渦漩流轉變成雙渦漩流 28

的臨界瑞里數愈小。 （4） 吾人在本文中所使用雙重級數展開的數值方法，經測試結果： 無法觀察到較高瑞里數時，圓柱間穩定流場之行為；當圓柱半 徑比為 0.8 以上時（狹窄寬距），無法觀察到圓柱間穩定流場 之分歧現象，因此，所使用雙重級數展開的數值方法，仍有 改進的空間。 29

參 考 文 獻

1. Beckmann, W. 1931 Die Warmeubertragung in zylindrischen Gasschichten bei naturlicher Konvektion. Forsch. Geb. d.

Ingenieurwesen, 2(5), 165-178.

2. Voigt H. and Krischer,O. 1932 Die Warmeubertragung in zylindrischen Luftschichten bei naturlicer Konvektion. Forsch. Geb. d.

Ingenieurwesen, 3(6), 303-306.

3. Crawford, L. and Lemlich, R. 1962 Natural convection in horizontal concentric annuli. I.E.C. Fund. 1, 260-264.

4. Abbott, M. R. 1964 A numerical method for solving the equations of natural convection in a narrow concentric cylindrical annulus with a horizontal axis. Q. Jl Mech. Appl. Math. 471-481.

5. Powe, R. E. Carley, C. T. and Bishop, E. H. 1969 Free Convection flow patterns in cylinders annuli. Journal of Heat Transfer 91, 310-314.

6. Powe, R. E. Carley, C. T. and Carruth, S. L. 1971 A numerical solution for natural convention in cylindrical annuli. Journal of Heat Transfer 93(12), 210-220.

7. Kuehn, T. H. Goldstein, R. J. 1976 An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders. Journal of Fluid Mechanics 74, 695-729.

8. Charrier-Mojtabi, M. C. Mojtabi, A. Caltagirone, J. P. 1979 Numerical solution of a flow due to natural convection in horizontal cylindrical annulus. Journal of Heat Transfer 101, 171-173.

9. Busse, F. H. 1981 Transition to turbulence in Rayleigh-benard convection, in: H. L. Swinney, J. P. Gollub(Eds.). Topics in

AppliedPhysics 45, Springer-Verlag, 97-137.

10. Tsui, Y. T. Tremblay, B. 1984 On transient natural convection heat transfer in the annulus between concentric horizontal cylinders withisothermal surfaces. International Journal of Heat and Mass

Transfer27, 103-111.

11. Rao, Y. S. Miki, Y. Fukuda, K. Takata, Y. Hasegawa, S. 1985 Flow patterns of natural convection in horizontal cylindrical annuli.

International Journal of Heat and Mass Transfer 28, 705-714.

12. Choi, J. Y. Kim, M. U. 1993 Three-dimensional liner stability of natural convective flow between concentric horizontal cylinders.

International Journal of Heat and Mass Transfer 36, 4173-4180.

13. Cadiou, P. Desrayaud, G. Lauriat, G. 1998 Natural convection in a narrow horizontal annulus：the effects of thermal and hydrodynamic instabilities. Journal of Heat Transfer, Transactions ASME 120, 1019-1025.

14. Yoo, J-S. 1998 Natural convection in a narrow horizontal cylindrical annulus: Pr≦0.3. International Journal of Heat and Mass Transfer 41, 3055-3074.

15. Yoo, J-S. 1999a. Transition and multiplicity of flows in natural convection in a narrow horizontal cylindrical annulus: Pr=0.4.

International Journal of Heat and Mass Transfer 42, 709-722.

16. Yoo, J-S. 1999b. Prandtl number effect on bifurcation and dual solutions in natural convection in a horizontal annulus. International

Journal of Heat and Mass Transfer 42, 3279-3290.

17. Yoo, J-S. 1999c. Prandtl number effect on transition of free -convective flows in a wide-gap horizontal annulus. International

Comm. Heat and Mass Transfer 26, 811-817.

18. Dyko, M. K. Vafai, K. Mojtabi, A. K. 1999 A numerical and

experimental investigation of stability of natural convective flows within a horizontal annulus. Journal of Fluid Mechanics 381, 27-61. 19. Mizushima, J. 2001 Transitions of natural convection in horizontal

annulus. International Journal of Heat and Mass Transfer 44, 1249-1257.

20. 韓瑞忠 1996 調制泰勒渦旋流的數值解. 國立交通大學碩士論文.

T1 > T2

θ

T2 T1 g R1 R2 圖一：基本物理模式示意圖 33

1.5 1.6 1.7 1.8 1.9 2 5×5 6×6 7×7 8×8 9×9 10×10 11×11 12×12 13×13 14×14 MxN Present Nu=1.68 Yoo(1999) Nu ) 1999 ( Yoo ~ 68 . 1 Nu= 圖二 η=0.2，Pr=0.5，Ra=5000 時，測試各展開項數所計算出的平均 紐賽數，並與Yoo（1999）比較。 34

Nu 1.3 1.4 1.5 1.6 1.7 1.8 7×7 8×8 9×9 10×10 11×11 12×12 13×13 14×14 MxN Present

Nu=1.52 Yoo(1999)1.52~Yoo(1999)

Nu=

圖三 η=0.2，Pr=0.3，Ra=4000 時，測試各展開項數所計算出的平均 紐賽數，並與Yoo（1999）比較。

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1000 1500 2000 2500 3000 3500 4000 Ra MxN=10x10 Yoo(1999) MxN=12x12 Nu 圖四 η=0.2，Pr=0.3 時，每個瑞里數所計算出的平均紐賽數，並與 Yoo（1999）比較。 36

Nu 1 1.2 1.4 1.6 1.8 2 2.2 1000 2000 3000 4000 5000 6000 7000 Ra MxN=10x10 Yoo (1999) MxN=12x12 圖五 η=0.2，Pr=1.0 時，每個瑞里數所計算出的平均紐賽數，並與 Yoo（1999）比較。 37

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 6×6 7×7 8×8 9×9 10×10 11×11 12×12 13×13 14×14 15×15 MxN Present

Ra=2850 Yoo(1996)2850~Yoo(1996)

Rac = c Ra 圖六 η=0.5，Pr=0.7 時，測試各展開項數所測的臨界瑞里數，並與 Yoo（1996）比較。 38

（a） （b） 圖七 當η=0.2，Pr=0.3 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=2000（c） Ra=2500（d）Ra=2700（e）Ra=2800（f）Ra=3600（g） Ra=4000 39

（c） （d） （e） 續 40

（f） （g） 續 41

（a） （b） 圖八 當η=0.2，Pr=0.5 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=3000（c） Ra=4000（d）Ra=6000（e）Ra=6800（f）Ra=6900（g） Ra=8000 42

（c） （d） （e） 續 43

（f） （g） 續 44

（a） （b） 圖九 當η=0.2，Pr=0.7 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=3000（c） Ra=4000（d）Ra=6000（e）Ra=8000（f）Ra=9000（g） Ra=9500 45

（c） （d） （e） 續 46

（f） （g） 續 47

（a） （b） 圖十 當η=0.5，Pr=0.3 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1500（b）Ra=1600（c） Ra=1900（d）Ra=2000（e）Ra=2900（f）Ra=3000（g） Ra=3100 48

（c） （d） （e） 續 49

（f） （g） 續 50

（a） （b） 圖十一 當η=0.5，Pr=0.5 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=1000（b）Ra=2000（c） Ra=2900（d）Ra=3000（e）Ra=3500（f）Ra=4000（g） Ra=6000 51

（c） （d） （e） 續 52

（f） （g） 續 53

（a） （b） 圖十二 當η=0.5，Pr=0.7 時，兩同心圓柱間自然對流流場之流線圖 （左）及其等溫線圖（右），（a）Ra=2000（b）Ra=2500（c） Ra=2800（d）Ra=2900（e）Ra=3500（f）Ra=4000（g） Ra=5000 54

（c） （d） （e） 續 55

（f） （g） 續 56

Nu 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 500 1000 1500 2000 2500 3000 3500 4000 Ra η=0.2 η=0.5 圖十三 當 Pr=0.3 時，不同瑞里數所計算出的平均紐賽數。 57

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Ra η=0.2 η=0.5 Nu 圖十四 當 Pr=0.5 時，不同瑞里數所計算出的平均紐賽數。 58

1 1.2 1.4 1.6 1.8 2 2.2 2.4 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Ra η=0.2 η=0.5 Nu 圖十五 當 Pr=0.7 時，不同瑞里數所計算出的平均紐賽數。 59

0 2000 4000 6000 8000 10000 0.3 0.4 0.5 0.6 0.7 Pr η=0.2 η=0.3 η=0.4 η=0.5 c Ra 圖十六 當η=0.2 至η=0.5，臨界瑞里數與 Pr 之關係。 60