高浮慣比之混合對流研究

(1)

國立交通大學

機械工程學系

碩士論文

高浮慣比之混合對流研究

Analysis of Mixed Convection

at High Richardson Number

研究生：劉冠蘭

指導教授：傅武雄博士

(2)

高浮慣比之混合對流研究

Analysis of Mixed Convection at High Richardson Number

研究生：劉冠蘭 Student: Kuan-Lan Liu

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

國立交通大學

機械工程學系

碩士論文

A Thesis

Submitted to Department of Mechanical Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master of Science in Mechanical Engineering

June 2011

(3)

i

高浮慣比之混合對流研究

研究生：劉冠蘭指導教授：傅武雄國立交通大學機械工程學系摘要本研究利用數值方法分析可壓縮流在三維垂直管道中的流動及熱傳機制。流場利用有限差分法進行計算，計算方法可分為兩部分：第一部份為非黏滯性項的尤拉方程式採用 Roe 方法計算通量，並且加入 Preconditioning 矩陣，讓程式在計算低速可壓縮流可獲得良好之收斂結果，而程式因為在使用 Preconditioning 時，加入 Artificial time term 時，已破壞了整個統御方程式，因此需使用 Dual time stepping 疊代使其在 Artificial domain 收斂時才能進入下一個真實時階；第二部份為黏滯性項的計算，採用二階中央插分法。在時間項方面則採用 LUSGS 隱式法，利用 LUSGS 疊代以求出下一時階物理量。出口設非反射性邊界條件避免可壓縮流中壓力波的干擾。在許多應用例子中，溫差常常大於

30K，因此 Boussinesq assumption 不適用。本文採用 OpenMP 方法提升運算速度。

由數值計算的結果得知，三維垂直管道之自然對流，因浮力效應往上推升，在出口端產生最大速度，雷諾數(Re=400)之熱傳效應較雷諾數(Re=100、200)之熱傳還差，主要是在出口端雷諾數(Re=100、200)所吸進的流量，大於雷諾數(Re=400)所吸進的流量，而雷諾數(Re=100)從出口端吸入的速度衝擊到出口速度，故熱傳較雷諾數(Re=200)小。雷諾數(Re=100、200)流場之迴流接近入口處，在流場中會有速度反曲點，造成流場的不穩定。

(4)

ii

Analysis of Mixed Convection at High Richardson Number

Student: Kuan-Lan Liu Advisor: Wu-Shung Fu Department of Mechanical Engineering

National Chao Tung University

Abstract

An investigation of heat transfer in a three-dimensional tapered chimney with consideration of the flow compressibility is studied numerically.The finite difference method is adopted and the computational approaches are divided into two parts. One is the Roe scheme applied for the flux of inviscid terms and the preconditioning matrix is added for the efficiency in all speed fields. The other one is the central difference method of second order utilized to solve viscous terms. The temporal term is solved by LUSGS. Non-reflection conditions at the outlet is derived in order to resolve reflections induced by acoustic waves. In many important natural convection problems, the temperature differences are often higher than 30K. Boussinesq assumption is unreasonable. Besides, the OpenMP method is also used to promote the computing efficiency.

By numerical results, there is the greatest flow speed near the outlet in the three-dimensional vertical natural convection pushed upward by buoyancy effect. The enhancement of heat transfer of Reynolds number 400 is worse than the enhancement of Reynolds number 100, 200. It is mainly due to the more flow rate sucked from exterior near the outlet in the case of Reynolds number 100, 200. And the flow sucked from the exterior impacted the flow exiting the outlet in the case of Reynolds number 100. The heat transfer is worse than the case of Reynolds number 200. Besides, the flow field is unstable due to the backflow near the outlet in the case of Reynolds number 100, 200.

(5)

iii 誌謝由衷的感謝指導老師傅武雄教授在這兩年來給予課業和論文上的指導，以及在生活各方面上的關心與照顧，在此謹致最高的謝忱與敬意。同時也感謝機械系諸位師長在課業方面的指導。另外要特別感謝博士班學長李崇綱、黃玠超、王威翔及黃耘在數值模擬計算以及軟體使用上的協助以及指導。還有實驗室同學黃上豪、林佑璁，范仕坤，學弟妹黃崑榕、鄭景木、羅啟修、李世豪在精神上的鼓勵讓我可以順利完成論文。更要感謝感謝父母親及兄弟姊妹在求學路上一路的支持與幫助，不斷的給予耐心與支持，讓我能繼續堅持下去。最後感謝朋友同學親友們的一路上的鼓勵與關心，今日才能順利完成學業。

(6)

iv 目錄中文摘要 ... i 英文摘要 ... ii 目錄 ... iii 表目錄 ... iv 圖目錄 ... v 符號表 ... vii 一、緒論 ... 1 二、物理模式 ... 9 2-1 物理尺寸與分析模式 ... 9 2-2 分析假設及統御方程式 ... 10 2-3 邊界條件 ... 12 三、數值計算模式 ... 14 3-1 統御方程式 ... 15 3-2 Roe scheme ... 17 3-3 MUSCL 法 ... 24 3-4 Preconditioning ... 26

3-5 Dual time stepping ... 32

3-6 LUSGS 法 ... 34 3-7 非反射性邊界 ... 36 四、結果與討論 ... 40 五、結論 ... 88 參考文獻 ... 89 附錄 ... 93

(7)

v

表目錄

表 3-1 精度係數值 ... 25 表 4-1 討論案例之資料 ... 50

(8)

vi 圖目錄圖 1-1 MPI Method 圖 ... 7 圖 1-2 OpenMP Method 圖 ... 8 圖 2-1 煙囪管道之物理模式圖 ... 13 圖 2-2(a~c) XY 及 YZ 平面物理模式圖 ... 14 圖 3-1 黎曼問題特徵值結構圖 ... 23 圖 3-2 差分示意圖 ... 31 圖 3-3 L1、L2、L3、L4與L5於管道兩端的方向性示意圖 ... 39 圖 4-1 網格測試圖 ... 47 圖 4-2 進口流量圖 ... 48 圖 4-3 隨 X 截面出口方向流量變化圖 ... 49 圖 4-4(a~h) 中央 XY 平截面流線圖(Re=0) ... 51 圖 4-5 中央 XY 截面等溫線圖(Re=0) ... 53 圖 4-6(a~d) YZ 截面速度場圖(Re=0) ... 54 圖 4-7 中央 XY 截面流線圖 (Re=950) ... 56 圖 4-8 中央 XY 截面等溫線圖 (Re=950) ... 57 圖 4-9(a~d) YZ 截面速度場圖(Re=950) ... 58 圖 4-10 中央 XY 截面流線圖 (Re=400) ... 60 圖 4-11 中央 XY 截面等溫線圖 (Re=400) ... 61 圖 4-12(a~d) YZ 截面速度場圖(Re=400) ... 62 圖 4-13(a~h) 中央 XY 截面流線圖 (Re=200) ... 64 圖 4-14(a~d) 中央 XY 截面等溫線圖 (Re=200) ... 66 圖 4-15(a~d) YZ 截面速度場圖(Re=200) ... 67 圖 4-16 四面沿 X 方向之紐塞數加總比較圖(Re=200) ... 69 圖 4-17(a~d) 中央 XY 截面流線圖 (Re=100) ... 70 圖 4-18(a~d) 中央 XY 截面等溫線圖 (Re=100) ... 71

(9)

