利用CUDA平行計算平台探討可壓縮流在三維煙囪管道的熱傳分析

(1)

國立交通大學

機械工程學系

碩士論文

利用 CUDA 平行計算平台探討可壓縮流在

三維煙囪管道的熱傳分析

Analysis of Heat Transfer in a Three Dimensional

Chimney with Consideration of the Flow

Compressibility by Paralleling the Program

in CUDA Platform

研究生：黃耘

指導教授：傅武雄博士

(2)

利用 CUDA 平行計算平台探討可壓縮流

在三維煙囪管道的熱傳分析

研究生：黃耘指導教授：傅武雄國立交通大學機械工程學系摘要本研究利用數值方法分析可壓縮流在三維煙囪漸縮管道中的流動及熱傳機制。流場利用有限差分法進行計算，計算方法可分為兩部分：第一部份為非黏滯性項的尤拉方程式採用 Roe 方法計算通量，並且加入 Preconditioning 矩陣，讓程式在計算低速可壓縮流可獲得良好之收斂結果；第二部份為黏滯性項的計算，採用二階中央插分法。在時間項方面則採用 LUSGS 隱式法。出口設非反射性邊界條件避免可壓縮流中壓力波的干擾。由數值計算結果得知，流體速度隨三維漸縮管道的主流方向截面積縮小而加快，使壁面中央區域熱傳效果大幅增加。但壁面與壁面間的夾角隨之減小，導致流經夾角附近 (壁面兩側)的流體，因壁面摩擦阻抗增加流速降低，反而使熱傳效果急遽劣化。此外雷諾數越大，熱傳效果越好，壁面兩側三維效應現象依舊存在。而利用 CUDA 平台作平行化運算可使計算速度較一般四核心中央處理器加快約 4.72 倍，可見其顯著的效率提升。

(3)

Analysis of heat transfer in three dimensional chimney

with consideration of the flow compressibility by paralleling

the program in CUDA platform

Student: Yun Huang Advisor:Wu-Shung Fu

Department of Mechanical Engineering National Chao Tung University

Abstract

An investigation of heat transfer in a three-dimensional chimney with consideration of the flow compressibility is studied numerically.The finite difference is adopted and the computational approaches could be divided into two parts. One is the inviscid terms The Roe scheme is utilized to deal with the flux of inviscid terms and the preconditioning matrix is added to let the scheme to be more effective for all speed filed. The other is viscous term and the central difference of second order is utilized to handle it. The temporal term is solved by LUSGS. Non-reflection conditions at the outlet is derived in order to resolve reflections induced by acoustic waves.

The result shows that the convergent chimney construction can accelerate flow velocities by changing the area of cross section ,and the superior enhancement of heat transfer rates of the center of the heat wall are achieved. With the decreasing of the angle between the walls, velocities of the fluids sides by wall decrease due to the friction and the heat transfer efficiency decrease. Moreover, the heat transfer efficiency increase with increasing Reynolds number. In order to improve the efficiency of CFD program,CUDA platform is employed for paralleled program in graphic hardware and shows significantly advantage over the CPU implementation.

(4)

致謝兩年來承蒙指導老師傅武雄教授的悉心指導與諄諄教誨，使本文方能順利完成，謹由衷致上最誠摯的敬意及感謝。此外特別感謝博士班學長李崇綱、連信宏、黃玠超、王威翔在研究上耐心的教導，同窗元盈、為陽、泳良的相互扶持與努力，還有冠藍、佑璁、上豪、仕坤在實驗室帶給我的歡笑及陪伴。最後要感謝父母及家人精神上的支持與照顧，讓我在研究上能心無旁騖。最後僅將此喜悅與求學過程中所有關心與幫助過的親友分享。

(5)

目錄中文摘要 ... i 英文摘要 ... ii 致謝 ... iii 目錄 ... iv 表目錄 ... v 圖目錄 ... vi 符號表 ... viii 一、緒論 ... 1 二、物理模式 ... 9 2-1 物理尺寸與分析模式 ... 9 2-2 分析假設及統御方程式 ... 10 2-3 邊界條件 ... 11 三、數值計算模式 ... 14 3-1 統御方程式 ... 15 3-2 Roe scheme ... 16 3-3 MUSCL 法 ... 23 3-4 Preconditioning ... 25 3-5 LUSGS 法 ... 31 3-6 非反射性邊界 ... 33 3-7 座標轉換 ... 35 3-8 CUDA 平行化方法 ... 42 四、結果與討論 ... 46 五、結論 ... 96 參考文獻 ... 97

(6)

表目錄

表 3-1 精度係數值 ... 24 表 4-1 平行化方法計算效率比較表 ... 58

(7)

圖目錄圖 1-1 GPU 及 CPU 浮點運算速度比較圖 ... 7 圖 1-2 顯示卡 NVIDIA Tesla C1060 安裝於本研究室之電腦示意圖 ... 8 圖 2-1 煙囪管道之物理模式圖 ... 12 圖 2-2(a~c) XY 及 YZ 平面物理模式圖 ... 13 圖 3-1 黎曼問題特徵值結構圖 ... 22 圖 3-2 差分示意圖 ... 30 圖 3-3 座標轉換示意圖 ... 39 圖 3-4(a~c) 漸縮管道示意圖 ... 40 圖 3-5 CUDA 平台架構圖 ... 45 圖 4-1 網格測試圖 ... 54 圖 4-2 完全發展流出口解析解對照圖 ... 55 圖 4-3(a~b) 完全發展流與非反射性邊界出口速度確認圖 ... 56 圖 4-4 使用 CUDA 平行化與未使用平行化紐塞數比較圖 ... 57 圖 4-5(a~b) 中央 XY 截面流線圖及等溫線圖(Re=200，φ =0) ... 59 圖 4-6(a~d) YZ 截面速度場圖(Re=200，φ =0) ... 60 圖 4-7(a~b) 局部紐塞數比較圖(Re=200，φ =0) ... 62 圖 4-8(a~b) 中央 XY 截面流線圖及等溫線圖(Re=200， 6 π φ= ) ... 63 圖 4-9(a~d) YZ 截面速度場圖(Re=200， 6 π φ= ) ... 64 圖 4-10(a~b) 局部紐塞數比較圖(Re=200， 6 π φ= ) ... 66 圖 4-11(a~b) 中央 XY 截面流線圖及等溫線圖(Re=200， 3 π φ= ) ... 67 圖 4-12(a~d) YZ 截面速度場圖(Re=200， 3 π φ= ) ... 68 圖 4-13(a~b) 局部紐塞數比較圖(Re=200， 3 π φ= ) ... 70

(8)

圖 4-14(a~b) 中央 XY 截面流線圖及等溫線圖(Re=200， 2 π φ= ) ... 71 圖 4-15(a~c) YZ 截面速度場圖(Re=200， 2 π φ= ) ... 72 圖 4-16(a~b) 局部紐塞數比較圖(Re=200， 2 π φ= ) ... 74 圖 4-17(a~b) 中央 XY 截面流線圖及等溫線圖(Re=400，φ=0) ... 75 圖 4-18(a~d) YZ 截面速度場圖(Re=400，φ =0) ... 76 圖 4-19(a~b) 局部紐塞數比較圖(Re=400，φ =0) ... 78 圖 4-20(a~b) 中央 XY 截面流線圖及等溫線圖(Re=400， 6 π φ= ) ... 79 圖 4-21(a~d) YZ 截面速度場圖(Re=400， 6 π φ= ) ... 80 圖 4-22(a~b) 局部紐塞數比較圖(Re=400， 6 π φ= ) ... 82 圖 4-23(a~b) 中央 XY 截面流線圖及等溫線圖(Re=400， 3 π φ= ) ... 83 圖 4-24(a~d) YZ 截面速度場圖(Re=400， 3 π φ= ) ... 84 圖 4-25(a~b) 局部紐塞數比較圖(Re=400， 3 π φ= ) ... 86 圖 4-26(a~b) 中央 XY 截面流線圖及等溫線圖(Re=400， 2 π φ= ) ... 87 圖 4-27(a~c) YZ 截面速度場圖(Re=400， 2 π φ= ) ... 88 圖 4-28(a~b) 局部紐塞數比較圖(Re=400， 2 π φ= ) ... 90 圖 4-29(a~b) 各雷諾數平均紐塞數比較圖(Re=200，400) ... 91 圖 4-30(a~d) 各角度平均紐塞數比較圖(Re=200) ... 92 圖 4-31(a~d) 各角度平均紐塞數比較圖(Re=400) ... 94

