垂直通道內振動對高溫面熱傳增益之研究

(1)

國立交通大學

機械工程學系

博士論文

垂直通道內振動對高溫面熱傳增益之研究

Enhancement of heat convection by an oscillating heat

surface in a vertical channel

研究生：黃建平

指導教授：傅武雄博士

(2)

垂直通道內振動對高溫面熱傳增益之研究

Enhancement of heat convection by an oscillating heat surface in a

vertical channel

研究生：黃建平 Student：Chien-Ping Huang

指導教授：傅武雄 Advisor：Wu-Shung Fu

國立交通大學

機械工程學系

博士論文

A Dissertation

Submitted to Department of Mechanical Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Doctor of Philosophy

in

Mechanical Engineering

June 2007

Hsinchu, Taiwan, Republic of China

(3)

垂直通道內振動對高溫面熱傳增益之研究研究生：黃建平指導教授：傅武雄國立交通大學機械工程學系摘要本研究利用數值方法分析垂直管道內高溫面加入水平振動的因子後，對高溫面散熱效率的影響。又採用的數值模擬為Arbitrary Lagrangian-Eulerian（ALE）座標描述系統，配合葛拉金有限元素法（FEM），處理此種流體與固體介面間相互影響的問題。首先探討在垂直管道內振動高溫面於自然對流領域內提升高溫面熱傳效率的機制，並且針對不同的振幅與雷利數求得加入振動條件後，高溫面散熱效果可以優於穩態非振動高溫面散熱效果的臨界振動頻率；進而探討當垂直管道流進入混合對流與強制對流範圍時，在不同的振幅與葛瑞秀夫數下，亦可發現加入振動條件後，振動高溫面之散熱效果可以優於穩態非振動高溫面散熱效果的臨界振動頻率。雖然一般的觀念大都認為存在振動的熱表面其熱傳的效果會比沒有振動的情形下要好，但是由研究結果得知，在不同的振幅與雷利數或葛瑞秀夫數搭配時，加入振動後的熱傳現象仍可能比沒有振動時的熱傳效率差，也因此可以定義且找出臨界振動頻率以提供研究振動熱傳者之參考。最後，有關臨界振動頻率的經驗公式在強制對流、混合對流以至於自然對流的形式皆不相同，本文皆已分別推導出來，且在三種的對流形式中都有相當平滑的曲線可以獲得。

中文摘要

(4)

Enhancement of heat convection by an oscillating heat surface in a vertical

channel

Student：Chien-Ping Huang Advisor：Wu-Shung Fu

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

The study investigates the heat transfer enhancement by a oscillation motion of heat surface in a vertical channel flow. In the numerical simulation analysis, a Galerkin finite element method (FEM) with arbitrary Langrangian-Eulerian coordinate (ALE) is adopted to solve the flow and thermal fields for moving boundary problem. Firstly, the mechanism of heat transfer enhancement by oscillation heat surface is discussed in natural convection scene. For different oscillation amplitude and Rayleigh number, the frequency that the heat transfer rate of oscillation heat surface is higher than non-oscillation condition is found and called critical oscillation frequency. The limitation enhancement of similar physical model is still researched by this study in mixed and forced convection scenes for different amplitude and Grashof number.

Generally, that a heat surface subject to a oscillation motion is advantageous to convective heat transfer is a well-known tuition. However, based upon the results indicated above the heat transfer rates with certain combinations of the oscillation amplitude and Rayleigh number or Grashof number are possibly smaller than that of the stationary condition. Then the critical oscillation frequency could be regarded as a criterion of whether the heat transfer rate of the heat surface subject to a oscillation motion being larger than that of the stationary state or not . Finally, the so called critical oscillation frequency can be expressed as a function of oscillation amplitude and Grashof number in all the natural, mixed, and forced convection scenes.

(5)

誌謝

衷心的感謝指導教授傅武雄博士多年來在學業、工作與論文上的悉心指導，以及在為人處事方面的諄諄教誨，在此謹致最高的謝忱與敬意。同時也要感謝機械系諸師長在課業上的教導，使本研究得以順利完成。更要感謝父母含辛茹苦的養育之恩，內人張憶如與岳父母的鼓勵與支持，還有眾多長官、同事與朋友的關懷，今日方能順利完成學業。最後，謹將此喜悅與所有關心我的人共同分享。

(6)

中文摘要...i 英文摘要...ii 誌謝...iii 目錄...iv 表目錄...vi 圖目錄...vii 符號說明...x 第一章緒論 ...1 第二章ALE 法之運動學理論與數學模式 ...8 2.1 網格描述與座標系統 ...8 2.2 座標定義域 ...8 2.3 速度定義 ...9 2.4 座標轉換 ...10 2.5 ALE 統御方程式...12 2.5.1 強制對流於濃度擴散之移動邊界問題 ...12 2.5.2 強制對流於熱傳增益之移動邊界問題 ...13 2.5.3 自然與混合對流於熱傳增益之移動邊界問題 ...13 2.5.4 加入多孔性介質於熱傳增益之移動邊界問題 ...14 2.6 格點速度 ...15 第三章自然對流垂直管道流中振動對高溫壁面熱傳增益之研究 ...22 3.1 前言 ...22 3.2 物理模式 ...23 3.3 統御方程式與邊界條件 ...23

(7)

3.5.3 振動頻率及雷利數對熱傳之影響 ...28 3.5.4 振動幅度對熱傳之影響 ...29 3.5.5 臨界振動頻率與經驗公式 ...29 3.6 結論 ...30 第四章混合與強制對流垂直管道流中振動對高溫壁面熱傳增益之研究 ...54 4-1 前言...54 4-2 物理模式...55 4.3 統御方程式與邊界條件 ...55 4.4 數值方法 ...57 4.5 結果與討論 ...57 4.5.1 格點與時間步進測試 ...57 4.5.2 振動高溫面對流場與溫度場之影響 ...58 4.5.3 振動頻率對熱傳之影響 ...59 4.5.4 葛瑞秀夫數對熱傳之影響 ...60 4.5.5 臨界振動頻率 ...61 4.5.6 經驗公式 ...62 4.5.7 熱傳增益原因之探討 ...63 4.5.8 實際應用範圍 ...63 4.6 結論 ...64 第五章結論與建議 ...95 5.1 總結 ...95 5.2 建議 ...95 參考文獻...97 附錄一. 自然對流模式之有限元素法離散化推導程序 ...102 附錄二. 混合對流模式之有限元素法離散化推導程序 ...108 作者簡介... 112

(8)

表目錄

表 1-1 主要冷卻方法比較表 [1] ...5 表 1-2 各種熱傳增益技術文獻數目統計表 [8] ...6 表 3-1 設計參數組合表...31 表 3-2 穩態非振動高溫面平均紐塞數與雷利數之關係表...32 表 4-1 設計參數組合表...65 表 4-2 穩態非振動高溫面平均紐塞數與葛瑞秀夫數之關係表...66 表 4-3 不同葛瑞秀夫數之臨界振動頻率F 列表（_cc L_c =0.1）...67 表 4-4 不同葛瑞秀夫數之臨界振動頻率F 列表（_cc L_c =0.2） ...68 表 4-5 不同葛瑞秀夫數之臨界振動頻率F 列表（_cc L_c =0.4） ...69

(9)

圖目錄

圖 1-1 各種熱傳模式與冷卻流體之熱傳量與可達成溫度差關係圖 [5]...7 圖 2-1 ALE 座標系統關係圖 ...17 圖 2-2 網格變化圖 ...18 圖 2-3 網格變化研究之物理模式圖[30] ...19 圖 2-4 兩組網格示意圖[30]...20 圖 2-5 兩組網格速度之等相對誤差線與局部紐塞數分佈圖[30] ...21 圖 3-1 物理模式圖 ...33 圖 3-2 計算流程圖 ...34 圖 3-3 網格測試(Ra=103) ...35 圖 3-4 時間步進測試(Ra ₌104,L_c ₌0.1,F_c ₌500_{) ...36}