vii 圖 4-19(a~d) YZ 截面速度場圖(Re=100) ... 72 圖 4-20 四面沿 X 方向之紐塞數加總比較圖(Re=100) ... 74 圖 4-21 進口速度分布圖 ... 75 圖 4-22 出口速度分布圖 ... 76 圖 4-23 各雷諾數壁面中央紐塞數比較圖(Re=0，100，200，400，950) ... 77 圖 4-24 各雷諾數角落紐塞數比較圖(Re=0，100，200，400，950) ... 78 圖 4-25 各雷諾數隨時態變化之平均紐塞數比較圖(Re=0，100，200，400，950) ... 79 圖 4-26 各雷諾數之總平均紐塞數比較圖(Re=0，100，200，400，950) ... 80 圖 4-27 中央截面等溫線比較圖(Re=100，200) ... 81 圖 4-28 t=15s時中央截面速度比較圖(Re=0，100，200，400，950) ... 82 圖 4-29 隨時間變化之平均紐塞數(Re=200、100) ... 84 圖 4-30 時間平均之中央截面速度分布比較圖(Re=100，200) ... 85 圖 4-31 時間平均之中央截面雷諾數分布比較圖(Re=100，200) ... 86 圖 4-32 時間平均之中央截面溫度數分布比較圖(Re=100，200) ... 87

(10)

viii

Nomenclature

a sound speed[m s⋅ −1]

p

C constant-pressure specific heat[J kg⋅ −1⋅k−1]

v

C constant-volume specific heat[ 1 1

J kg⋅ − ⋅k− ]

d caliber of the vertical square tube [m ]

e internal energy[ 1

J kg⋅ − ]

g acceleration of gravity[m s⋅ −2]

h enthalpy[J ]

k thermal diffusivity[W m⋅ −1⋅k−1]

l length of the vertical square tube [m ]

x

Nu

Nusselt number defined in Eq.(4-6)

0 0 ( ) ( ) x h d T Nu k T k T T z ∂   = _ _ −  ∂ 

Nu

average Nusselt number with area defined in Eq.(4-7)

0 0 1 1 ( ) ( _h ) A A d T Nu Nudxdy k T dxdy A A k T T z ∂   = = _ _ −  ∂ 

∫∫

( )

Nux t average Nusselt number with time defined in Eq.(4-8)

(

)

0 0 1 ( ) ( ) x t _t h d T Nu k T dt t k T T z ∂   = _ _ −  ∂ 

∫

( )

Nu t

average Nusselt number defined with area and time in Eq.(4-9)

( )

1 t t Nu Nudt t =

_∫

( )

_x y t Nu

∫

Total Nusselt number and average with time defined in Eq.(4-10)

(

)

0 0 1 ( ) ( ) y x t t x h d T Nu k T dxdt t k T T z ∂   = _ _ −  ∂ 

∫

∫ ∫

(11)

ix

2

Re Gr Ri=

Q Total flow in a section defined in Eq.(4-14)

2 0

Q

=

ρ

u d

R gas constant [J kg⋅ −1⋅k−1]

Ra Rayleigh number defined in Eq.(4-11)

2 3 0 0 2 0 ( ) Pr Pr ( ) h g T T d Ra Gr T T ρ µ − = ⋅ =

Re Reynolds number defined in Eq.(4-5)

0 0

0

Re

ρ

u d

µ

=

Ri Richardson number defined in Eq.(4-12)

T Kelvin temperature[K]

0

T Surrounding temperature[K]

h

T _{temperature of heat surface[}_K_]

t

∆ time difference[s ]

u _{velocity component in} _{x -direction[} 1

m s⋅ − ] v velocity component in y-direction[m s⋅ −1]

w _{velocity component in} _{z -direction[} 1

m s⋅ − ]

x Cartesian coordinate system x direction

y _{Cartesian coordinate system y direction}

z Cartesian coordinate system z direction

w

(12)

x Greek symbols ρ density[kg m⋅ −3] 0 ρ surrounding density[ 3 kg m⋅ − _] υ kinematics viscosity[m2⋅s−1] µ absolute viscosity[kg.m−1.s−1]

(13)

1

第一章

第一章、

、

、緒論

緒論

在科技日新月異及高污染的時代下，如何減少能源的損耗是重要的課題，為使能源永續發展，管理能源的使用，是必須要關注的；消費性電子產品蓬勃發展，高功率電子元件的發展與半導體製程技術的進步，促使消費性電子產品的使用與生活型態做結合，因而目前許多電子產品皆走向高效能、微小化的趨勢，致使電子元件舉凡電腦 CPU，顯示卡，筆記型電腦或智慧型手機等其封裝元件之發熱密度愈來愈高，其單位熱通量亦相對地不斷增加，對性能及可靠度等方面均造成不容忽視之影響。因此在未來電子產品的發展趨勢走向更輕薄短小之際，其性能及可靠度的提升將取決於其散熱技術。一般正常電腦 CPU 在執行程式下溫度大約為攝氏 50~60 度之間，超過攝氏八十度則電腦會基於保護狀態，自動關閉系統，長期下來內部零組件壽命將會減短。因此為了維持元件於額定溫度下運作，必須將此密集的熱量能有效散逸於系統外之環境。傳統的散熱片大多數是由鋁合金所製造，其熱傳導性只屬於中等程度，對於目前元件發熱功率越來越高的情況，逐漸有捉襟見肘的情形發生，目前電子元件的發熱量達到每平方公分數十瓦的等級，且接點可承受的溫度約在攝氏 150 度以下。因此如何改良空氣對流的方式以及增加對流所能散逸之熱量成為研究的主要課題之一。在工程應用上，熱對流一直是很重要的一環，其中以混合對流最接近生活中常見的對流形式；一般自然對流效果不足時，再加進強制對流使之達到預期之效果。以自然對流之浮力方向及強制對流流場之速度方向，可分為三種型式如下： 1. 助流(Aiding flow)：其浮力方向與流場之速度方向相同。 2. 逆流(Opposing flow)：其浮力方向與流場之速度方向相反。 3. 交錯流(Cross flow)：其浮力方向與流場之速度方向垂直。 Bergles[1, 2]將常見的增加熱傳效率方式，詳細地討論與整理，將其略分為兩大類：一為外作工的被動方式，另一類為需外加能量的主動式方法。討論混合對流方面， Ingham[3]等人研究二維垂直平板等溫壁面混合對流之逆流現象(flow reversal)，忽略軸向擴散，討論−300≤Gr Re≤70的區間，並與純強制對流作比

(14)

2

較，發現在較高的Gr Re 比中，會有逆流的現象，其設定為Boussinesq assumption。

Behzadmehr等人[4]探討垂直長管等熱通量之低雷諾數助流混合對流，在

8

Re=1000,1500 and Gr≤10 利用k−ε model，其出口之設定為完全展開流(fully

developed flow) ，利用無因次參數Nu及Gr界定層流到紊流的區域，在Re = 1000的設定

下，當 5 4 10 Gr> × ，其入口壓力會小於出口壓力(P₀ <P_L)，且在8 10× 5 <Gr< ×5 107之間屬於紊流；在Re = 1500的設定下，當 5 3 10 Gr> × ，其入口壓力會小於出口壓力(P₀ <P_L)，且在 6 8 2 10× <Gr<10 之間屬於紊流。Barletta等人[5, 6]利用數值計算無限長之長方管等溫及等熱通量之單面及多面壁面，使用Boussinesq assumption去探討層流混合對流在不同長方管比例的逆流現象(flow reversal)。Barletta等人 [7]討論在二維垂直平板混合對流之黏性耗散，左側為低溫壁面，右側為高溫壁面，使用Boussinesq assumption，downward flow之逆流現象較upward flow明顯，其黏性之耗散會增加浮力效應。Suastegui等人[8] 探討二維暫態層流有限管長之逆流混合對流，使用Boussinesq assumption，在一側中間段等溫，其餘為絕熱壁面，研究不同浮慣比在加熱附近的迴流現象。Yang等人[9]研究二維長直管，一面為等熱通量壁面，另一面為絕熱壁面，使用Boussinesq assumption探討浮慣比分別為加溫壁面Ri=100, 400及冷卻壁面Ri= − − −2, 3, 8，並討論Gr Re−Re的關