(9)

Nomenclature

x Cartesian coordinate system x direction

y _{Cartesian coordinate system y direction}

z Cartesian coordinate system z direction

ξ Curvilinear coordinate system ξ direction η _{Curvilinear coordinate system}_{η direction} ς _{Curvilinear coordinate system} _{ς direction}

1

l length of the front segment of the chimney [m] 2

l length of the convergent segment of the chimney [m] 3

l length of the back segment of the chimney [m]

1

d length of the square at inlet [m] 2

d length of the square at outlet [m]

φ convergent angle of the chimney ρ _density[ 3_] kg m⋅ − P pressure[N m⋅ −2] mj τ stress[N m⋅ −2] p

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

v

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

J transform matrix Jacobian

e internal energy[J kg⋅ −1]

T Kelvin temperature[K]

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

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

(10)

β _{pressure gradient[} 3 N m⋅ − ] A area[m2] t Δ time difference[ ] s a _{sound speed[} 1_] m s⋅ − R R H enthalpy[J kg⋅ −1⋅k−1] w

τ friction force acting per unit area on the surface

u_τ friction velocity

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

average Nusselt number difined in Eq.(4-8)

Greek symbols

v kinematics viscosity[m s⋅ −1] μ absolute viscosity[kg.m−1.s−1]

subscripts

i i= x,x-direction; i= y,y-direction; i= z, -direction z

(11)

第一章緒論

煙囪設計在國內外建築中十分常見，是房屋建築的一部分，作為用來排除由火引起的氣體或煙塵，能把煙氣排入高空的高聳結構，藉此減少聚積在下方鍋爐部分之煙氣，使一旁乾淨空氣得以進入，因而能改善燃燒條件，減輕煙氣對環境的污染，煙囪的設計原理主要是依靠自然對流，意即藉由溫度差異而造成空氣的對流而無須額外動力即能有效排氣，因而在生活中運用廣泛。近幾年來，隨著生活結構與建築型態的改變，住家的煙囪數量增加雖趨漸緩，但此種便利且少耗能的對流設計方式依然保留下來，並搭配太陽能板等綠色能源大量應用於大樓內的通風設計，或於建築物頂端開口通風排煙,可藉「煙囪效應」直接將熱氣、煙霧及火流向上排解出去,有助於侷限火勢用以於意外火災時濃煙能夠加速排出；除此之外，石化工業目前仍是台灣經濟命脈中最重要的一環，其相關產業涵蓋了化學肥料、農藥、清潔劑、成衣服飾、塑膠、橡膠、油漆等，產品種類繁多，尤其近幾年來石化相關產品更被延伸使用於電子材料與航太材料等高科技產品上，對於推動我國工業的轉型與升級，產生莫大的助益和貢獻，然而在其相關工廠如中油、台塑之煉油廠、煉鋼廠、各類化工廠，甚至於火力發電廠和垃圾焚化爐卻時常需在大量燃燒後排放高溫氣體以及汙染懸浮物，高聳的煙囪已然成為城市生活中最不可或缺的建築物之一。煙囪管道整體包含入口較大管徑、漸縮部分、以及出口小管徑之管道，其中能夠使煙囪發揮其排煙功效的重要設計，一為底部鍋爐的燃燒加熱，此時靠近最下方之壁面隨之升溫，空氣在此處密度降低，高熱氣體自然向上流動；另一則是漸縮管道藉由截面積的快速減小，可以加速流體的向上流動，並且在此增強壁面對空氣的熱傳，使自然對流能夠充分發揮效果，此一以漸縮或反向之漸擴管道設計增強熱傳效應或自然對流之設計應用廣泛，包含常見的熱交換器、火箭或飛機噴射引擎推進室、太陽能集熱器、特朗博牆(Thrombe walls)[1]、電子元件與薄膜散熱表面等；過去也有許多文獻提供像這樣與煙囪相似的漸縮或漸擴結構的管道模擬或實驗，可用以參考並進一步了解其流場及熱傳狀況

(12)

討論純自然對流部分，Sparrow 等[2]是第一份完整的探討傾斜漸縮平板中的自然對流與熱傳現象之實驗與數值模擬之研究，等溫絕熱之對稱壁面邊界條件下以水為流體，漸縮角度為 0 度、5 度、10 度和 15 度，得傾斜角度有助於熱傳效應率的增加，其實驗與模擬的數值亦可見類似趨勢。此研究被 Sparrow 等[3]延伸探討得更為完整，他們以最大間距估計(Maximum interwall spacing )計算垂直平行壁面之管道的紐塞數，並找出紐塞數與 Rayleigh number、截長比(Aspect ratio)之相關性方程式，且可以轉換運用於漸縮及漸擴管道。Su 和 Lin[4][5]使用 SIMPLEC 之數值方法來計算兩面平面、兩面斜板之流場，在長方形截面之水平漸縮管道中，對稱壁面以等溫絕熱條件，改變漸縮角度為 0 度、5 度和 10 度，以不同方向之進出口設定，形成漸縮及漸擴管道來觀察紐塞數及壓力變化，得漸縮管道在此範圍內當角度增加可增加熱傳並使壓降增高，但漸擴管則壓降漸低。此研究 Su 和 Lin[6]加以實驗測量方式探討長方型截面漸縮或漸擴管道的熱傳和壓降情況，截長比固定為 0.1，當漸縮角度固定時，熱傳效果隨雷諾數增加會使熱傳效果增加而壓降減少，漸擴管中亦有類似效果；而角度的增加結果則與模擬類似。Kihm 等[7]在等溫絕熱壁面的漸縮管道中探討空氣的自然對流現象，實驗方式是使用雷射散斑圖 (Laser Specklegram technique )觀察壁面溫度分佈，實驗以改變漸縮角度 0 度、15 度、30 度、45 度、60 度和八個不同縱橫比(Aspect ratio)來計算其對熱傳效果的影響，其實驗結果與理論分析趨勢接近。Mutama 和 Iacovides[8]以實驗探討較長的微漸縮(約 0.53 度)管道，觀察此二維結構的流線加速度及其速度和溫度邊界層狀態。Said [9]以數值計算模擬壁面等溫絕熱的條件下空氣為工作流體之二維漸縮管道的自然對流現象，漸縮角度為 0 度、2 度、5 度和 10 度，並計算紐塞數在低 Rayleigh number(小於 102 )時和變化參數之最佳相關性。此研究 Shalash[10]以與其相似之物理模式同時以數值和實驗延伸探討，實驗部分以 Mach-Zehnder 干涉儀搭配數位影像系統(Digital Image System)執行，數值部分則延伸計算管道的上游和下游的流動狀況，得在低 Rayleigh number 下若增加漸縮角度可使紐塞數增高，特別在管道出口處此一現象十分明顯；而在高 Rayleigh number 時，增加漸縮角度反而使紐塞數降低，此一現象則在入口處特別突出。Bianco 和 Nardini[11] 以 Fluent 模擬二維結構下兩對稱漸縮平板中，當壁面有固定厚度及熱傳導係數並以固定

(13)