圖 3-5 本研究與 McAdams [45]、Churchill 和 Chu[46]研究結果比較圖 ...37

圖 3-6 第 11 振動週期之流線與等溫線變化(Ra₌104,L_c ₌0.1,F_c ₌500_)...38 圖 3-7 平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.1_{) ...41} 圖 3-8 時平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.1_)...42 圖 3-9 時平均紐塞數Nu隨時間變化圖(Ra₌104,L_c ₌0.1_)...43 圖 3-10 時平均紐塞數Nu隨時間變化圖(Ra₌105,L_c ₌0.1_)...44 圖 3-11 平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.2_{) ...45} 圖 3-12 時平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.2_)...46 圖 3-13 時平均紐塞數Nu隨時間變化圖(Ra₌104,L_c ₌0.2_)...47 圖 3-14 時平均紐塞數Nu隨時間變化圖(Ra₌105,L_c ₌0.2_)...48 圖 3-15 平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.4_{) ...49} 圖 3-16 時平均紐塞數Nu隨時間變化圖(Ra₌103,L_c ₌0.4_)...50 圖 3-17 時平均紐塞數Nu隨時間變化圖(Ra₌104,L_c ₌0.4_)...51 圖 3-18 時平均紐塞數Nu隨時間變化圖(Ra₌105,L_c ₌0.4_)...52

(10)

圖 3-19 臨界振動頻率F 隨雷利數_cc Ra變化關係圖...53 圖 4-1 物理模式圖...70 圖 4-2 網格測試(Gr =104,Re=102) ...71 圖 4-3 時間步進測試(Gr ₌104,Re₌102,L_c ₌0.1,F_c ₌2_{) ...72} 圖 4-4 垂直管道流中強制、混合與自然對流在層流與紊流場的分布情形[56]...73 圖 4-5 第 16 振動週期之流線與等溫線變化(Gr ₌105,Re₌102,L_c ₌0.1,F_c ₌10_{) ...74} 圖 4-6 平均紐塞數Nu隨時間變化圖(Gr₌100,Re₌102,L_c ₌0.1_{) ...75} 圖 4-7 時平均紐塞數Nu隨時間變化圖(Gr₌100,Re₌102,L_c ₌0.1_)...76 圖 4-8 平均紐塞數Nu隨時間變化圖(Gr₌101,Re₌102,L_c ₌0.1_{) ...77} 圖 4-9 時平均紐塞數Nu隨時間變化圖(Gr₌101,Re₌102,L_c ₌0.1_)...78 圖 4-10 平均紐塞數Nu隨時間變化圖(Gr₌102,Re₌102,L_c ₌0.1_{) ...79} 圖 4-11 時平均紐塞數Nu隨時間變化圖(Gr ₌102,Re₌102,L_c ₌0.1_)...80 圖 4-12 平均紐塞數Nu隨時間變化圖(Gr₌103,Re₌102,L_c ₌0.1_{) ...81} 圖 4-13 時平均紐塞數Nu隨時間變化圖(Gr₌103,Re₌102,L_c ₌0.1_)...82 圖 4-14 平均紐塞數Nu隨時間變化圖(Gr₌104,Re₌102,L_c ₌0.1_{) ...83} 圖 4-15 時平均紐塞數Nu隨時間變化圖(Gr₌104,Re₌102,L_c ₌0.1_)...84 圖 4-16 平均紐塞數Nu隨時間變化圖(Gr₌105,Re₌102,L_c ₌0.1_{) ...85} 圖 4-17 時平均紐塞數Nu隨時間變化圖(Gr₌105,Re₌102,L_c ₌0.1_)...86 圖 4-18 高溫面振幅L_c =0.1時，不同Gr下振動頻率對Nust Nu_st 變化圖(Re=102)...87 圖 4-19 高溫面振幅L_c =0.2時，不同Gr下振動頻率對Nust Nu_st 變化圖(Re=102) ...88 圖 4-20 高溫面振幅L_c =0.4時，不同Gr下振動頻率對Nust Nu_st 變化圖(Re=102) ...89

(11)

圖 4-24 F_c =20時，第 25 個振動週期後之等溫線圖(Gr₌104,L_c ₌0.2_{) ...93}

(12)

符號說明

c 物質濃度 [kg m-3] cˆ 對流速度 [m s-1] C 無因次物質濃度 ab D 濃度擴散係數 [m2 s-1] f 任意函數 c f 高溫面振動頻率 [s-1] c F 無因次高溫面振動頻率 1 , c F 第一無因次高溫面振動頻率（F <_c_,₁ F ） _c_c 2 , c F 第二無因次高溫面振動頻率（F >_c_,₂ F ） _c_c c c F 無因次高溫面臨界振動頻率 Gr 葛瑞秀夫數 (Grashof number) g 重力加速度 [m s-2] k 熱傳導係數 [W m-2 K-1] 0 l 垂直管道長度 [m] 1 l 高溫面上緣至垂直管道頂部長度 [m] 2 l 高溫面長度 [m] 3 l 高溫面下緣至垂直管道底部長度 [m] c l 高溫面振幅 [m] 2 L 無因次高溫面長度 3 L 無因次高溫面下緣至垂直管道底部長度 c L 無因次高溫面振幅 P N 振動一週期的時間步進次數

(13)

1 i Nu + 第 i+1 個週期之時平均紐塞數 st Nu 穩態平均紐塞數 st Nu _{穩定時平均紐塞數} st,1 Nu 高溫面振動頻率為F 時之穩定時平均紐塞數 _c_,₁ st,2 Nu 高溫面振動頻率為F 時之穩定時平均紐塞數 _c_,₂ Y Nu 局部紐塞數 p 壓力 [N m-2] ∞ p 參考壓力 [N m-2_] P 無因次壓力 Pr 普朗特數 (Prandtl number，Pr =υ α) a R 雷利數 (Rayleigh number) e R 雷諾數 (Reynolds number) Sc 修密特數 (Schmidt number) t 時間 [s] T 溫度 [K] c T 進口溫度 [K] h T 高溫面溫度 [K] v u, x和y 方向速度 [m s-1] c u 高溫面振動速度 [m s-1] m u 高溫面最大振動速度 [m s-1] uˆ 格點速度 [m s-1] V U, x和y 方向無因次速度 c U 無因次高溫面振動速度 m U 無因次高溫面最大振動速度 Uˆ 無因次格點速度 0 v 進口速度 [m s-1] w 垂直管道寬度 [m]

(14)

y x, 卡式座標 [m] i x 參考定義域中之某點座標 [m] i y 空間定義域中之某點座標 [m] i z 物質定義域中之某點座標 [m] Y X, 無因次卡式座標希臘符號 α 熱擴散係數 [m2_s-1_] β 體積膨脹係數 [K-1] θ 無因次溫度 λ 處罰函數 ν 動黏滯係數 [m2 s-1] ρ 密度 [kg m-3] τ 無因次時間 P τ 無因次週期時間 Φ 計算變數 t Φ 物質定義域與空間定義域之座標轉換函數 t ˆ Φ 參考定義域與空間定義域之座標轉換函數 t Ψ 物質定義域與參考定義域之座標轉換函數 x Ω 參考定義域 y Ω 空間定義域 z Ω 物質定義域其他

(15)

第一章緒論