係，達到Fully Developed或Flow Reversal之區域。Boulama與Galanis[10]探討二維垂直平板，比較等溫及等熱通量壁面之計算，等溫(UTW)的影響參數較單一，等熱通量(UHF) 的影響參數為Gr_T Re 及 Gr_M Re與熱通量率。Desrayaud等人[11]探討二維垂直壁面高溫管道之流體逆流現象，其工作流體為空氣，高溫壁面(Th =60℃)與流體溫差50度，討論 4 6 4.71 10 ~ 1.24 10 Gr= × × 、Re=300 ~ 1300及管道長寬比，研究指出雷諾數(Reynolds number)的增加，會使上述現象消失，亦與高Peclet number (Re PrL⋅ )相關，此現象與管長無關。Zghal等人 [12]研究二維垂直管道，絕熱管壁間有一段加熱壁面，此研究設定

Boussinesq hypothesis ，比較不同加熱管長、雷諾數(Reynolds numbers)及浮慣比

(Richardson numbers)之流線及紐塞數(Nusselt number)，而是否會有迴流的現象主要是由 Peclet number及浮慣比的關係來決定。

(15)

3

Mai 等人[13]討論二維垂直圓柱之不可壓縮混合對流不穩定現象，在加溫部分，不

穩定的渦流會出現在壁面附近，逆流現象出現在波不穩定之下方；在降溫部分，低雷諾數會在頂部產生逆流之坡不穩定現象，其設定條件為 Boussinesq Assumption。Bhoite 等人[14]研究二維薄膜散熱，其條件設定亦是 Boussinesq assumption，文中提到強制對流速度( forced convection velocity υin, fc )是由風扇帶動的等溫流體，自然對流速度(natural

convection velocity υin, nc)則是由浮力效應所帶動，其混合對流進口速度設定則由強制對

流速度給定，且其壓力梯度(∆P)與雷諾數定義亦是與強制對流設定相同，在高雷諾數時會有較強的迴流產生，此時浮力效應可忽略。Nguyen[15]等人研究三維壁面高熱通量圓柱之暫態混合對流，其設定為 Boussinesq assumption，混合逆流(Opposing flow) 在

5 3 10 Gr= × ，出口壁面附近會先出現逆流現象；混合助流(Aiding flow)在Gr=106會先在軸中心出現逆流現象。Koizumi 等人[16]研究在三維空間中有一軸對稱發熱金屬球，其實驗設定為 5 10 3 . 3 × = Gr 與Re≤1900；其模擬是利用 Storm/CFD2000 software 設定條件為 Boussinesq Assumption，溫差 20 K，其浮慣比

(

2

)

Re Gr 在 0.33~0.83 之間，經由實驗及模擬對照可以相輔相成。Stiriba 等人[17]研究三維水平管道，入口處有一下凹處，下凹處有一垂直壁面為高溫壁面，其他為絕熱壁面，在 3 6 10 ~ 10 = Gr 範圍內及 Re= 100, 1000這兩個雷諾數下去比較其速度場、流線、溫度分佈與在高、低溫面上之

平均紐塞數(averaged Nusselt number)，可以得其結果在 Re=100 且(浮慣比)Ri≤100與

Re=1000 且Ri≤0.1時，流場結構會趨於穩定，此研究之條件亦是使用 Boussinesq

Assumption。以上之文獻多數使用 Boussinesq Assumption 作討論。

Laaroussi 等[18]研究二維垂直壁面高溫管道內密度變化，因高溫壁面附近受溫度影

響，密度會改變，應設定為理想氣體，並與 Boussinesq assumption 作討論，利用橢圓模型流允許入口區域有逆流產生，混合對流會因高溫壁面產生的密度改變而迫使流場向軸中心推進，造成逆流。

此外，亦有討論紊流入口之流場之混合對流，Azizi and Benhamou [19]討論在層流

混合對流中，浮力向上或向下對熱傳及質傳的影響，由研究中可得知，在重力與流動方向相反的流動中，存在壁面與進口流體的高溫差，會產生流體不穩定性及紊流產生。

(16)

4

Shahraeeni[20]研究垂直管中混合對流之紊流現象，以定熱通量的條件下觀察浮力效應對

熱傳係數之實驗研究，並利用

κ

−

ε

model 作模擬對照有用的實驗數據，由模擬可得向下

流動之熱傳與浮力效應都較大；對照向上流動之熱傳可增可減，則取決於浮力效應之強度；在實驗數據中顯示強浮力效應時，Nusselt number 會劇降，其結果與模擬相同。Balaji

[21]探討紊流計算時，以 direct numerical simulations (DNS)去計算紊流，而靠近壁面的區

域需要加入一個由流量分析的壁面函數，作討論與比較。綜合以上論述，前人在探討此類混合對流之垂直管道的研究，多數限制較多，其一是數值模擬多假設二維流場，實際三維流場之流動現象因而無法仔細探討；再者是為求符合 Boussinesq Assumption 必須用溫差 30 度以內[22]之限制且流體均視為不可壓縮流，壁面與流體溫度差設定受限，與實際應用之溫差相距甚大；而依目前大部分的計算方式，在計算流場部份依照流體速度將其區分為可壓縮流(大於 0.3 馬赫)與不可壓縮流(小於 0.3 馬赫)，此種區隔卻嚴重的影響其應用範圍。混合對流中給定的速度所帶給流場的質量流率是否能讓熱傳效益增加，是本文主要探討的問題，故此研究著重在高浮慣比

(

2

)

Re Gr 之助流(aiding flow)混合對流。本研究針對三維的流場作模擬且在自然對流流場中加入強制對流，考慮一垂直方型管之混合流場，其重力方向與強制對流流動方向相反，此時浮力效應方向與強制對流流動方向相同，觀察此時流場的熱傳增益及溫度場、速度場變化。因為密度會隨著壓力與溫度的變化而改變，對於實際工業界中的應用較為廣泛與實際，因此在本研究中的工作流體不論高速或低速皆視為可壓縮流，而不使用 Boussinesq Assumption，以增加應用範圍。本研究為三維垂直方形壁面加熱管道，流體在管道中因受到壁面加熱產生溫度梯度，探討此垂直方形管道對高溫壁面附近及出口之溫度及速度場的影響，並說明高溫壁面上之平均紐塞數隨著管道位置而變化，進而比較不同對流方式(自然對流及混合對流)與溫度場之影響及熱傳效率之增益。李[23]發展出黏性流場之全域速度場數值解法，此種方法最大的困難處在於計算低速流場時，由於可壓縮流必須遵守 Courant-Friedrichs-Lewy (CFL) 條件，因此在低速流體時，受限於流體變化傳遞速度(約等於聲速)，時階將會極小。在此種情況造成計算過程將耗費極大的時間與整體過低的計算效率。為了改善此缺點，本

(17)

5