熱通量加熱，出入口皆以延伸計算範圍之較大邊界設定來探討空氣的自然對流現象，此研究得由於在出口平板間隔很小時易出現阻流(chocked flow)現象，導致在低 Rayleigh number 時，漸縮角度對溫度分佈變化劇烈，然在較大時則反之。此研究在 Bainco 等[12] 延伸探討不同的結構設計和熱傳參數，並製成詳細之比較圖，可應用於漸縮角度 0 到 10 度，出入口寬度比小於 1，Rayleigh number 小於之範圍內。 Bainc 等[13]延續此研究，使用固定厚度之鋼板製成兩傾斜板及煙霧加以雷射光頁(Laser sheet)的流場可視化技術，亦計算紐塞數與 Rayleigh number 及漸縮角度之相關性，可應用在角度 10 度， Rayleigh number 在內。Kaiser 等[14]以 Fluent 和 Phoenics 作二維漸縮管道之模擬，在對稱壁面等溫之邊界條件探討自然對流及其熱傳效率，並比對前人文獻之實驗數據，計算紐塞數及 Rayleigh number、漸縮角、出口比例等參數之相關性，可應用於範圍低 Rayleigh number 時(小於 _{)，傾斜角度 0 度到 30 度，漸縮角 0 度到 60 度之間。} 7 1.2 10× 5 1.22 10× 106 Tao 等[15]以週期性三維流場模擬計算探討類似鰭管式熱交換器的波浪型板熱傳情況，以強制對流入口，雷諾數為 500~5000，在不降低效率的前提下為改善網格的精度而使用適體坐標系(Body-Fitted Coordinates, BFC)，其熱傳效率隨雷諾數與波數的增加而上升，而散熱片間距可找出一個最佳化參數而使熱傳效果最好。Castelloes 等[16]則探討層流中流體流經二維波浪型或平滑壁面之熱傳效率，其管道壁面用以不連續的邊界條件，並考慮入口及出口絕熱區域間界面，以廣義積分轉換(Generalized Integral Transform Technique,GITT)解此部分波浪變形的物理模式及暫態方程式，結果顯示熱傳效率隨貝克勒數(Peclet number)降低及波紋狀壁面而增加。Wang 和 Chen[17]模擬牛頓流體在二維波浪形渠道的強制對流熱傳分析，經由代數轉換、座標拉伸推導出流線函數與渦漩度的橢圓形方程式後，配合交換方向的三次樣線定置法(Cubic Spline Collocation Method)，分析其流場與熱傳效應，結果顯示在高雷諾數下，具大波長比及波紋狀壁面的管道為良好的熱傳裝置。Wang 和 Chen[18]延續此研究，使用座標轉換和交替方向隱式法

(Alternating-direction implicit method)解在等溫絕熱壁面之波浪板中的二維流場，並以此探討類似漸縮加漸擴管道中對加強熱傳效應之混合對流影響，得當波長比和雷諾數的增

(14)

加皆能有效的增強熱傳效果。於混合對流方面，亦有 Haung 等[19]研究垂直二維漸縮管道的混合對流實驗，其對稱壁面的一邊以固定熱通量加熱，另一則為絕熱條件，漸縮角度固定為 3 度，以流場可視化之技術顯現流場結構，並探討紐塞數變化及其與其他參數之相關性。Gau 等[20]以類似物理模式而使用樹脂玻璃管道作底部加熱改為水平的平行和對稱漸縮平行管道，流體為空氣並以添加煙霧作流場可視化，用以探討此混合對流的二次流結構及其熱傳的增高情況，在可視化流場中可以清楚的看到渦流因側向延伸而由較圓的渦旋轉變形成縱向對流捲狀，並且其捲狀物數量較水為流體時少。Liu 和 Gau[21]延伸本篇研究，在管道內的流體結構藉由注入煙霧可見在壁面底部流動，觀察漸縮、漸擴、平行板此三種流道之二次流的結構，此渦流肇始於由底部的加熱層造成橫向的不穩定波動，實驗結果可清晰觀察到不同於平行板，在漸擴道中由於減速的擾動影響使二次流渦旋較早出現，並明顯增強熱傳效果，漸縮管道則反之。此外，亦有討論紊流入口之流場，Yang 等[22]以單邊傾斜成漸縮管道之平板實驗，以高紊流強度之入口條件模擬薄膜冷卻的效果，比較於不同紊流強度、漸縮角度 0 度到 20 度之壁面的溫度分佈及其熱傳效率，研究得此小角度傾斜下冷卻效率隨角度增加，但角度太大則否，且紊流強度增加有助於熱傳效果。Zamora 等[23]延續 Kaiser 等[14]之研究，入口加以不同紊流強度，使用低雷諾數下的k−ε 紊流模式，探討過渡層到紊流層的範圍內之流場的質量流率及變動參數之相關性，可應用於高 Rayleigh number，約在，傾斜角 0 度到 30 度，漸縮角 0 度到 60 度之間。 2 10 ~ 10− 12 綜觀以上的文獻，為了簡化管道流的研究，除了文獻[4][15]外，都以二維模式為研究對象。為了配合 Bossinesq assumption 的使用，流體均視為不可壓縮流，因此加熱物與流體間的溫度差必須限於 30°C 以內[24]。然而真實煙囪的形狀為一三維漸縮管道，雖然內部流體之流速未達壓縮流流場之 0.3 馬赫(Mach number)的流速，但加熱壁爐與流體間的溫差遠高於 30°C，因此內部流場居於低速可壓縮流狀態，Bossinesq assumption 無法成立。流體必須視為可壓縮性流體，方可符合真實物理現象。

(15)

基於以上原因，本文利用數值方法探討可壓縮流在三維漸縮煙囪管道內流動與熱傳機構。在計算具進出口條件之管道內低速可壓縮流場時，最棘手的問題是流體的流動速度和壓力波的速度相差過大，導致在出口反射回彈的壓力波會干擾流體流動，因此必須在出口處設置非反射性(non-reflecting)邊界條件[25]，避免干擾發生。其次採用慣用的中央插分法處理黏性項部分；而非黏性項則必須採用 Roe 法[26]，以解決偏微分方程式的不連續問題；為了處理在低速可壓縮流時不易收斂的困擾，Preconditioning 法[27]的加入，可有效提高收斂效率，使其在高速與低速可壓縮流場均可適用；再利用