高溫壁面或電子元件表面的散熱問題在工業界或電腦資訊設備研發中非常常見，如何有效增加高溫壁面與電子元件的熱傳效益，一直是業界與學術界致力探索的主題。冷卻高溫壁面的模式主要有空氣自然對流(Natural air convection)、空氣強制對流(Forced air convection)、空氣強制對流加水冷式熱交換器(Forced air convection plus water-cooled heat exchanger)以及包括蒸發冷凝的液體冷卻(Liquid-cooled)與熱管(Heat pipe)等方法[1]。上述各種模式的特性比較如表1-1[1]所示。除了上述方法之外，近年來也有許多新的散熱觀念不斷被提出，Iida 等[2]在液體中加入鋁金屬微粒，當液體發生沸騰現象時，鋁金屬微粒也被液化，可將熱傳效率提高十倍。Fusegi [3]利用數值分析方法來模擬存在肋狀突起管道中之振動對向流，當提高振動對黏性的相對強度時，可顯著的提高熱傳效率。Sitter 等[4]則利用實驗的方法來探討高強度音場對重力場與微重力場下池沸騰的影響，發現音場有助於提高熱傳效率。儘管如此，空氣對流具有結構簡單、冷卻劑（空氣）取得容易、成本較低等特性，因而被廣泛的應用在各種電子元件的散熱中。然而，近年來隨著半導體製程進入奈米等級的尺寸以及電腦高速運算的要求，不同功能的電路經由整合後直接製作在同一晶片上，也因此單一晶片內的電子元件數目逐漸增加，電子構裝也有顯著的進步並朝向輕薄短小的構裝趨勢邁進，這也造成熱量集中造成現今電子元件的發熱量已經達到每平方公分數十瓦的等級，相較於圖1-1 [5]之各種冷卻高溫壁面方法對壁面熱傳量與溫度差示意圖，此等級的熱傳量已經達到空氣強制對流的極限。 Bergles [6,7]詳細地討論與整理各種增加熱傳效率的方法，並且將常見的增加熱傳效率方法略分為兩大類：一為不須另外作功的被動式(passive)方法，比如在壁面加一層塗料、將壁面作表面粗糙處理、增加熱傳表面積、加裝混合流體或使流體產生旋轉流動的機構、與在流體中加入添加物等；另一類為須外加能量的主動式(active)方法，包括有以機構增加流體的擾動、使壁面振動或振動流場、增加電磁場、噴入流體、移除覆蓋熱傳面的冷凝水膜、或加裝噴嘴直接對高溫壁面噴出流體等。此外，Bergles 與 Webb [8]進一步將過去有關提高熱傳效率的文獻，依據上述分類整理製表，如表1-2 所示。比較上述各種增加熱傳量的方法，不論利用主動式或被動式的方法來增加熱傳效率都會面臨到極限的問題。主要原因在於流場熱傳面上形成的速度與溫度邊界層妨礙熱的

(16)

的附近，受壁面影響的區域。在此區域中，剪應力使流體的流速緩慢的變化至與壁面速度相同。至於溫度邊界層則是流體溫度受壁面熱傳的影響，其溫度逐漸變化至與壁面溫度相同的區域。在溫度邊界層中，流體的溫度變化和緩，溫度梯度較小。根據熱傳導的傅利葉定律，在溫度梯度小的情況下，熱傳量也較小。因此，溫度邊界層的存在將限制壁面所能傳出之熱傳量。若想要大幅提高高溫壁面的散熱效率，必須破壞掉高溫面上的邊界層，使壁面直接與低溫流體接觸，進而提高壁面與流體間的溫度梯度，達成增加壁面熱傳量的效果。又轉動設備的運轉而造成發熱裝置振動之現象極為常見，例如許多機械或電子設備含有馬達或風機等轉動元件，而當該些元件運作時，設備上原有會發熱的電子零件如IC 等即會受到該些轉動馬達或風機的影響而產生振動。至於該些 IC 的散熱現象有可能在自然對流、強制對流或是混合對流下進行。以電腦為例，雖然電腦機殼多有一至三顆的風機運轉，但由於電路板的使用頻繁，包括主機板、顯示卡、網路卡、數據卡、音效卡以及硬碟、光碟機等等，除了發熱量極大的 CPU 或部分顯示晶片會加裝散熱風扇以強制對流的模式加強熱傳效率外，常常造成機殼內部分 IC 並非處於理想的強制對流散熱模式，而較為類似自然對流或混合對流的熱傳現象，此時又在機殼或 CPU 等散熱風扇的運轉下，即會有固定振動的現象產生。因此，存在振動現象的發熱電子元件之散熱問題相當常見，本文的主要目的即在探討此一熱傳現象，而此流體與振動物體間相互影響的問題，被歸類於動態的移動邊界問題(moving boundary problem)。

現今描述流體運動問題的座標系統可分為兩種：拉格朗日(Lagrangian)座標系統與尤拉法(Eulerian)座標系統。利用拉格朗日座標系統描述流場時，觀測流場的焦點集中在各個流體質點上；在數值計算的過程中，計算網格將隨著流體質點的運動而移動。因此，拉格朗日座標系統可以準確的描述移動邊界與固體液體介面間相互影響的情況。但是流體質點的運動方向與速度大小皆不相同，當流體質點間速度變化較大時，網格在計算的過程中將過度扭曲變形而影響計算。反之，尤拉座標系統將觀測流場的焦點固定在流場的空間中，描述此固定範圍中的流體質點運動情況，因此計算網格是固定在流場的空間中，不隨流體質點運動，避免了前述的網格扭曲問題。也使得尤拉座標系統在處理移動

(17)

系統之混合形式用來計算包含可動液體邊界之二維液動問題，且稱此種混合形式為「Coupled Eulerian-Lagrangian Code」，其證明此方法可適切地描述移動邊界問題。Hirt 等[10]則採用了「Arbitrary Lagrangian-Eulerian Method」的名稱來定義此種混合座標的數值方法，並利用有限差分(finite difference)形式處理可壓縮流體，對於 ALE 法的使用有初步的描述，且對該法的穩定性、準確性與網格劃分方式有詳細的探討。Hughes 等[11] 則率先將ALE 法應用在有限元素(finite element)的分析上來解析黏性不可壓縮流之自由表面問題，對於此法所牽涉到的三種座標系統與座標間的轉換關係有著深入的說明，並且驗證了ALE 法使用於自由表面的波動問題有較佳的收斂解。Belytschko 等[12]利用此一方法來計算流體與結構物表面間的振動問題，使得此一類問題得以獲得較為正確的數值模擬方法。

隨著ALE 法逐漸發展完整，以 ALE 法解決相關工程問題的文獻也逐漸增加。在自由液面的流體力學問題方面，Huerta 與 Liu [13]以 ALE 法配合網格重建求解較大位移量的自由邊界問題，結果證明 ALE 法可適切地模擬自由液面產生的流體激濺(sloshing)現象。在固體與流體間相互影響的問題上，Donea 等[14]、Huerta 與 Liu [15]、Nomura 與 Hughes [16]、及 Nomura [17]等皆採用 ALE 法模擬流體流經圓柱之流場，結果發現圓柱表面產生的流體漩渦剝離(vortex shedding)現象，將造成圓柱振動。Masud 與 Hughes [18] 以ALE 法模擬流體中移動區域(moving domain)問題，並且以靜止流體中運動圓柱與在海面下發射飛彈之潛水艇為例，探討固體與流體間的相互作用。在材料成型與固體力學方面，Liu 等[19]採用葛拉金(Galerkin)之有限元素法於 ALE 中來計算鍛造成形的過程，