研究在計算可壓縮流時加入 Preconditioning 法，藉此讓流體即使在低速時，也可有較高的效率與良好的收斂性。在計算此種低馬赫數流體的方面，目前有密度基底法

(Density-based method)與壓力基底法(Pressure-based method)。本研究以密度基底法為

主。而在密度基底法中又以 Turkel[24]提出 Preconditioning 法最為廣泛應用，不僅可同時應用在可壓縮流與不可壓縮流中，更可以讓程式的收斂性增加。在邊界設定最棘手的問題是流體的流動速度和壓力波的速度相差過大，導致在出口反射回彈的壓力波會干擾流體流動，故在入口及出口處設置非反射性(non-reflecting)邊界條件[25]，避免干擾發生，若設定為完全展開流，則會與實際流場有所差異。

故本文需要求解完整的 Navier-Stokes 方程式與理想氣體方程式(ideal gas equation) 以得此可壓縮流中密度的變化，以期能同時考慮密度與壓力變化之效應符合實際物理現象。採用中央插分法處理黏性項部分；而非黏性項則必須採用 Roe 法[26]，以解決偏微分方程式的不連續問題；本研究在作數值計算時，主要是利用網格之間的物理量，因此使用 MUSCL(Monotone Upstream-centered Schemes for Conservation Laws)法[27]來計算

ROE 法中網格與網格間的物理量；為了處理在低速可壓縮流時不易收斂的困擾，計算時

加入了 Preconditioning 法[28]，可有效提高收斂效率，使其在高速與低速可壓縮流場均可適用，但卻破壞原統御方程式，為了彌補此一缺點，故需要加入 Artificial time term 來修正；並使用 Dual time stepping[29]來計算暫態的物理量；最後使用 LUSGS(implicit

lower-upper symmetric Gauss-Seidel algorithm)法[30]做迭代運算。

此外，在處理複雜的流體力學問題時需要大量的計算過程，利用多核心處理器來提升運算速度已為目前的發展主流。多核心處理器對於單一執行緒在平行運算方面，並無法提升計算速度；若利用多執行緒的程式架構，可透過不同核心來同時計算，達到提升計算效率、節省計算時間之目的；但是多執行緒的程式在撰寫、編輯上，也都比單一執行緒的程式架構複雜。常見的平行運算方法有 MPI 和 OpenMP 兩種方法；本文採用

OpenMP 方法提升運算速度。MPI (Message Passing Interface)是一種分散式記憶體

(distributing memory)的觀念，而 OpenMP (Open Multi-Processing)為共享式記憶體(sharing memory)的觀念，而兩種方法各有優缺點。MPI 在撰寫程式上面相較於 OpenMP 較困難，

(18)

6 且計算速度較慢且會受限於網路效率，在設備擴充方面較昂貴；反之 OpenMP 在程式撰寫上較簡單，計算速度快且不會受限於網路效率，在設備擴充方面需較少的經費。因此本程式利用 OpenMP 來進行平行化運算，效率為原程式的四倍，減少計算時間及成本。結果發現，三維高溫壁面之自然對流，其隨時間從逆流產生至穩定狀態之變化，因浮力效應往上推升，在出口端產生最大速度。若控制入口流量與自然對流流量相同，因已強制入口之流量，浮力效應無法發揮到最大，其熱傳效果會較自然對流略低一些；若控制入口流量在 1/2 之自然對流流量，在出口端會吸進流量造成逆流，以彌補流場不足的部分，因出口之逆流縮減出口流量，故其熱傳效果較差；若控制入口之流量低於自然對流可吸進量之 1/4，其浮慣比，其出口會大量吸入不足之流量並與強制之速度造成碰撞，造成內部流場之不穩定，其熱傳效果高於自然對流可吸進量 1/4 流量之混合對流。

(19)

7

(20)

8

(21)

9

第二章

第二章物理模式

物理模式

2-1、、、、物理尺寸與分析模式物理尺寸與分析模式物理尺寸與分析模式物理尺寸與分析模式：：：：圖(2-1)為本文之物理模式圖，正方形截面之三維垂直方形管道，正方形截面之邊長為d ，其管道總長度為l，全長皆為等溫加熱壁面。初始溫度T₀的冷卻空氣以均勻流等速度u₀由管道入口(ABCD)進入，在管道壁面為高溫壁面Th，Th對管道內的空氣加溫，進行熱交換；探討自然對流現象時，入口速度u₀為零，流體在管道中因受到壁面加熱溫度上升，產生效應。其中

x

方向為流動方向， y 、皆為垂直壁面方向。z

u

、

v

、

w

分別為其所對應的速度。為簡化管道流研究，以往大多將出口設為速度完全發展流條件及壓力為大氣壓力邊界，但此方式在可壓縮流中壓力波在出口處易會反射回計算區域而影響收斂性，較不適用於此低速可壓縮流流場；且以往為使熱管道內的熱流場更容易成為完全發展流，在管道的長度以及寬度上需要一定的比例，但此一方式會使網格需求增加，計算時間較長，因此本研究在出口條件部分設為非反射性(non-reflecting)邊界[25] 條件，此時可視其出口為遠處的邊界條件，大幅減少網格數量。此研究分別探討純自然對流及高浮慣比之混合對流，討論其質量流率、流場方向及密度分布作分析，並以與自然對流進口流量相同之混合對流、1/2 自然對流流量之混合對流及 1/4 自然對流流量之混合對流做分析比較。

(22)

10 2-2、、、、分析假設及統御方程式分析假設及統御方程式分析假設及統御方程式分析假設及統御方程式：：：：本研究選擇層流流場作為模擬流場，設定如下： 1. 可壓縮流，空氣密度會隨溫度與壓力而改變。 2. 工作流體為空氣，假設為理想氣體。流體性質為牛頓流體(Newtonian fluid)，黏滯係數為等方向性。 3. 有重力效應影響。 4. 進出口條件皆為完全非反射條件。 5. 考慮溫度變化對流體所造成的影響。

6. 所有壁面均為不可滑移(No slip condition)。

統御方程式分別為連續方程式(Continuity equation)、動量方程式(Momentum equation)與能量方程式(Energy equation)，壓力方面則假設流體為理想氣體，利用理想氣體狀態方程式定義。連續方程式、動量方程式、能量方程式如下： S z H y G x F t U ₌ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (2-1) 其中

(

)

T U = ρ ρu ρv ρw ρe (2-2) 2 2

2

xx xy xz xx xy xz

u

P

uv

F

uw

V

T

e

u

Pu

k

u

v

w

x

ρ

τ

ρ

τ

ρ

τ

ρ

τ







_{+ −}









₋



=





−







_

_

_∂



+

−



_

_





_∂











(2-3)

(23)

11                     − − − ∂ ∂ − +       + − − + − = yz yy yx yz yy yx w v u y T k Pv v V e vw P v vu v G τ τ τ ρ τ ρ τ ρ τ ρ ρ 2 2 2

(

)

(

)

               − − − − = gu g S 0 0 0 0 0 ρ ρ ρ ρ

ρ

為密度，p為壓力。u 、v 、w 分別為x、y、z方向的速度。k為 thermal diffusivity。

2 2 2 1 ( ) 2 v e=C T+ u + +v w ，C_v為等容比熱。理想氣體狀態方程式： P=ρRT (2-7) 2 2 2 zx zy zz zx zy zz w wu wv H w P V T e w Pw k u v w z ρ ρ τ ρ τ ρ τ ρ τ τ τ    ₋     ₋  =  + −    _ _ _∂   _ ₊ _ ₊ ₋ ₋ ₋ ₋   _∂      (2-4) (2-5) (2-6)