MUSCL(Monotone Upstream-centered Schemes for Conservation Laws)法[28]計算網格與網格間的物理量；最後使用 LUSGS(implicit lower-upper symmetric Gauss-Seidel algorithm) 法[29]做迭代運算。通常計算不規則形狀熱流場時，採用下列三種方法。一為非結構性網格法(Unstructured Grid )，網格排列依不規則形狀而編排，不易有秩序排列，較不利於使用後述本研究所開發的平行運算法；二為偏微分方程式格點產生法(partial differential equation grid generation)[15]，程式撰寫上較複雜且生成速度較慢；三為代數格點產生法 (algebraic grid generation )[16][17]，程式撰寫較簡單，生成速度較快，並且容易控制格點疏密排列，較為適合後述本研究所開發的平行運算法。因此本研究採用代數格點產生法 [30]將三維漸縮煙囪管道轉換成方形管道，加速計算效率。眾所周知在三維流場的計算極為複雜耗時，電腦必須有龐大的儲存容量與有效的計算方法方以應付。平行運算為一種有效的方法，目前最常見可分為兩大類，一為 OpenMP(Open Multi-Processing)，一台電腦擁有多核心(multicores)的中央處理器(Central Processing Unit，CPU)，目前最多為十二核心的中央處理器，方便節省空間與金錢，但共同執行的處理器核心有限。另一為 MPI(The Message Passing Interface,)利用網路連結各電腦，平行運算界網路連結之電腦群同時執行，可連結極多數，但效率會逐漸降低，且價格昂貴須有廣大空間供電腦置放。為解決上述缺失，本研究自力開發現現有電腦與英偉達(NVIDIA)公司的硬體顯示卡(價格約四萬台幣)之間的連結與程式編輯，將平行運算於 CUDA((Compute Unified Device Architecture)計算平台上執行。亦即將顯示卡內原為繪圖功能而具有利於浮點運算、數量極多執行緒(multithread)與高記憶體頻寬(memory bandwidth)的特性，轉換成所需的計

(16)

算功能。過去以顯示卡(NVIDIA 8800GTX )作為平行運算工具，在 CUDA 平台上作數值模擬的有 Brandvik 和 Pullan[31]，以有限體積法離散三維尤拉方程式，解渦輪葉片上的非黏性可壓縮流，得其計算效能約為中央處理器(Intel Core 2 Duo 2.33 GHz)的 16 倍，測試網格較少時 CPU 計算速度則會相較有些微增加；Corrigan 等[32]用非結構性網格，以 NACA0012 機翼及一飛彈為例，解此超音速可壓縮流流場，比較中央處理器(Intel Core 2 Q9450)及顯示卡(NVIDIA Tesla)在不同案例下單精準及雙精準的計算效能，結果在單精準浮點數下使用顯示卡計算速度約為四核心中央處理器的 9.4~9.9 倍，而雙精準則為 1.56~2.5 倍。圖(1-1)為英特爾(Intel)公司的中央處理器與英偉達(NVIDIA)公司的硬體顯示卡單浮點運算速度的比較圖[33]，可得知兩者的差距。圖(1-2)則是英偉達公司的硬體顯示卡安裝於本研究室之電腦示意圖，利用 CUDA 計算平台於平行運算，不但可大幅縮短計算時間，且可降低電腦成本與使用空間。利用上述計算方法與工具執行數值計算，結果得知，三維漸縮管道隨主流方向截面積逐漸縮小，流速加快，導致壁面中央區域熱傳效果大幅增加。但壁面與壁面間的夾角隨之減小，導致流經夾角附近(壁面兩側)的流體，因壁面摩擦阻抗增加流速降低，而使熱傳效果急遽劣化，三維效應的此現象過去從未發現。此外探討雷諾數的影響，雷諾數越大，熱傳效果越好，壁面兩側三維效應現象依舊存在，而利用 CUDA 平台作平行化運算可使計算速度較一般四核心中央處理器加快四倍以上，可見其顯著的效率提升。

(17)

(18)

NVIDIA Tesla C1060

(19)

第二章物理模式

2-1、物理尺寸與分析模式： 圖(2-1)為一正方形截面之煙囪管道圖，其詳細長度及溫度邊界條件如圖(2-2)所示，前段管道長度為，漸縮管道長度為，後段管道長度為；漸縮角度為l₁ l2 l3 φ，其中入口大正方形截面之邊長為，出口小正方形截面之邊長為。前段管道與漸縮管道為高溫壁面，其他部分則為絕熱壁面。初始溫度的冷卻空氣以均勻流等速度由管道入口噴入，高溫壁面溫度為。其中 1 d h T 2 d 0 T u₀ x方向為 streamwise 方向，y為垂直壁面方向，最後為 spanwise 方向。u、、分別為其所對應的速度。為簡化管道流研究，以往大多將出口設為速度完全發展流條件及壓力為大氣壓力邊界，但此方式在可壓縮流中壓力波在出口處易會反射回計算區域而影響收斂性，較不適用於此低速可壓縮流流場；且以往為使熱管道內的熱流場更容易成為完全發展流，在管道的長度以及寬度上需要一定的比例，但此一方式會使網格需求增加，計算時間較長，因此本研究在出口條件部分設為非反射性(non-reflecting)邊界[25]條件，此時可視其出口為遠處的邊界條件，選定後段管道長度約為最小寬度的兩倍長，大幅減少網格數量。以往研究中煙囪出口大小因受漸縮角度影響，出口太小時流體速度較快，此時出口處阻力較大，可能因此使煙囪管道中的煙灰無法順利排出；若出口太大，則流體速度降低，然因此可能造成管道壁面之散熱面積加大，而造成流體降溫，排煙效果不良，此外因流體煙氣在尾端停留過久亦較有機會對管道造成腐蝕，故而實用上漸縮角度多在 45 度以內。本文擴大研究範圍，漸縮的角度則選定為 0 度、30 度、60 度及 90 度。 z v w

(20)

2-2、分析假設及統御方程式： 本研究選擇層流流場作為模擬流場，設定如下： 1. 可壓縮流，空氣密度會隨溫度與壓力而改變。 2. 工作流體為空氣，假設為理想氣體。流體性質為牛頓流體(Newtonian fluid)，黏滯係數為等方向性。 3. 忽略重力效應影響。

4. 所有壁面均為不可滑移(No slip condition)

統御方程式分別為連續方程式(Continuity equation)、動量方程式(Momentum equations)、能量方程式(Energy equations)，與理想氣體狀態(Equation of State of Ideal Gas)方程式。

連續方程式： ( _j) j u t x 0 ρ _ρ ∂ ₊ ∂ ∂ ∂ = (2-1) 動量方程式： ( ) i j i ij j i u p u u t x x xj ρ _ρ _σ ∂ ₊ ∂ _{= −}∂ ₊ ∂ ∂ ∂ ∂ ∂ 其中 [ 2 ( 3 j i k ij jf ij j i k u u u )] x x x σ =μ_∈ ∂ +∂ − δ ∂ ∂ ∂ ∂ (2-2) 能量方程式： ( _j _j) _j ( _i j j E u E pu q u t x x x_j ij) ρ _ρ _σ ∂ ₊ ∂ ₊ _{= −} ∂ ₊ ∂ ∂ ∂ ∂ ∂ (2-3) 理想氣體狀態方程式： p=ρRT (2-4)

(21)

2-3、邊界條件：

本研究所採用的統御方程式為可壓縮 Navier-stokes 方程式，因此需要給定的邊界條件有：初始狀態(Initial condition)、入口條件(Inlet condition)、出口條件(Outlet condition)、壁面邊界(Boundary condition)。 2-3.1 初始狀態：初始速度、初始壓力、初始密度初始速度 u： 0 /m s 初始速度 v： 0 /m s 初始速度 w： 0 /m s 初始壓力 p：一大氣壓力(101300Pa) 初始密度ρ ：空氣密度(1.28kg m/ 3) 2-3.2 入口條件：入口速度 u：均勻流u₀ 入口速度 v： 0 /m s 入口速度 w： 0 /m s 2-3.3 出口條件：出口速度 u：非反射性邊界條件出口溫度 T：非反射性邊界條件出口壓力 p：非反射性邊界條件 2-3.4 壁面邊界：邊界速度：不可滑移條件，u=v=w=0m/s 邊界溫度：高溫壁面溫度T_h 邊界密度：在垂直壁面方向，梯度為零， 0 n ρ ∂ = ∂ 邊界壓力：壓力設為 P 0 n ∂ ₌ ∂

(22)

x

y z

(23)

2 d 1

d

z y φ 2 l x y 圖 2-2(a) XY 平面漸縮部分截面圖圖 2-2(b) YZ 平面物理模式圖

0 T

y

∂

₌

∂

h

T

=

T

u

0

T

3

l

2

l

1

l

x y 圖 2-2(c) XY 平面物理模式圖

(24)

第三章數值計算模式

本章主旨在說明本論文的數值計算所使用的所有模式。第一節首先要知道所求解的方程式為 Navier-Stokes 方程式。參考 Fu 和 Li 等人[34]之數值方法，本研究將方程式拆解為非黏滯項與黏滯項。其中在計算黏滯性項則採用二階精度的中央差分法。第二節介紹的為黎曼解中的 ROE 法(Roe scheme)[26]，利用 ROE 法來求出非黏滯項的通量。接著第三節介紹 MUSCL 法(Monotone Upstream-centered Schemes for Conservation

Laws)[28]，此法是為了要解出 ROE 法中使用的網格之間的物理量，然後為了防止在高階插分時產生震盪現象，在 MUSCL 法插分的結果方程式中加入 Minmod limiter 以確保程式不會發散。第四節為介紹 Preconditioning 法，因為當計算低速可壓縮流時，因速度和音速的數量及上差距過大，在數值分析時不好計算，所以為彌補此一缺點須使用 Preconditioning 法。而第五節為 LUSGS 法(implicit lower-upper symmetric Gauss-Seidel algorithm)[29]，程式因為在使用 Preconditioning 法，故而在統御方程式中修改了計算時階，故需要加入 LUSGS 法，待時階項收斂時才能達到真實穩態。第六節則詳細介紹本文在出口條件部分採用的非反射性邊界條件，用以避免低速可壓縮流中，出口處壓力波易反彈而造成的干擾；第七節為座標的轉換，為了要觀察三維漸縮煙囪管道的流場及熱場的變化及不規則邊界的影響，因此對於此三維模擬建構系統之邊界採用代數格點產生法將座標轉換的方法，使物理平面之漸縮處拉伸展開成一規則的計算平面，因此當程式作計算時使用曲線座標系做計算，當差分時於直角座標系做差分的動作，使格點排列方便後述的平行計算方法，對於方程式座標的轉換在此做說明。第八節為本研究所使用的硬體顯示卡平行化方法，在 CUDA 平台上做平行化運算可大幅提升運算效率，在此一小節詳細的說明。綜合上述，本文的計算過程為利用 MUSCL 法算出 ROE 法所需的網格間物理量，求出非黏滯的通量，並在計算通量時加入 Preconditioning 法，以拉近與音速的 Order。接下來使用 2 階中央插分法對黏滯項做差分進而求出黏滯項；然後再與 ROE 法求出的通量結合得到真正的物理通量。最後使用 LUSGS 法疊代以求出正確時階的物理量。

(25)

3-1、統御方程式： 本研究在計算流場的方面其統御方程式分可為兩大部分，第一部份為非黏滯性項的尤拉方程式，第二部份為黏滯性項。 3 1 2 F ₀ F F U t x y z ∂ ∂ ∂ ∂ ₊ ₊ ₊ ₌ ∂ ∂ ∂ ∂ (3-1) 其中， 1 2 3 ( , , , , )T U= ρ ρ ρ ρ ρu u u e 。 1 1 1 2 2 2 3 3 3 ( ) m m m m m m m m m m m m m m u u u P F u u P u u P T e P u u x ρ ρ δ τ ρ δ τ ρ δ τ ρ λ τ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − ⎟ ⎜ =_⎜ + − ⎜ ₊ ₋ ⎟ ⎜ ⎟ ∂ ⎜ ₊ ₋ ₋ ⎟ ⎜ _∂ ⎟ ⎝ j j⎠ ⎟ ⎟ 。 (3-2) ρ 為密度，P為壓力。、、分別為u1 u2 u3 x、y、z方向的速度。λ ρ= C kp ，k 為 thermal diffusivity。 2 1 2 ( +u2 u + 2 3 1 2 v C T u = + ) C_v e ，為等容比熱。τ_ij為 stress tensor。上式可拆解為黏滯性項與非黏滯性項： 1 1 1 2 2 2 3 3 3 0 ( ) m m m m m inviscid viscid m m m m m m mj j m i u u u P F F F u u P u u P T u e P u x ρ ρ δ τ ρ δ τ ρ δ τ τ ρ λ ⎛ ⎞ ⎜ _{⎟ ⎛} _⎞ ⎜ _{⎟ ⎜} _⎟ ⎜ + ⎟ _⎜ _⎟ ⎜ _{⎟ ⎜} _⎟ = + =_⎜ + _{⎟ ⎜}− ⎟ ⎜ ₊ ⎟ _⎜ _⎟ ⎜ _{⎟ ⎜}_⎜ _⎟_⎟ ∂ ⎜ ₊ ₋ _{⎟ ⎝} _⎠ ⎜ _∂ ⎟ ⎝ ⎠ ∀ =m 1, 2, 3。 (3-3) 左式由非黏滯項F_inviscid 組成的方程式即稱為尤拉方程式。

(26)

3-2、Roe scheme： 在雙曲線的守恆形式方程式中，若其初始條件包含有不連續的片段連續(piecewise) 常數，此類型的問題通稱為黎曼(Riemann)問題。因為其包含有不連續解，因此在流體計 算上有著相當廣泛的應用。一維線性黎曼方程式如下： 0 U U A t x ∂ ∂ + = ∂ ∂ ，其中A為一常數 Jacobian 矩陣。 (3-4) 初始條件為 (0) (0) (0) 。左上方括號代表時間 t 為 0。 1 ( , , _m )T U = u " u 求出A之特徵值矩陣以及特徵向量。 1 A= ΛK K− ，其中Λ為特徵值矩陣： 1

0

m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ " # % # " 。 (1) ( ) , , m T K = ⎣⎡K " K ⎤⎦ 為特徵向量，故AK( )i =λ_iK( )i 。 接著定義特徵變數W 〈characteristic variables〉，其定義如下： ( , ) W =W t x ， ( 1) 或。因此 W =K − U U =KW U K W t t ∂ ∂ = ∂ ∂ 且 U W K x x ∂ ∂ = ∂ ∂ ，將此結果代入 (3-4)式中可得： = 0 t x KW +AKW ，可再繼續簡化成： 0 t x W + ΛW = (3-5)

方程式(3-5)稱為 canonical form 或 characteristic form。將以上的結果簡單整理如下： 0 i i i W W t λ x ∂ ₊ ∂ ∂ ∂ = 2 ₌ ，或 (3-6) 1 1 1 2 3 3 ... 0 0 ... 0 0 : : : : : 0 ... _m t x w w w w w w λ λ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ₊⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (3-6)可由特徵曲線法求得其解為：

(27)

(0) ( , ) ( ) i i i w x t =w x−λt =α_i x−λ_it< 0 (0) ( , ) ( ) i i i w x t =w x−λt = β_i x−λ_it> (3-7) 0 其中，α_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 p U x t α K βK = = + =

∑

+

∑

i i (3-8) 除此之外，還可決定出U x t( , )中的 jump Δ ： U ( ) 1 m i R L i i U U U αK = Δ = − =

∑

(3-9) 其中αi =β αi− i。在一維線性黎曼問題中，雖然有 exact solution，但在非線性問題裡需利用疊代等方法，這些動作將耗費大量的時間，因此在實際應用上並不廣泛。為了解決此問題，一般皆求解近似黎曼問題〈approximation Riemann problem〉解而不直接求其 exact solution。在求解近似黎曼問題中最被廣泛應用的方法為 Roe 所提出，亦即為 Roe scheme，其內容如下：假設一維尤拉方程式： 0 U F t x ∂ ∂ + = ∂ ∂ (3-10) 根據 chain rule，可將方程式(3-10)改寫如下： 0 U F U t U x ∂ ₊ ∂ ∂ ∂ ∂ ∂ = 再令 ( ) F A U U ∂ = ∂ ，於是方程式(3-10)可以表示成： ( ) 0 U U A U t x ∂ ₊ ∂ ∂ ∂ = (3-11) 其中，A U( )就稱為 Jacobian 矩陣。

而 Roe scheme 將原本的 Jacobian 矩陣A U( )用一常數 Jacobian 矩陣〈constant Jacobian matrix〉 (A U U _L, _R)代替，因此本來的黎曼問題可以改寫成近似黎曼問題： ( ) 0 U U A U t x ∂ ₊ ∂ ₌ ∂ ∂

(28)

( , 0)

U x =U L x< (3-12) 0

( , 0)

U x =U R x>0

於是前述方法可以得到(3-10)的近似解。由以上的原理可得知，在近似黎曼問題上， Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性轉變成線性，但是初始條件並沒有改變，因此可以得到方程式(3-10)的近似解。為了要求得合理的常數 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-7)與(3-8)式得到， ₁ 2 ( / ) i U x t + 的解可以利用下面的方程式計算： ( ) 1 0 2 ( / ) i i L i i U x t U K λ α + = +