並將所得之結果與採用拉格朗日座標系統所得之結果相較，發現採用 ALE 法時，節點隨時間改變後的分布較為平滑，且所獲得的結果也較接近實際的材料澆鑄過程。基於上述的討論得知，ALE 法可適切地處理移動邊界問題。ALE 法將統御方程式定義於另一參考座標上，並且根據網格的移動定義網格移動速度，因此網格速度與流體質點的速度無關，可以避免網格過度扭曲。本研究針對振動發熱源對自然對流與混合對流的影響作一深入的探討，期望藉由本文的分析有助於 IC 等發熱電子元件在振動中的熱傳現象有更清楚的認知，並藉以找出自然對流、混合對流與強制對流模式下，當振動高溫面之振動頻率達到某一大小時，其熱傳效果會高於無振動高溫面熱傳效率的臨界振動頻率F 。 _cc 據此，除本章外，本文尚分為：第二章、ALE 法之運動學理論與數學模式；第三章、

(18)

垂直管道流中振動對高溫壁面熱傳增益之限制研究；第五章、總結與建議。文末另包括自然對流與混合對流模式之有限元素法詳細的離散化推導程序。

(19)

(20)

(21)

(22)

第二章 ALE 法之運動學理論與數學模式

Arbitrary Lagrangian-Eulerian (ALE)是結合拉格朗日(Lagrangian)與尤拉(Eulerian)座標系統所衍生出的座標系統，能有效的處理移動邊界問題，並且在計算的過程中控制網

格的變形。本章節將詳細說明 ALE 法的運動學理論、數學模式、格點速度以及應用方

式等。

2.1 網格描述與座標系統

現今數值模擬的技術中，不論利用有限差分法（finite difference method）或有限元素法（finite element method）離散計算區域，描述流體運動的參考座標有兩種，一為拉格朗日(Lagrangian)座標系統，另一為尤拉(Eulerian)座標系統。其中，採用拉格朗日 (Lagrangian)座標系統所計算之流場，其觀察流場的焦點會集中在各個流體質點上，因此數值計算過程中，其網格會隨著節點上的流體速度而移動；另外，若是採用尤拉(Eulerian) 座標系統來計算流場速度時，則觀察流場的焦點會固定在空間中，所以網格不會隨著流體流動而改變位置。通常模擬本文所探討的移動邊界等類似問題，多會採用拉格朗日 (Lagrangian)座標系統，但是當移動邊界所造成的流場之外另有其他流場存在時，則容易產生誤差乃至於無法求得收斂解，因此利用 ALE 的方法將可更加有效的處理此一移動邊界問題。

2.2 座標定義域

在使用 ALE（Arbitrary Lagrangian-Eulerian）法計算移動邊界問題的座標中，有三種座標定義域的關係存在，分別敘述如下所示，亦可參考圖2-1 之描繪[11,16]： 1. 空間定義域（spatial domain，Ω_y）代表固定不動的定義域，因此在數值方法中為固定的座標系統，其為一般流力問題

(23)

計算中，此定義域即為流體本身，其座標系統稱為拉格朗日(Lagrangian)座標系統。若是存在某座標點z ，則其與空間定義域之座標點_i y 的關係如下： _i

( )

i t

( )

i i z ,t z y =Φ =Φ (2-1) 3. 參考定義域（referential domain，Ω_x）在探討移動邊界問題的領域中，此定義域會移動，通常稱之為參考（Referential）座摽系統。若是存在某座標點x ，則其與空間定義域之座標點_i y 的關係如下： _i

( )

i t

( )

i i ˆ x ,t ˆ x y =Φ =Φ (2-2) 而物質座標點z 與參考座標點_i x 的關係則為 _i

( )

i t

( )

i i z ,t z x =Ψ =Ψ (2-3) 且 r -1 t t Φˆ Φ Ψ = • (2-4)

2.3 速度定義

根據上列所述的三種座標系統與空間定義域，ALE 法亦因此而產生了三種不同的速度定義，分別描述如下： 1. 流體速度（material velocity，u） t z i t y u ∂ ∂ = (2-5)

2. 節點速度（mesh velocity，uˆ）

t x i t y uˆ ∂ ∂ = (2-6) 3. 對流速度（convective velocity，cˆ） t z i t x cˆ ∂ ∂ = (2-7) 且 uˆ -u cˆ= (2-8)

(24)

2.4 座標轉換

另外，由三種不同的座標定義域可以獲得如下的座標轉換關係[20]： 1. 參考(Referential)座標與拉格朗日(Lagrangian)座標之轉換根據微分定理可知，任一函數 f 可以滿足 i t t z i i x z t x x f t f t f ∂ ∂ ∂ ∂ + ∂ ∂ = ∂ ∂ (2-9) 若是把函數 f 以座標函數 y 帶入，即可得 i t t z i i i x i z i t x x y t y t y ∂ ∂ ∂ ∂ + ∂ ∂ = ∂ ∂ (2-10) 亦即 i z i i i t x x y u u ∂ ∂ ∂ ∂ + = ) (2-11) 則

(

)

i i z i y x uˆ -u t x i ∂ ∂ = ∂ ∂ (2-12) 再將其帶入式(2-9)即可獲得最後之座標轉換式

(

)

i x z y f uˆ -u t f t f t t ∂ ∂ + ∂ ∂ = ∂ ∂ (2-13) 2. 尤拉(Eulerian)座標與拉格朗日(Lagrangian)座標之轉換同理，根據微分定理可知，任一函數 f 可以滿足 i t t z i i y z t y y f t f t f ∂ ∂ ∂ ∂ + ∂ ∂ = ∂ ∂ (2-14) 而上式最後之微分項由前面之流體速度的定義式(2-5)可知其恰好為u，故可直接獲得尤拉(Eulerian)與拉格朗日(Lagrangian)兩座標間的轉換關係式如下： i y z y f u t f t f t t ∂ ∂ + ∂ ∂ = ∂ ∂ (2-15)

(25)

再把函數 f 以尤拉(Eulerian)座標函數 y 帶入，即可得 i t t y i i i x i y i t x x y t y t y ∂ ∂ ∂ ∂ + ∂ ∂ = ∂ ∂ (2-17) 亦即 i y i i i t x x y u 0 ∂ ∂ ∂ ∂ + = ) (2-18) 則 i i y i y x uˆ -t x i ∂ ∂ = ∂ ∂ (2-19) 再將上式所得之結果帶入式(2-17)即可獲得參考(Referential)座標與尤拉(Eulerian)座標之轉換關係式如下： i x y y f uˆ -t f t f t t ∂ ∂ ∂ ∂ = ∂ ∂ (2-20) ALE 綜合上述座標系統的概念，將觀測流場的焦點定在新定義的參考定義域 (reference domain)上，此參考定義域以任意速度uˆ在計算區域中運動；計算網格則固定在參考定義域上，與參考定義域同以uˆ的速度移動。ALE 與拉格朗日(Lagrangian)、尤拉 (Eulerian)座標系統間的關係，由uˆ決定：

1. 當uˆ=0時：

參考定義域(reference domain)與計算網格皆靜止不動，此時 ALE 描述流場的觀點與尤拉(Eulerian)座標系統相同。

2. 當uˆ=u時：

此時參考定義域(reference domain)的移動速度與流體質點的速度相同，ALE 描述流場的觀點則與拉格朗日(Lagrangian)座標系統相同。

3. 當uˆ≠0且uˆ≠u時：

此時參考定義域(reference domain)描述流場的觀點即為 ALE 法。

又在實際的數值模擬運用上，無論uˆ=0、uˆ=u或uˆ≠0且uˆ≠u的各種情形，都會同時

發生在同一計算空間內的不同網格節點上，如此彈性的變化，正是利用 ALE 法模擬移

(26)

2.5 ALE 統御方程式