(24)

12 2-3、、、、邊界條件邊界條件邊界條件邊界條件：：：：本研究所採用的統御方程式為可壓縮 Navier-stokes 方程式，因此需要給定的邊界條件有：初始狀態、入口條件、出口條件、壁面邊界。 2-3.1 初始狀態：初始速度、初始壓力、初始密度初始速度 u： 0 /m s 初始速度 v： 0 /m s 初始速度 w： 0 /m s 初始壓力 p：一大氣壓力(101300Pa) 初始密度ρ ：空氣密度( 3 1842 . 1 kg m ) 2-3.2 入口條件：入口速度 u： u₀ (混合對流)、非反射性邊界(自然對流) 入口速度 v： 0 /m s_{(混合對流)、非反射性邊界(自然對流)} 入口速度 w： 0 /m s_{(混合對流)、非反射性邊界(自然對流)} 入口壓力 p：非反射性邊界入口溫度 T：非反射性邊界 2-3.3 出口條件：出口速度 u：非反射性邊界出口速度 v：非反射性邊界出口速度 w：非反射性邊界出口壓力 p：非反射性邊界入口溫度 T：非反射性邊界 2-3.4 壁面邊界：邊界速度：不可滑移條件，u=v=w=0m/s 邊界溫度：加熱壁面溫度Th 邊界壓力：在垂直壁面方向，梯度為零， p 0 n ∂ = ∂

(25)

13

K

T

Pa

P

0592

.

298 101300

0

=

圖 2-1 垂直管道之物理模式圖

l

d

A

B

C

D

E

F

H

G

0 u

= =

v

w

=

h

T

=

T

0

u

g

(26)

14

第三章

第三章數值計算模式

數值計算模式

本章主旨在說明本論文的數值計算所使用的所有模式。第一節首先要知道所求解的方程式為 Navier-Stokes 方程式。參考 Fu 和 Li 等人[22]之數值方法，本研究將方程式拆解為非黏滯項與黏滯項。其中在計算黏滯性項則採用二階精度的中央差分法。第二節介紹的為黎曼解中的 ROE 法(Roe Scheme)[26]，利用 ROE 法來求出非黏滯項的通量。接著第三節介紹 MUSCL 法 (Monotone Upstream-centered Schemes for Conservation

Laws)[27]，此法是為了要解出 ROE 法中使用的網格之間的物理量，然後為了防止在高

階插分時產生震盪現象，在 MUSCL 法插分的結果方程式中加入 Minmod limiter 以確保程式不會發散。第四節為介紹 Preconditioning 法[28]，因為當計算低速可壓縮流時，因速度和音速的數量級(order)上差距過大，在數值分析時不好計算，所以為彌補此一缺點須使用 Preconditioning 法。第五節為 Dual time stepping，因使用 Preconditioning 法時，加入 Artificial time term 而破壞了完整的統御方程式，為了使計算暫態結果較準確，因此需使用 Dual time stepping 疊代，使其在 Artificial domain 收斂時才能進入下一個真實時階，更提高程式的效率。第六節為 LUSGS 法(implicit lower-upper symmetric Gauss-Seidel

algorithm)[30]，程式因為在使用 Preconditioning 法，故而在統御方程式中修改了計算時

階，故需要加入 Artificial time term，待時階項收斂時才能達到真實穩態。第七節為完全非反射邊界之介紹，為了使自然對流之邊界可以更接近真實狀態，在程式的出入口邊界使用完全非反射邊界。綜合上述，本論文在數值上的計算過程為，利用 MUSCL 法算出

ROE 法所需要的網格間物理量，求出非黏滯的通量，並且在計算通量時加入 Preconditioning 法，以拉近與音速的數量級(order)。接下來使用二階中央插分法對黏滯

項做插分進而求出黏滯項；然後再與 ROE 法求出的通量結合得到真正的物理通量。最後使用 Dual time stepping 及 LUSGS 法疊代以求出正確時階的物理量。

(27)

15 3-1、、、、統御方程式統御方程式統御方程式統御方程式：：：：本研究在計算流場的方面其統御方程式分可為兩大部分，第一部份為非黏滯性項的尤拉方程式，第二部份為黏滯性項。 S z H y G x F t U = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (3-1) 其中

(

)

T U = ρ ρu ρv ρw ρe (3-2) 2 2

2

xx xy xz xx xy xz

u

P

uv

F

uw

V

T

e

u

Pu

k

u

v

w

x

ρ

τ

ρ

τ

ρ

τ

ρ

τ







_{+ −}









₋



=





−







_

_

_∂



+

−



_

_





_∂











                    − − − ∂ ∂ − +       + − − + − = yz yy yx yz yy yx w v u y T k Pv v V e vw P v vu v G τ τ τ ρ τ ρ τ ρ τ ρ ρ 2 2 2

(

)

(

)

               − − − − = gu g S 0 0 0 0 0 ρ ρ ρ ρ 2 2 2 zx zy zz zx zy zz w wu wv H w P V T e w Pw k u v w z ρ ρ τ ρ τ ρ τ ρ τ τ τ    ₋     ₋  =  + −    _ _ _∂   _ ₊ _ ₊ ₋ ₋ ₋ ₋   _∂      (3-3) (3-4) (3-5) (3-6)

(28)

16

ρ

為密度，p為壓力。u 、v 、w 分別為x、y、z方向的速度。k為 thermal diffusivity。

2 2 2 1 ( ) 2 v e=C T+ u + +v w ，Cv為等容比熱。上式可拆解為黏滯性項與非黏滯性項：                 + + −                   ∂ ∂ − +       + + = + = xz xy xx xz xy xx viscid inviscid w v u x T k Pu u V e uw uv P u u F F F τ τ τ τ τ τ ρ ρ ρ ρ ρ 0 2 2 2 (3-7) (3-8)                 + + −                   ∂ ∂ − +       + + = + = zz zy zx zz zy zx viscid inviscid w v u z T k Pw w V e P w wv wu w H H H τ τ τ τ τ τ ρ ρ ρ ρ ρ 0 2 2 2 左式由非黏滯項組成的方程式即稱為尤拉方程式。 2 2 0 2 yx yy inviscid viscid yz yx yy yz v vu v P G G G vw V T u v w e v Pv k y ρ ρ τ ρ τ ρ τ τ τ τ ρ              ₊    = + = −        _ _ _∂     _ ₊ _ ₊ ₋  _ ₊ ₊ _  _∂      (3-9)

(29)

17 3-2、、、、Roe scheme：：：：在雙曲線的守恆形式方程式中，若其初始條件包含有不連續的片段連續(piecewise) 常數，此類型的問題通稱為黎曼(Riemann)問題。因為其包含有不連續解，因此在流體計算上有著相當廣泛的應用。一維線性黎曼方程式如下： 0 u u a t x ∂ ₊ ∂ ₌ ∂ ∂ (3-10) 其中a為一常數Jacobian矩陣。初始條件為

( )

, 0 ₀

( )

0 0 L R u x u x u x u x <  = = >  從新改寫(3-10) 0 U U A t x ∂ ₊ ∂ ₌ ∂ ∂ 求出A之特徵值矩陣以及特徵向量。 1 A= ΛK K− ，其中Λ為特徵值矩陣： 1

0