∑

_< (3-13) 或 ( ) 1 0 2 ( / ) i i R i i U x t U K λ α + = −

∑

_> (3-14) 其中 1 2 i+ 表示網格與網格之間的交界面(face)。而黎曼問題的近似解，則須從解近似黎曼問題著手： ( ) 0 U F U t x ∂ ₊∂ ∂ ∂ = ) ，根據(3-12)式可得知F = AU 為了符合守恆的條件，因此下式必須成立： ( _R) ( _L) ( _R) ( _L F U −F U =F U −F U (3-15)

(29)

接著在固定體積的條件下，積分近似解 ₁ 2 (0) i U + ，可得到通量(flux)的數值公式： 1 1 2 2 ( (0)) ( _R) ( i i F F U F U F U + = + − − ₎ R _(3-16) 再從F = AU的關係中可進一步求得： 1 1 2 2 (0) ( _R) i i F AU F U AU + = + − − R _(3-17) 再根據(3-13)式與(3-14)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) m i i R i R i i i i F F U A K F U K λ α + + = −

∑

_> = −

∑

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

∑

_> = +

∑

₌ i i λ α (3-19) (3-18)與(3-19)所指的λ 與_i− λ 分別是代表負的特徵值與正的特徵值，接著再利用平均的_i+ 方法將 ₁ 2 i F + 更進一步表示成： ( ) 1 1 1 ( ) ( ) 2 m i R L i i i i F F U F U λ αK + ₌ ⎡ = _⎢⎣ + − ⎦

∑

2 ⎤ ⎥ (3-20) 再由(3-8)式可再次改變 ₁ 2 i F + 的形式如下： 1 2 1 ( ) ( ) 2 R L i F F U F U A + = ⎡⎣ + − U ⎤Δ ⎦ (3-21) 其中Δ =U U_R −U_L、 1 A =A+−A− =K Λ K− ， 1

0

m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ " _# _% _# " 。接下來需找出 A中所需的物理量，必須利用下列方法：現考慮一維等溫尤拉方程式： ( ) 0 t x U +F U = (3-22) 其中 1 2 u U u u ρ ρ ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥ ⎢}= _⎥ ⎣ ⎦ ⎣ ⎦ ； 1 2 2 2 f u F f u a p ρ ρ ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥ ⎢}= _⎥ + ⎣ ⎦ ⎣ ⎦ ， a 為聲速