一般考慮浮力項後之二維流場的質量方程式（連續方程式）、動量方程式與能量方程式或濃度方程式分別可以表示為：連續方程式 0 y v x u ₊ ₌ ∂ ∂ ∂ ∂ （2-21）動量方程式       + + − = + + ∂ ∂ y u x u x p 1 y u v x u u t u 2 2 2 2 ∂ ∂ ∂ ∂ ν ∂ ∂ ρ ∂ ∂ ∂ ∂ （2-22） g y v x v y p 1 y v v x v u t v 2 2 2 2 −       + + − = + + ∂ ∂ ∂ ∂ ∂ ∂ ν ∂ ∂ ρ ∂ ∂ ∂ ∂ （2-23）能量方程式       + = + + ∂ ∂ y T x T y T v x T u t T 2 2 2 2 ∂ ∂ ∂ ∂ α ∂ ∂ ∂ ∂ （2-24）濃度方程式       + = + + y c x c D y c v x c u t c 2 2 2 2 ab _∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ （2-25） 2.5.1 強制對流於濃度擴散之移動邊界問題 針對移動邊界問題中，強制對流模式下流場與濃度場擴散的數值模擬，作者於1995 年[21]即已進行此一議題的研究，其變換後的連續方程式依舊相同，如式(2-21)，但動量方程式與濃度方程式則改變如下：動量方程式

(

)

( )

_      + + − = + − + y u x u x p 1 y u vˆ -v x u uˆ u t u 2 2 2 2 ∂ ∂ ∂ ∂ ν ∂ ∂ ρ ∂ ∂ ∂ ∂ ∂ ∂ （2-26）   v v p 1 v v v _∂ _∂ _∂ _∂ 2 _∂ 2 ∂

(27)

(

)

(

)

_      + = + − + y c x c D y c vˆ -v x c uˆ u t c 2 2 2 2 ab _∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ （2-28）又若採用特定之無因次參數後，其無因次方程式則可表示為：動量方程式

(

)

(

)

_      + + − = − + + Y U X U Re 1 X P Y U Vˆ V X U Uˆ -U U 2 2 2 2 ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂τ ∂ （2-29）

(

)

(

)

_      + + − = − + Y V X V Re 1 Y P Y V Vˆ V X V Uˆ -U + V 2 2 2 2 ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂τ ∂ （2-30）濃度方程式

(

)

(

)

_      + = − + − + 2 ₂ 2₂ Y C X C Sc Re 1 Y C Vˆ V X C Uˆ U C ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂τ ∂ （2-31） 2.5.2 強制對流於熱傳增益之移動邊界問題 針對移動邊界問題於強制對流模式所進行增加熱傳效益的研究，則在後續的揚 [22-26]、王[27-30]以及童[31-35]等都有極為詳細的分析，揚為利用振動的方式來增加矩形柱的散熱效益，王則利用往復運動的薄塊搭配噴流的作用，以獲得對高溫壁面更好的熱傳效果，童利用相同的解法來探討振動圓柱對渠道內高溫壁面的熱傳效益，以及高溫圓柱加入振動後的散熱效果等。綜合上述問題，其所採用的 ALE 方程式中，質量方程式與動量方程式和前節所述相同，如式(2-21)、式(2-26)及式(2-27)，而濃度方程式則由能量方程式替換如下：

(

)

(

)

_      + = + − + y T x T y T vˆ -v x T uˆ u t T 2 2 2 2 ∂ ∂ ∂ ∂ α ∂ ∂ ∂ ∂ ∂ ∂ （2-32）若採用特定之無因次參數後，其無因次能量方程式則可表示為