m λ λ     Λ =      ⋯ ⋮ ⋱ ⋮ ⋯ 。 (1) ( ) , , m T K =_K ⋯ K _ _{為特徵向量，故} ( )i ( )i i AK =λK 。接著定義特徵變數W〈characteristic variables〉，其定義如下： ( , ) W =W t x ，W =K( 1)− U 或U=KW。因此 U K W t t ∂ ₌ ∂ ∂ ∂ 且 U W K x x ∂ ₌ ∂ ∂ ∂ ，將此結果代入(3-10)式中可得： 0 t x KW +AKW = ，可再繼續簡化成： 0 t x W + ΛW = (3-11)

方程式(3-11)稱為canonical form或characteristic form。將以上的結果簡單整理如下：

(30)

18 0 i i i W W t λ x ∂ ₊ ∂ ₌ ∂ ∂ ，或 1 1 1 2 2 3 3 ... 0 0 ... 0 0 : : : : : 0 ... _m t x w w w w w w λ λ               ₊    ₌                         (3-12) 上式可由特徵曲線法求得其解為： (0) 0 ( , ) ( ) 0 i i i i i i i x t w x t w x t x t α λ λ β λ − <  = − = − >  (3-13) 其中，α_i與β_i為初始值的特徵變數。由於U=KW，可以得到 (0) ( ) ( , ) ( ) m i i i i u x t =

∑

w x−

λ

t K 參照圖(3-1)，可以進一步推導出 ( ) ( ) 1 1 ( , ) p m i i i i i i p U x t αK βK = = + =

_∑

+

_∑

(3-14) 除此之外，還可決定出U x t( , )中的 jump ∆U： ( ) 1 m i R L i i U U U

α

K = ∆ = − =

_∑

ɶ _(3-15) 其中α β αɶ_i = −_i _i_。在一維線性黎曼問題中，雖然有 exact solution，但在非線性問題裡需利用疊代等方法，這些動作將耗費大量的時間，因此在實際應用上並不廣泛。為了解決此問題，一般皆求解近似黎曼問題〈approximation Riemann problem〉解而不直接求其 exact solution。在求解近似黎曼問題中最被廣泛應用的方法為 Roe 所提出，亦即為 Roe scheme，其內容如下：假設一維尤拉方程式： 0 U F t x ∂ ₊∂ ₌ ∂ ∂ (3-16) 根據 chain rule，可將方程式(3-16)改寫如下：

(31)

19 0 U F U t U x ∂ ₊ ∂ ∂ ₌ ∂ ∂ ∂ 再令 ( )A U F U ∂ = ∂ ，於是方程式(3-16)可以表示成： ( ) 0 U U A U t x ∂ ₊ ∂ ₌ ∂ ∂ (3-17) 其中，A U( )就稱為 Jacobian 矩陣。

而 Roe scheme 將原本的 Jacobian 矩陣A U( )用常數 Jacobian 矩陣〈constant Jacobian

matrix〉A U U( L, R) ɶ _{代替，因此本來的黎曼問題可以改寫成近似黎曼問題：} ( ) 0 U U A U t x ∂ ₊ ∂ ₌ ∂ ɶ ∂ 0 ( , 0) 0 L R U x U x U x <  = <  (3-18) 於是前述方法可以得到(3-18)的近似解。由以上的原理可得知，在近似黎曼問題上，

Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性轉變成線性，但

是初始條件並沒有改變，因此可以得到方程式(3-16)的近似解。為了要求得合理的常數 Jacobian 矩陣，須合乎 Roe 所提出的四項條件： 1. U 與F 之間，存在著線性轉換的關係。 2. 當UR−UL→U ，則A U U( L, R)→A U( ) ɶ _，此處 F A U ∂ = ∂ 。 3. A Uɶ( _L−U_R)= −F_L F_R。 4. 矩陣Aɶ的特徵向量必須線性獨立。這四項條件都是雙曲線方程式所需具備的，這同時也說明了 Roe 所提出的常數 Jacobian 矩陣必須有實數特徵值，其所對應的特徵向量必須線性獨立。除此之外，條件 3.則是為了符合守恆定律（conservation law）與 Rankine-Hugoniot 條件。

線性黎曼問題的解析解，可以直接從(3-13)與(3-15)式得到， 1 2 ( / ) i U x t + 的解可以利用下面的方程式計算： ( ) 1 0 2 ( / ) i i L i i U x t U K λ α + = +

∑

_< ɶ _(3-19)

(32)

20 或 ( ) 1 0 2 ( / ) i i R i i U x t U K λ α + = −

∑

_> ɶ _(3-20) 其中 1 2 i+ 表示網格與網格之間的交界面(face)。而黎曼問題的近似解，則須從解近似黎曼問題著手： ( ) 0 U F U t x ∂ ₊∂ ₌ ∂ ∂ ɶ ，根據(3-18)式可得知Fɶ = AUɶ 為了符合守恆的條件，因此下式必須成立： ( _R) ( _L) ( _R) ( _L) F Uɶ −F Uɶ =F U −F U (3-21) 接著在固定體積的條件下，積分近似解 ₁ 2 (0) i U + ，可得到通量(flux)的數值公式： 1 1 2 2 ( (0)) ( _R) ( _R) i i F F U F U F U + = + − − ɶ ɶ _(3-22) 再從Fɶ = AUɶ 的關係中可進一步求得： 1 1 2 2 (0) ( _R) _R i i F AU F U AU + = + − − ɶ ɶ _(3-23) 再根據(3-19)式與(3-20)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i R i i i i F F U A K F U K λ α λ α + + = −

∑

_> = −

∑

₌ ɶ _ɶ ɶ _ɶ ɶ _(3-24) 或 ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i L i i i i F F U A K F U K λ α λ α − + = +

∑

_> = +

∑

₌ ɶ _ɶ ɶ _ɶ ɶ _(3-25) (3-24)與(3-25)所指的

λ

i − ɶ _與 i

λ

ɶ+_{分別是代表負的特徵值與正的特徵值，接著再利用平均的} 方法將 ₁ 2 i F + 更進一步表示成： ( ) 1 1 2 1 ( ) ( ) 2 m i R L i i i i F F U F U λ αK + ₌   = _ + − _ 

∑

 ɶ ɶ _(3-26) 再由(3-14)式可再次改變 ₁ 2 i F + 的形式如下： 1 2 1 ( ) ( ) 2 R L i F F U F U A U +   = _ + − ∆ɶ _ (3-27) 其中∆ =U U_R−U_L、 Aɶ = Aɶ+−Aɶ− =Kɶ ɶΛ Kɶ−1， 1

0

m λ λ     Λ =        ⋯ ɶ _⋮ _⋱ _⋮ ⋯ 。

(33)

21 接下來需找出 Aɶ 中所需的物理量，必須利用下列方法：現考慮一維等溫尤拉方程式： ( ) 0 t x U +F U = (3-28) 其中 1 2 u U u u ρ ρ     = =      ； 1 2 2 2 f u F f u a p ρ ρ     = = ₊      ，a為聲速方程式(3-28)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示： 2 2 0 1 ( ) 2 F A U a u u U   ∂ = =_ _ − ∂   (3-29) 特徵值：λ₁= −u a，λ₂ = +u a 特徵向量： (1) 1 K u a   =_ _ −  ， (2) 1 K u a   =_ _ +   接著選定 parameter vector Q 1 2 q U Q q _u ρ ρ ρ     =_{ }=     _ _ (3-30) 再將F 與U 利用Q表示： 2 1 1 1 2 1 2 u q U q Q u q q     = = =      (3-31) 1 1 2 2 2 2 2 2 1 f q q F f q a q     = = ₊      (3-32) 為了表示出∆U與∆F 需在定義 averaged vectorQɶ： 1 2 1 1 ( ) 2 2 L R L R L L R R q Q Q Q q _u _u ρ ρ ρ ρ  ₊    =_{ }= + =   +     _ _ ɶ ɶ ɶ (3-33) 再找出Bɶ =B Qɶ( )ɶ 與Cɶ =C Qɶ ɶ( )使得 U B Q ∆ = ∆ɶ _； F C Q ∆ = ∆ɶ _(3-34) 將(3-34)結合可得