(30)

方程式(3-22)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示： 2 2 0 1 ( ) 2 F A U a u u U ⎡ ∂ = _{= ⎢} − ∂ _⎣ _⎦ ⎤ ⎥ (3-23) 特徵值：λ₁= − ，u a λ₂ = + u a 特徵向量： (1) 1 K u a ⎡ ⎤ = ⎢ ₋ ⎥ ⎣ ⎦， (2) 1 K u a ⎡ ⎤ = ⎢ ₊ ⎥ ⎣ ⎦ 接著選定 parameter vector Q 1 2 q U Q q _u ρ ρ ρ ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥}= ⎢ ⎣ ⎦ ⎢_⎣ ⎥_⎦⎥ ⎥ 1 (3-24) 再將F與U 利用Q表示： 2 1 1 1 2 1 2 u q U q Q u q q ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥}= _{= ⎢} ⎣ ⎦ ⎣ ⎦ (3-25) 1 1 2 2 2 2 2 2 f q q F f q a q ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥ ⎢}= _⎥ + ⎣ ⎦ ⎣ ⎦ (3-26) 為了表示出ΔU與ΔF 需在定義 averaged vector Q ： 1 2 1 1 ( ) 2 2 L R L R L L R q Q Q Q q _u _u ρ ρ ρ ρ ⎡ + ⎤ ⎡ ⎤ =_{⎢ ⎥}= + = ⎢ + ⎢ ⎥ ⎣ ⎦ _⎣ _⎦ R ⎥ (3-27) 再找出B=B Q( ) 與C=C Q ( )使得 U B Q Δ = Δ ； FΔ = ΔC Q (3-28) 將(3-28)結合可得 1 ( ) F CB− U Δ = Δ (3-29) 再根據上述條件 3 求出近似 Jcaobian 矩陣 1 A=CB − (3-30) 為了滿足(3-28)，可以求得

(31)

2 ⎤ ⎥ ⎤ ⎥ 1 2 1 2q 0 B q q ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ ； (3-31) 2 1 2 1 2 2 q q C a q q ⎡ = ⎢ ⎣ ⎦ 再帶入(3-30)可得 2 2 0 1 2 A a u u ⎡ = ⎢ ₋ ⎣ ⎦ (3-32)

u 為 Roe averaged velocity

L L R R L R u u u ρ ρ ρ ρ + = + (3-33) 因此可以用同樣方法得到以下物理量： L L R R L R v v v ρ ρ ρ ρ + = + (3-34) L L R R L R w w w ρ ρ ρ ρ + = + (3-35) L L R R L R H H H ρ ρ ρ ρ + = + _(3-36) 1/ 2 [( 1)( 1/ 2 )] a= γ − H − V (3-37) 其中 u 、、 v w 分別代表x 方向、y方向、方向的速度。z H、則分別為焓和音速。 a (3-33)～(3-36)式中的U_L以及U_R則是利用MUSCL法求出。

(32)

y p λ i λ 2 λ 1 λ λ_m 1 m λ ₋ i λ 1 p λ ₊ x< 0 x= 0 x> x 0 圖 3-1 黎曼問題特徵值結構圖

(33)

3-3、Monotonic Upstream-Centered Scheme for Conservation Laws(MUSCL)： 本論文使用的是採用Abalakin等[28]中所使用的插分法。其方程式如下： 1/ 2 1/ 2 1/ 2 L i i i u₊ = +u ΔuL₊ R + 2) + 1) 2) + 3) (3-38) 1/ 2 1/ 2 1/ 2 R i i i u₊ = −u Δu (3-39) 1 / 2 (1 )( 1 ) ( 1) L i i i i u₊ β u₊ u β u Δ = − − + −u_i₋ 1 1 ( 3 3 c i i i i u u u u θ − + + − + − + 2 1 ( 3 3 d i i i i u u u u θ ₋ ₋ ₊ + − + − + (3-40) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i u₊ β u₊ u β u₊ Δ = − − + −u_i₊ +θc(−u_i₋₁+3u_i−3u_i₊₁+u_i 1 2 ( 3 3 d i i i i u u u u θ ₊ ₊ ₊ + − + − + (3-41) 其中(3-40)、(3-41)式中的β 、θ 、c θ 值可由表(3-1)中查得。代入不同的值可以得到不d 同的精度。本論文則是使用三階精度，以減少數值計算的消散性。在程式中，高次項的插分法在不連續的情況下，容易使震盪變大，為了降低震盪，本研究在 MUSCL 法插分出來的方程式中加入 minmod limiter，用來確保程式不會發散。因此(3-38)與(3-39)式需改寫如下： 1/ 2 1/ 2 min mod( 1/ 2) L i i i u+ = +u ΔuL+ R + (3-42) 1/ 2 1/ 2 min mod( 1/ 2) R i i i u+ = −u Δu (3-43)

(34)

表 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

(35)

3-4、Preconditioning 法： 為了提高程式的應用範圍，在 Navier-Stokes 方程式中加入 preconditioning 法，讓程式不論在高速或低速流內計算可壓縮流，皆可獲得精確的結果。 0 U F G H t x y z ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-44) (3-44)為原始方程式，接著將保守形式（conserved variables）轉變成主要變數形式（primitive variables），其形式如下： 0 p U _F _G _H M t x y z ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-45) 其中U_p =[p u v w T]T，M為轉換矩陣： 0 0 0 0 0 0 0 0 0 1 p T p T p p p T p T u u U v M U w w H u v w H C ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ∂ _⎢ _⎥ = = ∂ ⎢ ⎥ ⎢ ⎥ ⎢ ₋ ₊ ⎥ ⎣ ⎦ T p v ρ (3-46) 其中 _p p ρ ρ =∂ ∂ ； T T ρ ρ = ∂ ∂ 接著將(3-45)式的方程式乘上矩陣K 2 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ( ) u v K w H V u v w ⎡ ⎤ ⎢ ₋ ⎥ ⎢ ⎥ ⎢ − = ⎢ ₋ ⎥ ⎢ ⎥ ⎢₋ ₋ ₋ ₋ ₋ ⎥ ⎢ ⎥ ⎣ ⎦ 0 1 ⎥_(3-47) 再將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-48)

(36)

將(3-48)式帶入(3-45)式，連續方程式： ( ) 0 p p u v w t x y z ρ ρ ρ ρ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-49) 在理想氣體中可將(3-49)再表示成 2( ) 0 p u v w C t x y z γ ∂ ₊∂ρ ₊∂ρ ₊∂ρ ₌ ∂ ∂ ∂ ∂ (3-50) 其中 C 為聲速 由(3-50)式可以看出，在等密度條件下，由於ρ_p為零，(3-49)式將變成 0 u v w x y z ρ ρ ρ ∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ (3-51) 上式即為不可壓流的連續方程式。綜上所述，可以得知只要改變(3-48)式中的ρ_p項，利用當地流場速度（local velocity）的倒數取代，即可轉換系統中的特徵值，藉此改變低速情況下流場的聲速，使聲速與流場速度冪次級數(order)相同，系統不再受到 CFL(Courant-Friedrichs-Lewy Condition)條件的限制，提高程式的計算效率。利用θ取代ρ_p項： 2 1 1 ( r p U TC θ = − ) (3-52) max U ε× if u < × ε C r U = u if ε× <C u < (3-53) C C if u > C 其中ε為一極小的值，約等於 5 10− ，其主要是用來防止停滯點（stagnation point）在計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，必須大於流體的當地擴散速度（local diffusion velocity），因此還需加入下列限制：