( )

) Y θ X θ ( PrRe 1 Y θ Vˆ -V X θ ) Uˆ U ( τ θ 2 2 2 2 ∂ ∂ + ∂ ∂ = ∂ ∂ + ∂ ∂ − + ∂ ∂ （2-33） 2.5.3 自然與混合對流於熱傳增益之移動邊界問題 在經過楊、王及童等人之研究後，有關移動邊界中強制對流模式的熱傳增益之相關問題已有相當深入的探討，且獲得的結果相當有助益，然而針對移動邊界於自然對流與

(28)

混合對流等模式之問題則尚未有所觸及，因此本文乃針對自然對流與混合對流之問題加以分析，其變換後的質量方程式、x方向動量方程式與能量方程式仍然相同，如式 (2-21)、式(2-26)及式(2-32)，但y 方向（重力方向）之動量方程式在考慮浮力項後，則 可變化為：

(

)

( )

g y v x v y P 1 y v vˆ -v x v uˆ u t v 2 2 2 2 −       + + − = + − + ∂ ∂ ∂ ∂ ∂ ∂ ν ∂ ∂ ρ ∂ ∂ ∂ ∂ （2-34）而其無因次化後的動量與能量方程式則會因為流場是處於自然對流或混合對流的模式，而有不同的形式產生，通常在自然對流的模式下，因為流場不存在入口速度，亦即不採用雷洛數Re為無因次參數下，多以雷利數Ra或葛瑞秀夫數Gr為無因次參數。而在混合對流模式時，即會加入雷洛數Re之定義，而改變無因次方程式之形式。有關自然對流或混合對流下，詳細的無因次參數使用與無因次化後的動量與能量方程式會在後續的章節中加以說明。 2.5.4 加入多孔性介質於熱傳增益之移動邊界問題 連[36]將移動邊界的問題擴展到加入多孔性介質後熱傳增益的影響之研究，其變換後在沒有多孔性介質存在空間（外流場）其連續方程式、動量方程式與能量方程式依舊相同，如式(2-21)、式(2-26)、式(2-34)及式(2-32)，但存在多孔性介質的空間（內流場）其動量方程式與能量方程式則改變如下：動量方程式 u u K F u K y u x u x p 1 u y ) vˆ v ( u x ) uˆ u ( t u 2 2 2 2 v ε νε ν ρ ε ε ε ε − −       ∂ ∂ + ∂ ∂ + ∂ ∂ − =       ∂ ∂ − +       ∂ ∂ − + ∂ ∂ （2-35） ε ε νε ν ρ ε ε ε ε g v u K F v K y v x v y p 1 v y ) vˆ v ( v x ) uˆ u ( t v 2 2 2 2 + − −       ∂ ∂ + ∂ ∂ + ∂ ∂ − =       ∂ ∂ − +       ∂ ∂ − + ∂ ∂ v （2-36）

(29)

其中，uv 為多孔性介質內的平均流速，ε為多孔性介質之等效孔隙率，且

(

₁₅₀ _ε1.5

)

1.75 F = （2-38）

(

)

2 2 p 3 ε 1 150 d ε K − = （2-39） Cp k_e e _ρ α = （2-40）又 d 為顆粒直徑，_p * t 0 e e k k k = + ，其中k 為停滯熱傳導係數（_e0 Stagnant conductivity）， * t

k 為熱散逸傳導係數（thermal dispersion conductivity）。同

樣地，其無因次化後的動量與能量方程式會因為流場處於自然對流或混合對流的模式而有不同的形式產生，本文則不加以探討。

2.6 格點速度

在 ALE 法中，參考座標系統的移動速度是可以任意給定的，因此計算網格的移動速度與流體速度間彼此獨立，可以根據計算的需要給定適當的網格速度。本研究依據不同的物理邊界條件，將網格速度分為三種區域： 1. 固定邊界上，網格速度為 0。 2. 移動邊界上，網格速度等於移動面速度，亦與移動邊界上的流體等速。 3. 在固定邊界與移動邊界間，網格速度可以依需求設定，以避免網格過度扭曲導致計算發散。又本文根據格點與壁面間的距離，使網格速度成線性分布，有關網格的分布與高溫面振動後的網格變化情形可以參考圖2-2 所示，在振動高溫面及其相鄰的上下兩個振動緩衝區的高度內，為網格變化區域，除此區域外，其於的網格不會隨著振動而改變。另不同的網格分布對相同的問題所獲得的數值解之差異性，在王[30]的論文中有極為詳細的探討，其以兩組不同之網格速度，計算振動矩形發熱方柱對管道流之溫度與速度場的影響，藉以探討任意選取網格速度，對計算結果的影響。其物理模式如圖2-3 所示，兩組網格的變化如圖2-4 所示，最後的誤差如圖 2-5 所示。結果顯示，兩網格 U 值差異皆在0.05%以下，V 值的相對誤差值範圍在-0.0405%至 0.0232%之間，又網格 A 與

(30)

此得知，以 ALE 法求解移動邊界問題時，不同的網格移動速度對結果的影響有限，但是仍需選取適當的網格速度，以避免網格在移動的過程中過度扭曲，影響結果。

(31)

(32)

(33)

圖2-3 網格變化研究之物理模式圖[30]

3.8

0.4

insulated wall heated moving block flow inlet A B C D E F G H 0.2 10 0.4 0.2 flow outlet x y

(34)

(a) 網格 A（大區域變化）

(b) 網格 B（小區域變化）

(35)

(a) U 值之等相對誤差線分佈 (b) V 值之等相對誤差線分佈 NuY (c)局部紐塞數分佈圖2-5 兩組網格速度之等相對誤差線與局部紐塞數分佈圖[30] position 0 40 80 120 160 G E F H case A case B

(36)

第三章自然對流垂直管道流中振動對高溫壁面熱傳增益之研究

3.1 前言

過去，針對自然對流之振動分析的研究雖然常有所見，但是大多屬於密閉環境（enclosure）的分析，Fu 和 Shieh [37,38]即針對密閉環境中變動加速度和振動對自然對流的影響加以探討，並發現簡明的修正公式以預測共振頻率及紐塞數的變化。Ishida 和 Kimoto [39]也研究了密閉環境中垂直振動對自然對流的影響，並且發現了振動頻率與紐塞數的關係。Forbes 等 [40]則利用實驗的方法來探討密閉環境中振動對自然對流的影響，而在其實驗中利用振動的方法最多可以增加50%的熱傳速率，此外，由於物體振動所引起的流體波動等幾個類似的自然對流與強制對流的研究亦在該文獻中有所呈現。 Ichioka 等 [41]成功的使用 CFD 技術來模擬兩個振動圓柱與多個振動圓柱列所引起的流體振動問題，其所造成的渦流剝離與流場不穩定現象皆可清楚的觀察到。Baxi 和 Ramachandran [42]利用實驗的方法來分析振動圓柱對自然與強制對流的影響，當雷諾數小於200 時，其對熱傳速率的影響變動不大，但是當雷諾數大於 200 時，振動對於熱傳速率的提升將有很明顯的影響，甚至可以比沒有振動情形下的自然對流高達七倍的熱傳速率。Krishnan 與 Rao [43]利用實驗的方法來探討振動對兩根管路間熱交換現象的影響，當振動頻率達 780-840cpm 時，振動頻率與熱傳效果沒有明顯關聯，但是當頻率繼續增加時，則熱傳速率會有微量的增加。Ivanova 與 Kozlov [44]亦針對水平圓柱之垂直振動對葛雷修數（Grashof numbers）的影響作一深入研究。Fu 與 Tong [31,34]則針對管道中振動圓柱對強制對流流場的影響作一探討，並發現了振動對於熱傳速率的效果明顯增加。本研究的主要熱傳機構乃是基於一個振動的高溫面對於垂直管道流的影響，因為高溫面的振動方向為水平方向且位於垂直管道的左側壁面，因此會形成一個凸起與凹下的反覆運動。當高溫面凸起時，高溫的流體很容易從壁面上分離，也導致熱傳效率提高，但相反地當高溫面凹下時，則高溫的流體將因為流體的停滯而不利於分離，除非當振動的頻率或振幅大到足以克服這一屏障（凹下狀態）。為了避免計算的複雜性，當設備存

(37)

3.2 物理模式

本研究所採用的物理模式如圖3-1 所示，主要為一個高與寬 w 的二維垂直管道。其中高溫振動壁面高度為且位於左側壁面，振動面距離上方與下方開口面高度分別為與，又該一振動面的溫度固定為 Th。在垂直管道中除了高溫振動面外，其餘兩側的壁面皆為絕熱狀態，而外面的環境溫度亦固定為 Tc。當振動開始時（t>0），高溫面將以垂直於重力的方向作一水平來回振動，其振動頻率為且振幅為。由於高溫面於垂直管道流中振動的交互作用影響，流場與溫度場的變化屬於隨時間而變化的移動邊界問題，因此可以ALE 方法適當的描述。

3.3 統御方程式與邊界條件

本研究利用以下的假設條件簡化問題。 1. 流場為二維的層流流場。 2. 流體為不可壓縮之牛頓流體。 3. 流體的各種性質除了密度外皆為常數。 4. 流體與壁面間符合無滑動條件(no-slip condition)。根據_l₂、Th、Tc為參考變數而定義出以下之無因次參數，並將 ALE 統馭方程式簡化如下：

(

)

_, , , 2 2 2 2 l l l l α τ = = ⋅ = = t Ra t t y Y x X f , , 2 2 l l α α = ⋅ ⋅ = Ra v V Ra u U

(

)

, , ˆ ˆ 2 2 2 l l ρ α α = _⋅ ⋅ = Ra p P Ra u U , 2 l lc c L = , , F 2 2 c c h c c T T T T Ra f − − = ⋅ = θ α _l l _Pr _, _r

(

)

_, 2 3 2 ν β α ν _G ₌ g Th −Tc l =

(

)

αν β 3 2 Pr g Th Tc l Gr Ra= ⋅ = − （3-1）

其中uˆ為格點速度（mesh velocity），_{l 與}_c f 分為高溫面振動幅度與振動頻率。_c

(38)

0 Y V X +∂ = ∂ ∂ ∂U _（3-2）動量方程式

(

)

Y U V X U U U U ∂ ∂ ∂ ∂ τ + − + ∂ ∂ _ˆ       + + − = Pr 2₂ 2₂ Y U X U Ra X P ∂ ∂ ∂ ∂ ∂ ∂ _（3-3）

(

ˆ

)

Y V V X V U U V ∂ ∂ ∂ ∂ τ + − + ∂ ∂ _θ ∂ ∂ ∂ ∂ ∂ ∂ ₊ _⋅       + + − = Pr ₂ Pr 2 2 2 Y V X V Ra Y P _（3-4）與能量方程式

(

)

_      + = + − + ∂ ∂ 2 2 2 2 1 Y X ˆ Y X Ra V U U ∂ θ ∂ ∂ θ ∂ ∂ ∂θ ∂ ∂θ τ θ _（3-5）又無因次高溫壁面振動速度U_c 的計算方式如下所示：

(

π _cτ

)

m c U F U = cos 2 （3-6）其中無因次的最大振動速度U 則可以由下式獲得：_m c c m F L U =2π （3-7）初始條件以高溫壁面靜止不動之穩態速度與溫度場為初始條件（固定狀態）。邊界條件在垂直管道下方開口面AB 為 0 , 0 , 0 = = ∂ ∂ = ∂ ∂ _θ Y V Y U （3-8）除了高溫振動面EF 外，其餘的壁面為 0 , 0 V , 0 = ∂ ∂ = = X U θ （3-9）在垂直管道上方開口面CD 為 0 , 0 , 0 = ∂ ∂ = ∂ ∂ = ∂ ∂ Y Y V Y U θ （3-10）

(39)

3.4 數值方法

本研究的數值方法採用葛拉金有限元素法(Galerkin finite element method)。所有元素均為八節點二次等參元素，並配合處罰函數(penalty function)與連續方程式消去壓力項，利用後項差分隱式法(backward different implicit method) 處理時間微分項， Newton-Raphson 迭代法處理非線性項。為縮短計算所需的時間，本研究採用鋒面法 (frontal method)，配合 Gauss Jordan 消去法與高斯積分法，求解經過上述步驟產生之聯

立方程組，詳細的離散方程式與矩陣可參考附錄A 所示。又數值計算流程圖如圖 3-2，詳細步驟如下： 1. 根據所需的元素數目、網格數目與分佈情況，求出各網格點的位置，並指定邊界條件。 2. 求高溫振動面靜止時的速度與溫度場，作為其振動時的初始條件。 3. 計算各網格點的網格速度，並檢查初始條件與邊界條件有無錯誤。 4. 求所需的參數(Jacobian 矩陣、形狀函數等)。 5. 反覆解聯立方程組求出速度與溫度場，直到每一網格點的速度與溫度值滿足下列收斂條件。 3 1 m m 1 m 10− + + < − ϕ ϕ ϕ （3-12）式中ϕ表示U、V 與 θ。 6. 在每一個計算時間裡，均需檢查每一元素和整個計算域的連續方程式的殘值（ Y V X U sidual Re ∂ ∂ + ∂ ∂ = ）以確保在整個計算過程中均能滿足質量守恆。 7. 重複步驟5,6 至達到所需的無因次時間。

3.5 結果與討論

本研究所採用的工作流體為空氣，其普朗特常數為Pr=0.71。主要的研究參數包含雷利數Ra、振幅Lc以及振動頻率F ，本文所採用的參數數值如表 3-1 所示。 _c

高溫振動面的局部紐塞數（local Nusselt number）定義如下：

θ

(40)

高溫振動面的平均紐塞數（average Nusselt number）定義如下：

∫

+ = 3 2 3 2 1 L L L NuYdY L Nu （3-14）

高溫振動面振動一週期之時平均紐塞數（time-averaged Nusselt number）定義如下：

∫

= P 0 P Nud 1 Nu τ τ τ （3-15）其中τ_p代表振動一個週期的時間亦可等於1 F_c。 3.5.1 格點與時間步進測試 為了達到合理的邊界條件，在經由數值測試的結果於高溫面距離頂部與底部的長度分別採用了L₁ =81與L₃ =21，管道之寬度W =2、高溫面長度L₂ =1、高溫振動面與上下固定管道壁面之振動緩衝區長度各為0.5。另為了獲得較準確的速度 U 和 V 以及溫度 θ ，在雷利數_Ra₌₁₀3_{之穩定流場下，採用四種非均勻分佈的計算網格進行網格測試，} 其元素個數分別3840、6848、10640 和 18000，沿著高溫面中心平行 X 軸之 U、V 和θ 分布如圖 3-3 所示。根據網格測試的結果，本研究採用 3840 個元素之計算格點。又在時間步進∆τ 大小方面，在雷利數Ra=104、高溫面振幅L_c =0.1以及高溫面振動頻率 500 F_c = 的條件下（τ_p =1/F_c =1/500=2.0×10−3），分別以∆τ =2.0×10−4（N_P =10）、 1.0×10-4（N_P =20）、0.5×10-4（N_P =40）以及 0.25×10-4（N_P =80）四種不同的時間步進大小測試，結果如圖 3-4 所示，針對此一條件下採用_∆_τ ₌₀_.₅_×₁₀−4_{為時間步進大} 小。但實際的時間步進值並不一定，而會因為雷利數Ra 與振動幅度 Lc的不同而有所改變，基於本測試結果在其他不同條件下的週期步進次數N_P（振動一週期的時間步進次數）皆採用40，而其相對應的時間步進大小可計算如下： c c P F F N × = × = ∆ 40 1 1 τ （3-16）

(41)

Churchill 和 Chu [46]之研究，其與本研究的數值分析結果比較如圖 3-5 所示，而本研究與兩者間會有些微差距的原因在於本題目的高溫範圍僅限於垂直壁面中的振動段內，但另兩者的分析中則以整面垂直壁面皆為高溫段所致，又不同雷利數下的高溫面平均紐塞數分布如表3-2 所示。又為了更加清楚的顯示流場與溫度場的變化情形，所以有關流線與等溫線的顯示範圍僅以高溫面附近為主，並利用箭頭符號來表示振動高溫面的運動方向，而其無因次流線函數Ψ 定義如下： X V , Y U ∂ ∂ = ∂ ∂ = Ψ Ψ （3-17）圖3-6 所示即為雷利數_Ra₌₁₀4_{、高溫面振幅} ₌₀_.₁ c L 、高溫面振動頻率F_c =500下第十一個振動週期的流線與等溫線隨時間變化情形。圖 3-6(a)為第十一振動週期的初始狀況（亦為第十振動週期的結束狀況），無因次時間為τ =0.02，此時高溫面位於週期振動範圍的中間且運動方向朝右，該一時間的運動速度最大，靠近高溫面的流體因為擠壓而向右運動，所以流線皆由高溫面向右延伸，又因為自然對流所引起熱浮力的影響而向上。圖3-6(b)，高溫面已經移動到最右邊的位置，此時高溫面為靜止狀態亦即速度為零，鄰近高溫面的區域其流線仍保持至高溫面向右上方延伸的形式，但在鄰近區域以外則與穩態流場下的流線較為類似。圖 3-6(c)，高溫面又回到位於週期振動範圍的中間但運動方向朝左，該一時間的運動速度亦為最大，不同於圖 3-6(a)的擠壓現象，流體會被高溫面所吸引過來，導致流線直接從下方連接到高溫面。圖3-6(d)，高溫面移動到水平振動的最左邊位置，此時高溫面亦為靜止狀態，在其遠離高溫面區域外的流線與圖3-6(b)或穩態流場下的流線類似，直接由下方邊界連接到上方邊界，但在鄰近高溫面的區域其流線則與圖 3-6(b)相反，並非由壁面外推而是由下方邊界連接至高溫面的吸引模式。圖 3-6(e)所示則為與圖 3-6(a)所示之高溫面位置相同，為第十一振動週期的結束狀況，也因此其流線的分布形式幾乎與圖 3-6(a)相同。在等溫線的分布變化方面，其分布多集中在靠近振動高溫面的區域內，亦即該範圍內的溫度變化較大，又由圖3-6 內的各時間等溫線分布可知，其溫度分布變化不大，並不會隨著高溫面的振動位置而有大幅度的改變，單以等溫線的分布圖無法判別其差異性。

(42)

3.5.3 振動頻率及雷利數對熱傳之影響 圖3-7 所示為當雷利數_Ra₌₁₀3_{且高溫面振幅} ₌₀_.₁ c L 時，平均紐塞數Nu 隨時間τ 的變化情形。由於高溫面的來回振動現象，導致平均紐塞數與會有類似的週期變化發生，且當振動頻率F 越大時紐塞數也會隨之變大。當振動頻率_c F_c =50時，大部分時間下的平均紐塞數都會比固定狀態下（高溫面無振動情形，stationary state）的紐塞數要小，但是當振動頻率F_c =200時，平均紐塞數的大小就完全會比固定狀態下要來的大。另由曲線圖中可以發現，當高溫面凸起時，平均紐塞數會位於曲線的波峰部分，相反地，當高溫面向左運動造成凹陷時，平均紐塞數會位於曲線的波谷段。同樣地在雷利數Ra=103且高溫面振幅L_c =0.1下，時平均紐塞數Nu隨時間的變化情形則如圖3-8 所示，由該圖可以更加清楚的發現，當振動頻率F 為_c 200 時，時平均紐塞數Nu的數值已經大於固定狀態下的平均紐塞數。但是當振動頻率F 等於_c 100 時，時平均紐塞數的大小仍會隨著時間遞減而趨近於固定狀態下的數值，又當F_c =50時，除了初始的幾個振動週期外，其時平均紐塞數皆明顯的小於固定狀態下的平均紐塞數。圖 3-9 所示為雷利數_Ra₌₁₀4_{且高溫面振幅} ₌₀_.₁ c L 時，時平均紐塞數Nu隨時間τ 的變化情形。在任一相同振動頻率F 下，當雷利數_c Ra 增加時，時平均紐塞數Nu亦會隨之增加，但必須補充說明的是，即使在相同的固定狀態下（stationary state），當雷利數Ra 增大時，平均紐塞數 Nu 也會隨之增加（參考圖 3-5 所示）。因為振動高溫面的運動方向與自然對流方向垂直，所以當高溫面處於凹陷情形時，將不利於熱傳的發展。類似上面圖 3-8 所述，當振動頻率F 增加時，時平均紐塞數_c Nu整體趨勢亦會增加，但仍有少數例外的狀況，如F_c =500且τ >0.005時，其時平均紐塞數Nu 就比振動頻率 1000 = c F 時要來的大。在相同的振動幅度L_c

(

=0.1

)

下如果雷利數 Ra 增加到_{10 （如圖 3-10 所示），則振}5 動頻率F 必須超過_c 3000 以上才可使時平均紐塞數Nu大於固定狀態下的平均紐塞數 Nu，這比先前兩個雷利數_Ra ₌₁₀3_與_Ra₌₁₀4_{的需求振動頻率要高上許多。综合上述各}

(43)

約為200、1500 以及 4000。 3.5.4 振動幅度對熱傳之影響 圖3-11 所示為當雷利數_Ra₌₁₀3_且振幅 ₌₀_.₂ c L 時之平均紐塞數Nu 隨時間τ 的變化情形。此時的振動幅度L 等於_c 0.2，這代表在相同的振動頻率F 下，高溫面的運動速度_c 與位移幅度將比前述的分析範例要增加，這樣的結果也導致平均紐塞數的變動範圍更大，這可以由圖3-9 當_Ra₌₁₀3_且 ₌₀_.₁ c L 的結果比較後可知。相同的情形，參考圖3-12 所示，當振動頻率F 等於_c 70 時其時平均紐塞數Nu才會大於固定狀態下的數值。而該一振動頻率（F_c =70）比圖3-8 所示的振動頻率（F_c =200）要小。當雷利數Ra增加到104_與₁₀5_{時，類似的分析在圖}_{3-13 與圖 3-14 中有清楚的展現，} 此一造成時平均紐塞數會大於固定狀態下的振動頻率F 仍然會隨著雷利數的增加而變_c 大，但是該一振動頻率依舊會小於L_c =0.1時的需求頻率（如圖3-9 與圖 3-10 所示）。圖3-15 為當高溫面振幅L 增加到_c 0.4 時，當雷利數_Ra₌₁₀3_{之平均紐塞數}_{Nu 隨時} 間τ 的變化情形，此時的高溫面振動速度與振動幅度將大幅增加，而平均紐塞數的變動性亦會隨著高溫面的振動而更加劇烈地震盪。同樣地隨著振動幅度L 的增加，如圖_c 3-16，當雷利數_Ra₌₁₀3_{時，振動頻率} c F 只要 等於60 即可使時平均紐塞數Nu大小值超過固定狀態下的數值，也更明顯小於圖3-8 所示之需求振動頻率F_c =200。類似的情況在雷利數Ra=104和10 亦可發現如圖 3-17 與5 圖3-18 所示，其大小分別為F_c= 300 與 900，這可相對應於圖 3-13 與圖 3-14 所示之需求振動頻率F_c= 600 與 1500。 3.5.5 臨界振動頻率與經驗公式 一般來說，當高溫表面受到振動影響後都會以為熱傳速率必定會增加，但是經過本文的研究後可以發現，不同的振動頻率F 與振幅_c L 組合後，其熱傳速率仍然會有低於_c 固定狀態之熱傳速率的可能。基於上述的研究在不同的雷利數Ra 與振幅L 下都可以發_c

(44)

只要高溫面的振動頻率超過該一振動頻率時，當振動頻率愈大，得到的熱傳效果也會愈好，而此一振動頻率即可定義為臨界振動頻率F （_cc critical oscillation frequency）。

圖3-19 中的每一個點狀符號代表每一個不同雷利數 Ra 與振幅L 組合後所取得的臨_c 界振動頻率數值，而利用這些點狀符號的位置與數值的方法歸納後，即可獲得一經驗公式如式3-17 來預測不同情形下的臨界振動頻率F 數值。_cc

( ) (

0.44 0.75

) ( )

log

(

4.9 2.7 0.4

)

log F_cc = − L_c+ Ra + L_c2− L_c+ （3-18）圖3-19 所示之直線可以稱為臨界振動頻率預測線，此一直線可由式（3-17）推導而得，其中雷利數的變動範圍可由10 到1 _{10 。由圖 3-19 即可獲知，當設定一振幅}5 c L 後， 在任何雷利數Ra下只要振動頻率會落在該一臨界振動頻率預測線上方，則代表該一狀況下的熱傳速率必定會大於固定狀態的熱傳速率，換句話說，在該區域內的高溫面熱傳速率會比沒有振動情下的高溫面之熱傳速率增加。

3.6 結論

本章利用數值分析的方法來探討垂直管道流中高溫振動面對於自然對流的影響，藉由振動頻率、振幅以及雷利數的變化以找出相對應的結果。其主要的結論可歸納如下所述： 1. 在相同的雷利數下，高溫面在不同的振動頻率與振動幅度之組合下，其熱傳速率有可能會小於固定不動的狀態。 2. 當振動幅度與振動頻率增加時，高溫面的熱傳速率亦會相對的增加。 3. 臨界振動頻率可由雷利數與振動幅度兩項訂定，當雷利數與振動幅度相同時，只要振動頻率大於臨界振動頻率時即可獲得更好的熱傳效果。

(45)

表3-1 設計參數組合表

參數數值範圍

雷利數Ra 101、102、103、104、105 高溫面振幅L_c 0.1、0.2、0.4

(46)

表3-2 穩態非振動高溫面平均紐塞數與雷利數之關係表雷利數 Ra 高溫面平均紐塞數 st Nu 101 1.19105 102 2.14185 103 3.48328 104 5.61508 105 9.61796

(47)

(48)

開始建網格邊界條件初始條件網格點位置網格是否正確否結束計算形狀函數與 Jacobian 矩陣解速度與溫度場速度與溫度場是否收斂否是是輸出資料是否繼續計算移動網格點是否需要網格重建是計算新網格點之速度溫度值否否是

(49)

0 0.4 0.8 1.2 1.6 2

X

-2 -1.5 -1 -0.5 0

U

Elem=18000 Elem=10640 Elem=6848 Elem=3840 0 0.4 0.8 1.2 1.6 2

X

0 10 20 30 40

V

Elem=18000 Elem=10640 Elem=6848 Elem=3840 0 0.4 0.8 1.2 1.6 2

X

0 0.2 0.4 0.6 0.8 1

T

Elem=18000 Elem=10640 Elem=6848 Elem=3840 圖3-3 網格測試(_Ra ₌₁₀3₎

垂直通道內振動對高溫面熱傳增益之研究

國 立 交 通 大 學

機械工程學系

博 士 論 文