(34)

22 1 ( ) F CB− U ∆ = ɶ ɶ ∆ _(3-35) 再根據上述條件 3 求出近似 Jcaobian 矩陣 1 Aɶ=CBɶ ɶ− (3-36) 為了滿足(3-34)，可以求得 1 2 1 2q 0 B q q   =    ɶ ɶ ɶ ɶ ； 2 1 2 2 1 2 2 q q C a q q   =    ɶ ɶ ɶ ɶ ɶ (3-37) 再帶入(3-36)可得 2 2 0 1 2 A a u u   =_ _ −   ɶ ɶ ɶ (3-38)

uɶ為 Roe averaged velocity

L L R R L R u u u ρ ρ ρ ρ + = + ɶ _(3-39) 因此可以用同樣方法得到以下物理量： L L R R L R v v v ρ ρ ρ ρ + = + ɶ _(3-40) L L R R L R w w w ρ ρ ρ ρ + = + ɶ _(3-41) L L R R L R H H H ρ ρ ρ ρ + = + ɶ (3-42) 1/ 2 [( 1)( 1/ 2 )] aɶ= γ − Hɶ − Vɶ _(3-43) 其中uɶ、vɶ、wɶ 分別代表x方向、y方向、z方向的速度。Hɶ 、aɶ則分別為焓和音速。(3-39) ～(3-42)式中的UL以及UR則是利用 MUSCL 法求出。

(35)

23

y

x<0 x=0 x>0 x

(36)

24

3-3、、、、Monotonic Upstream-Centered Scheme for Conservation Laws(MUSCL)：：：：本論文使用的是採用I. Abalakin等[27]中所使用的插分法。其方程式如下： 1/ 2 1/ 2 1/ 2 L L i i i u₊ = +u ∆u₊ (3-44) 1/ 2 1/ 2 1/ 2 R R i i i u₊ = −u ∆u₊ (3-45) 1 / 2 (1 )( 1 ) ( 1) L i i i i i u ₊ β u ₊ u β u u₋ ∆ = − − + − 1 1 2 ( 3 3 ) c i i i i u u u u θ − + + + − + − + 2 1 1 ( 3 3 ) d i i i i u u u u θ − − + + − + − + (3-46) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i i u₊ β u₊ u β u₊ u₊ ∆ = − − + − ( 1 3 3 1 2) c i i i i u u u u θ − + + + − + − + 1 2 3 ( 3 3 ) d i i i i u u u u θ + + + + − + − + (3-47) 其中(3-46)、(3-47)式中的β 、θc_、θd 值可由表(3-1)中查得。代入不同的值可以得到不同的精度。本論文則是使用三階精度，以減少數值計算的消散性。在程式中，高次項的插分法在不連續的情況下，容易使震盪變大，為了降低震盪，本研究在 MUSCL 法插分出來的方程式中加入 minmod limiter，用來確保程式不會發散。因此(3-44)與(3-45)式需改寫如下： 1/ 2 1/ 2 min mod( 1/ 2) L L i i i u₊ = +u ∆u₊ (3-48) 1/ 2 1/ 2 min mod( 1/ 2) R R i i i u₊ = −u ∆u₊ (3-49)

(37)

25 表 3-1 精度係數值 β θc θd Order 1/3 1/3 1/3 1/3 0 -1/6 0 -1/10 0 0 -1/6 -1/15 3 4 4 5

(38)

26 3-4、、、、Preconditioning 法法法法：：：：為了提高程式的應用範圍，在 Navier-Stokes 方程式中加入 preconditioning 法，讓程式不論在高速或低速流內計算可壓縮流，皆可獲得精確的結果。 S U F G H t x y z ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-50) 上式為原始方程式，接著將保守形式（conserved variables）轉變成主要變數形式（primitive variables），其形式如下： S p U F G H M t x y z ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-51) 其中Up =[p u v w T]T，M為轉換矩陣： 0 0 0 0 0 0 0 0 0 1 p T p T p T p p T p T p u u U v v M U w w H u v w H C ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ       ∂ _ _ = = ∂      ₋ ₊    (3-52) 其中 _p p ρ ρ =∂ ∂ ； T T ρ ρ =∂ ∂ 接著將(3-51)式的方程式乘上矩陣K 2 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ( ) 1 u v K w H V u v w    ₋     ₋  =_ _ −   ₋ ₋ ₋ ₋ ₋      (3-53) 再將K與M相乘 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 p T p KM C ρ ρ ρ ρ ρ ρ         =     ₋    (3-54) 將(3-54)式帶入(3-51)式，連續方程式：

(39)

27 ( ) S p p u v w t x y z ρ ρ ρ ρ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-55) 在理想氣體中可將(3-55)再表示成 2 ( ) S p u v w C t x y z γ ∂ ₊∂ρ ₊∂ρ ₊∂ρ ₌ ∂ ∂ ∂ ∂ (3-56) 其中C為聲速由(3-56)式可以看出，在等密度條件下，由於ρ_p為零，(3-55)式將變成 S u v w x y z ρ ρ ρ ∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ (3-57) 上式即為不可壓流的連續方程式。綜上所述，可以得知只要改變(3-54)式中的ρp項，利用當地流場速度（local velocity）的倒數取代，即可轉換系統中的特徵值，藉此改變低速情況下流場的聲速，使聲速與流場速度冪次級數(order)相同，系統不再受到 CFL(Courant-Friedrichs-Lewy Condition)條件的限制，提高程式的計算效率。利用

θ

取代ρ_p項： 2 1 1 ( ) r p U TC θ = − (3-58) max if if if r U u C U u C u C C u C ε ε ε  _× _{< ×}  = × < <  _>  (3-59) 其中

ε

為一極小的值，約等於 5 10− ，其主要是用來防止停滯點（stagnation point）在計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，U_r必須大於流體的當地擴散速度（local diffusion velocity），因此U_r還需加入下列限制：

max( , ) r r U U x ν = ∆ 將

θ

帶入(3-54)式後，可得到一新矩陣Γnc

(40)

28 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T nc p C θ ρ ρ ρ ρ ρ         Γ =     ₋    (3-60) 經過上述推導之後，方程式從(3-51)式轉變如下： ( ) S p nc U F G H K t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-61) 為了讓(3-61)式中的通量項再度轉換成保守形式，在乘上 1 K− 1 (K _nc) Up F G H S t x y z −_Γ ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-62) 根據(3-62)式，定義 1 0 0 0 0 0 0 0 0 0 1 nc p T u u T v K v T w w T H u v w H C T ρ θ ρ θ ρ ρ θ ρ ρ θ ρ ρ θ ρ ρ ρ ρ − −       −      ₋    Γ = Γ =    ₋       ₋  − +       最後方程式簡化成如下形式： S p U F G H t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-63)

其中Γ為 preconditioning 矩陣，U_p為 primitive form[ , , , , ]P u v w T T

由於方程式在時間項經過改變，因此必須重新推導 Roe 所提出的近似黎曼解。在(3-27) 式中，可以觀察到 ₁ 2 i F + 項，是由 1 ( ( ) ( )) 2 F UR +F UL 的中央差分法加上為了解決不連續面