r U r U max( , ) r r U U x ν = Δ

(37)

將θ帶入(3-48)式後，可得到一新矩陣Γ nc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T nc p C θ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ Γ = ⎢ ⎥ ⎢ ⎥ ⎢₋ ⎥ ⎣ ⎦ ⎥ (3-54) 經過上述推導之後，方程式從(3-44)式轉變如下： ( p nc U _F _G _H K t x y z ∂ ∂ ∂ ∂ Γ + + + ∂ ∂ ∂ ∂ )=0 (3-55) 為了讓(3-55)式中的通量項再度轉換成保守形式，在乘上 1 K− 1 (K _nc) Up F G H t x y z −_Γ ∂ ₊∂ ₊∂ ₊∂ ∂ ∂ ∂ ∂ =0 (3-56) 根據(3-56)式，定義 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 ρ θ ρ θ ρ ρ θ ρ ρ θ ρ ρ θ ρ ρ ρ ρ − − ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ₋ ⎢ Γ = Γ = ⎢ ⎥ ⎢ ₋ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ − ⎥ − + ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎥ ⎥ 最後方程式簡化成如下形式： 0 p U F G H t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-57)

其中為 preconditioning 矩陣，為 primitive form

由於方程式在時間項經過改變，因此必須重新推導 Roe 所提出的近似黎曼解。在(3-21) 式中，可以觀察到 Γ U_p [ , , , , ] /P u v w T T _J 1 2 + 項，是由 1 ( ( ) ( )) 2 F UR +F UL i F 的中央差分法加上為了解決不連續面

問題的 artificial viscosity term 1

2 A UΔ

_{所組成。加入 preconditioning 的方程式只需在}

(38)

1 1 1 0 ( ) 0 ( ) 0 ( ) 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 − − − ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + Γ + + = ∂ ∂ ∂ ∂ ∂ _{+ Γ} ∂ ₊ ∂ ₊ ∂ ₌ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + Γ + + = ∂ ∂ ∂ ∂ 0 (3-58) 其中 p U M U ∂ = ∂

所以 artificial viscosity terms 改寫如下： 1 1 2 1 1 ( ) 2 R L 2 i F F F − AM + = + − Γ ΔU (3-59) P 其中 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-60) X V V X V X V Δ − = Δ Δ = ∂ ∂ (9) (7) (3-61)

(39)

X W W X W X W Δ − = Δ Δ = ∂ ∂ (9) (7) (3-62) Y U U Y U Y U Δ − = Δ Δ = ∂ ∂ 2 ) 14 ( ) 2 ( 2 (3-63) 其中 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-64) 同理 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-65) 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-66) Z U U Z U Z U Δ − = Δ Δ = ∂ ∂ 2 ) 5 ( ) 11 ( 2 (3-67) 其中 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-68) 同理 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-69) 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-70)

(40)

(41)

3-5、LUSGS 法: 方程式中(3-58)中的 Navier-Stokes 方程式在時間項方面遭到修改，本文選用收斂速度較快的 LUSGS 法來計算修改後的方程式。原方程式: 0 p U F G H t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-71) 首先，對時間項採二階的後項差分離散，可得 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 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 F F G G H H t x y z + − + + + + + + + − + − + − − + Γ + − + − + − Δ Δ Δ Δ )=0 (3-72) 接著整理(3-60)式，先將其線性化 1 3( ) 4 ( ) ( ) ( 2 k k k p k k k k k k x p p y p p x p p U M U U U F A U G B U H C U t δ δ δ − + Δ − + Γ + + Δ + + Δ + Δ + Δ )=0 (3-73) 其中 k1 k p p p U U + U Δ = − ， ₁ ₁ , , , , 2 2 1 ( ) k k k x i j k i j k F F F x δ + − = − Δ ，Ap=AM

( )

1 , , 2 1 , , 2 1 i j k k k k xAp Ap _i _{j k} Ap x δ − + ⎡ ⎤ = ⎢ − ⎥ Δ ⎢_⎣ ⎥_⎦。再將ΔU_p項放置在等號左邊，其餘則在右邊: 1 3 ( ) 2 k k k x p y p x p p k M A B C U t δ δ δ − + Γ + + Δ = Δ R (3-74) 此處 (3 4 1) ( 2 k k k k k x y x U U U R F G t δ δ δ − − + = − − + + = Δ ) 0 k k H 。現在令 1 p p A = Γ− A 、 1 p p B = Γ− B 、 1 p p C = Γ−C 並將其分為兩部分: p p p A =A+ +A 、− B_p =B+_p +B 、−_p C_p =C+_p+C −_p 其中 1( | | 2 p p A ) A ± = A+± λ I 、 1 ( | | 2 p p B ) B ± = B+± λ I 、 1 ( | | 2 p p C C ± = C+± λ I) A λ、λ_B、λ_C分別為A 、p B 、p C 中最大的特徵值。 p

(42)

將上述帶入(3-56)得 ~ ~ ~ ~ , , 1 , 1 , 1 1 ~ ~ ~ ~ ~ ~ ~ ~ , , 1 , 1 , 1 , , 1 , 1 , 1 1 3 [ 2 ] p i p i p i p i p i p i p i p i p i p i p i p i k p A A A A M t x x B B B B A A C C U R y y z z + + − + − + − − + + − + + + − + − + − − + − − − − Γ + + + Δ Δ Δ − − − − + + + Δ = Δ Δ Δ Δ Γ _(3-75) 其中 ~ ~ ~ ~ ~ ~ ~ ~ , , 1 , 1 , ( ) p i p i p i p i x p p x p x p A A A A A A A A x x δ δ δ 1 + + − − + +₊ − ₌ − +₊ + − ₌ − ₊ − Δ Δ + − 可將(3-75)式整理成 1 (L+ +D U)ΔU_p = Γ− Rk 其中 1, , , 1, , , 1 1 1 1 ( _{p i}) _{j k} ( _{p i j}) _k ( _{p i j k}) L A B C x y z + + + − − ⎡ ⎤ = −_⎢ + + _⎥ Δ Δ Δ ⎣ ⎦ − 1 , , , , , , , , , , , , 3 1 1 1 ( ) ( ) ( ) ( ) ( ) ( ) 2 p i j k p i j k p i j k p i j k p i j k p i j k I D M A A B B C C t x y z τ − + − + − + − ⎧ _⎡ _⎤ _⎡ _⎤ _⎡ ⎫ = + Γ +_⎨ _⎣ − _⎦+ _⎣ − _⎦+ _⎣ − ⎤_⎬ Δ Δ _⎩Δ Δ Δ ⎦_⎭ 1, , , 1, , , 1 1 1 1 ( _{p i}) _{j k} ( _{p i j}) _k ( _{p i j k}) U A B C x y z − − − + + + ⎡ ⎤ =_⎢ + + _⎥ Δ Δ Δ ⎣ ⎦ 最後整理為 1 1 (L+D D) − (D U+ )ΔUk_p = Γ− Rk 上式可以用以下的步驟解出： 1.(L+D)ΔU*_p = Γ−1Rk 其中 * 1 ( ) k p p U D− D U U Δ = + Δ 2. ( ) k * p p D U+ ΔU = ΔD U * 1 k k p p U U D U U− Δ = Δ − Δ p p 3. k 1 k p p U + =U + ΔU 4.重複步驟一至三，直到 1 0 k k p p U U t + ₋ Γ = Δ _{則下一個物理時間項可求得。}