問題的artificial viscosity term 1 2 A∆U

ɶ _{所組成。加入}

preconditioning的方程式只需在

(41)

29 1 -1 1 -1 1 -1 S ( ) S ( ) S ( ) S p p p p p p p U F G H t x y z U F G H t x y z U _U _U _U A B C t x y z U U U U AM BM CM t x y z − − − ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ ∂ _{+ Γ} ∂ ₊∂ ₊∂ _{= Γ} ∂ ∂ ∂ ∂ ∂ _{+ Γ} ∂ ₊ ∂ ₊ ∂ _{= Γ} ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + Γ + + = Γ ∂ ∂ ∂ ∂ (3-64) 其中 p U M U ∂ = ∂

所以 artificial viscosity terms 改寫如下：

1 1 2 1 1 ( ) 2 R L 2 P i F F F − AM U + = + − Γ ∆ (3-65) 其中 1 1 AM KA DA KA − − Γ = × × 黏滯性項；在黏滯性項方面，採用二階中央差分法。由於在尤拉方程式中計算的範圍皆為網格與網格之間的通量項，因此在黏滯性項方面，所需要得到的速度梯度項也必須是網格之間的通量項。下列以三維的 X 方向為例，圖 3-2 為其示意圖。圖 3-2 中各編號所代表的位置分別為： 1 (i, j+1,k)；2 , , ) 2 1 (i+ j k ；3 (i+1, j+1,k)； 4 (i, j,k−1)；5 , , 1) 2 1 (i+ j k− ；6 (i+1, j,k−1； 7 (i, j,k)；8 , , ) 2 1 (i+ j k ；9 (i+1, j,k)； 10 (i, j,k+1)；11 , , 1) 2 1 (i+ j k+ ；12 (i+1, j,k+1)； 13 (i, j−1,k)；14 , 1, ) 2 1 (i+ j− k ；15 (i+1, j−1,k)；其各速度梯度差分分別如下表示： X U U X U X U ∆ − = ∆ ∆ = ∂ ∂ (9) (7) (3-66)

(42)

30 X V V X V X V ∆ − = ∆ ∆ = ∂ ∂ (9) (7) (3-67) X W W X W X W ∆ − = ∆ ∆ = ∂ ∂ (9) (7) (3-68) Y U U Y U Y U ∆ − = ∆ ∆ = ∂ ∂ 2 ) 14 ( ) 2 ( 2 (3-69) 其中 2 ) 1 ( ) 3 ( ) 2 ( U U U = + ； 2 ) 15 ( ) 13 ( ) 14 ( U U U = + 所以 Y U U U U Y U U Y U U Y U ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 15 ( ) 13 ( ) 1 ( ) 3 ( 2 ) 2 ) 15 ( ) 13 ( ( 2 ) 2 ) 1 ( ) 3 ( ( (3-70) 同理 Y V V V V Y V V Y V V Y V ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 15 ( ) 13 ( ) 1 ( ) 3 ( 2 ) 2 ) 15 ( ) 13 ( ( 2 ) 2 ) 1 ( ) 3 ( ( (3-71) Y W W W W Y W W Y W W Y W ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 15 ( ) 13 ( ) 1 ( ) 3 ( 2 ) 2 ) 15 ( ) 13 ( ( 2 ) 2 ) 1 ( ) 3 ( ( (3-72) Z U U Z U Z U ∆ − = ∆ ∆ = ∂ ∂ 2 ) 5 ( ) 11 ( 2 (3-73) 其中 2 ) 12 ( ) 10 ( ) 11 ( U U U = + ； 2 ) 6 ( ) 4 ( ) 5 ( U U U = + 所以 Z U U U U Z U U Z U U z U ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 6 ( ) 4 ( ) 12 ( ) 10 ( 2 ) 2 ) 6 ( ) 4 ( ( 2 ) 2 ) 12 ( ) 10 ( ( (3-74) 同理 Z V V V V Z V V Z V V z V ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 6 ( ) 4 ( ) 12 ( ) 10 ( 2 ) 2 ) 6 ( ) 4 ( ( 2 ) 2 ) 12 ( ) 10 ( ( (3-75) Z W W W W Z W W Z W W z W ∆ + − + = ∆ + − ∆ + = ∂ ∂ 4 ) 6 ( ) 4 ( ) 12 ( ) 10 ( 2 ) 2 ) 6 ( ) 4 ( ( 2 ) 2 ) 12 ( ) 10 ( ( (3-76)

(43)

31

(44)

32

3-5、、、、Dual time stepping :

方程式(3-51)中的 Navier-Stokes 方程式在時間項方面遭到修改，因此利用修改後的方程式來計算暫態結果並不恰當，因此本程式再加入 dual time stepping[29]的模組，不僅讓程式在計算暫態結果方面較準確，更提高程式的效率，縮短計算時間。

首先，先在原始 Navier-Stokes 方程式加入一虛擬時間項，稱為 artificial time term。

方程式改變如下： S U U F G H t x y z τ ∂ ₊∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ ∂ (3-77)

其中τ 即為 artificial time，t為 physical time

U 為 conservative form

(

ρ ρ ρ ρ ρ, u, v, w, e

)

T

接著在 artificial time term 加入 preconditioning method：

p U U F G H S t x y z τ ∂ ∂ ∂ ∂ ∂ Γ + + + + = ∂ ∂ ∂ ∂ ∂ (3-78)

最後對 artificial time term 採一階的有限差分離散，對 physical time term 採二階的後項差分離散， F x ∂ ∂ 、 G y ∂ ∂ 、 H z ∂ ∂ 利用中央插分法可得 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 , , , , , , , , , , , , 2 2 2 2 2 2 3 4 1 1 1 ( ) ( ) ( ) 2 k k k n n p p k k k k k k i j k i j k i j k i j k i j k i j k U U U U U F F G G H H S t x y z τ + + − + + + + + + + − + − + − − − + Γ + + − + − + − = ∆ ∆ ∆ ∆ ∆ (3-79) 接著整理上式，先將其線性化 1 3( ) 4 ( ) ( ) ( ) 2 n n n p p k k k k k k x p p y p p x p p U U M U U U F A U G B U H C U S t δ δ δ τ − ∆ + ∆ − + Γ + + + ∆ + + ∆ + + ∆ = ∆ ∆ (3-80) 其中∆U_p =Uk_p+1−U_pk， p U M U ∂ = ∂ 1 k n p U + =U + ∆M U ， k 1 k p p F + =F + ∆A U ， k k p p F A AM U ∂ = = ∂

(45)

33 再將∆U_p項放置在等號左邊，其餘則在右邊: 1 3 1 1 ( ) 2 k k k k x p y p x p p I M A B C U R t

δ

τ

− − −   + Γ + Γ + + ∆ = Γ _∆ _∆    (3-81) 此處 1 3 4 ( ) ( ) 2 k n n k k k k x y x U U U R S F G H t δ δ δ − − + = − − + + ∆ ，I為單位矩陣，其中 p U M U ∂ = ∂ ， p p F A U ∂ = ∂ 、 p p G B U ∂ = ∂ 與 p p H C U ∂ = ∂ 為 flux Jacobian。

k為 artificial time 中的疊帶次數， n 為 physical time 的計算階數。上述方程式，當

artificial time term 收斂時， Upk1 Ukp ₀ τ + ₋

Γ =

∆ ，方程式即會回復到原始的 Navier-Stokes 方程

高浮慣比之混合對流研究

國 立 交 通 大 學

機械工程學系

碩士論文