(43)

3-6、非反射性邊界(Non-reflecting Boundary)： 本研究主要探討可壓縮流在三維漸縮管道內流動與熱傳機構，過去在模擬低速可壓縮流流場時，出口條件多利用速度完全發展流條件及大氣壓力邊界，然而此上述的出口條件設定方式，由於壓力波易會在出口處反射回計算區域而影響收斂，因此不適用於低速可壓縮流；且完全發展流出口在物理模式的設定上需有一定長寬比例，將使網格需求量增加。故本文在出口條件部分採用 Fu 和 Li 等人[24]改善 Poinsot 和 Lele[25]的非反射性邊界條件(non-reflecting boundary)，可適用於極低速可壓縮流，且大幅減少網格總數及計算時間。如上述在 3-4 節中，經過 Preconditioning 法修正的 Navier-Stokes 方程式如下: 0 p U F G H t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ _(3-57) 其中Γ為 preconditioning 矩陣，U_p為[ ,P u v w T, , , ] /T J

根據 Poinsot 等人[25]所提出的 LODI(The local one-dimensional invsicid relations)將(3-57) 簡化成在出口邊界的一維 Navier-Stokes 方程式: 0 p U F t x ∂ ∂ Γ + = ∂ ∂ (3-76 ) 為了將 F x ∂ ∂ 轉換成主要變數形式(primitive form)，將(3-76)的左邊乘上 1 − Γ 1 0 p U F t x − ∂ _{+ Γ} ∂ ₌ ∂ ∂ (3-77) 其中 1 1 p 1 p p U U F F A p x U x x − ∂ − ∂ ∂ − ∂ Γ = Γ = Γ ∂ ∂ ∂ ∂ (3-78) 將(3-78)代入方程式(3-76)以主要變數形式為變數(primitive variables)的方程式如下所示： 1 0 p p U U A t x − ∂ ∂ + Γ = ∂ ∂ p (3-79) 再將 1 p A − Γ 做相似轉換以得到進出口的特徵速度：

(44)

1 p 1 A K Kλ − Γ = − (3-80) 此處K為特徵向量 (Eigen vector)，λ為 1 p A − Γ 的特徵值(Eigen values)，即為出口處的特 徵速度。根據 Dennis 等人[35]轉換低速情況下原流體速度 u 及原音速，將其調整使修 正過後的流場速度 u 與聲速 ' 冪次級數(order)接近。 c ′ c 特徵值矩陣為 1 2 3 4 5 u u u u c u c λ λ λ λ λ λ ⎞ ⎛ ⎛ ⎞ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ = = ⎟ ⎜ _{⎜ ′ ′}₊ ⎟ ⎟ ⎜ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ _⎝ _′₋ _′_⎠ ⎝ ⎠ (3-81) 其中 ( 1) 2 u u′ = β+ 與 2 2 ( 1) 4 2 u c c 2 β− + β ′ = 再令 1 Up L K x λ − ∂ = ∂ (3-82) 展開L可得 1 2 3 4 5 1 ( ) ( ) ( )[ ( ) ( )[ ( ) T P P u x x x w L _u x L v L L _u x L P u u c u c u L ] ] x x P u u c u c u x x κ ρκ ρ ρ ∂ ∂ ∂ ⎞ ⎛ ₊ ₋ ⎟ ⎜ _∂ _∂ _∂ ⎟ ⎜ ∂ ⎟ ⎜ ⎞ ⎛ ⎟ ⎜ ⎟ ⎜ _∂ ⎟ ⎜ ⎟ ⎜ _∂ ⎟ ⎜ ⎟ ⎜ = = ₋ ⎟ ⎜ ⎟ ⎜ ∂ ⎟ ⎜ ⎟ ⎜ _∂ _∂ ⎜ ⎟ ⎜ ′+ ′ − ′− −′ ⎟ ⎝ ⎠ _⎜ _∂ _{∂ ⎟} ⎟ ⎜ _∂ _∂ ′− ′ − ′+ −′ ⎟ ⎜ ∂ ∂ ⎝ ⎠ (3-83) L即為出口處，波強度隨時間變化量。為了得到L與 Up t ∂ ∂ 的關係，將(3-82)代回(3-79)並整理可得 0 p U KL t ∂ + = ∂ (3-84)

(45)

展開(3-84)可得壓力與垂直進出口的速度方程式 4 5 1 [ ( ) ( )] 0 2 P L u c u L u c u t c ∂ ₊ _′_{+ −}_′ ₋ _′_{− −}_′ ′ ∂ = (3-85) 4 5 1 ( ) 2 u L L t ρc ∂ ₊ ₋ ′ ∂ = 0 (3-86) 在出口條件部份，從出口部份往計算範圍傳遞的波強度隨時間變化量為，為了讓計算範圍不被出口部份回傳的變化量所影響，因此設定 4 L 5 0 L = 。(3-85)與(3-86)則如下式： 4 1 [ ( )] 0 2 P L u c u t c ∂ ₊ _′_{+ −}_′ ′ ∂ = (3-87) 4 1 0 2 u L t ρc ∂ ₊ ₌ ′ ∂ (3-88) 由(3-88)可得 4 2 u L c t ρ ′∂ = − ∂ (3-89) 將(3-89)帶回(3-87)可得 ( ) P u u c u t ρ t ∂ ₋ _′_{+ −}_′ ∂ ₌ ∂ ∂ 0 1 ) k + (3-90) 方程式(3-90)表示出口部份壓力與速度對時間的變化量，因此可離散(3-90)當作出口部份的壓力邊界條件。 1 ( )( k k k P + =P −ρ u′+ −c′ u u −u (3-91)

(46)

3-7、座標轉換：

本研究為了觀察此三維漸縮煙囪管道的流場及溫度場的變化及不規則邊界的影響，採用代數格點的網格建置方法(algebraic grid generation)，以座標轉換的方式，將物理平面之漸縮處拉伸展開成一規則的計算平面，如圖 3-3 所示，對於方程式座標的轉換在此節做說明。首先，原座標系統，也就是卡式座標系統，稱作 Physical domain( , , )x y z ，轉換過後的曲線座標系統則稱為 Computational domain( , , )ξ η ς 。其關係式如下： ( , , , ) ( , , , ) ( , , , ) t t x y z t x y z t x y z τ ξ ξ η η ς ς = = = = (3-92) 經由 Chain rule 對偏微分作計算，可得兩座標間的偏微分轉換公式如下： t t t x x x y y y z z z t x y z ξ η ς τ ξ η ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ∂ ₌ ∂ ₊ ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ₌ ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ ∂ ₌ ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∂ ς (3-93)

由上式可得在 physical domain 和 computational domain 之間存在一個轉換矩陣，稱為 Jacobian 轉換矩陣，三維的轉換矩陣如下： ( , , ) 1 ( , , ) ( ) ( ) ( ) J x y z x y z_ξ _{η ς} y z_{ς η} x y z_η _{ξ ς} y z_{ς ξ} x y z_ς _{ξ η} y z_{η ξ} ξ η ς ∂ = = ∂ − − − + − (3-94) 本研究主要在漸縮管道處做網格座標轉換的動作，首先為一個正 Y 方向及正 Z 方向的四分之一漸縮管道，如圖 3-4(a)所示，轉換方程式如下： 2 1 1 2 2 1 1 2 2 2 ( 2 2 2 ( ) 2 x d d d y l d d d z l ) ξ ξ η ξ ς = − = + − = + (3-95)