國 立 交 通 大 學
機械工程學系
碩士論文
三維波形管之熱傳分析
Analysis of Heat Transfer in Three Dimensional Wavy Duct
研 究 生：羅啟修
指導教授：傅武雄 博士
三維波形管之熱傳分析
Analysis of Heat Transfer in Three Dimensional Wavy Duct
研 究 生： 羅啟修 Student: Chi-Xiu Lo
指導教授： 傅武雄 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 2012
i
三維波形管之熱傳分析
研究生：羅啟修 指導教授：傅武雄 國立交通大學機械工程學系 摘要 本研究利用數值方法分析可壓縮流在三維垂直波形管道中的流動及熱傳機制。 網格採用橢圓偏微分方程法生成。流場利用有限差分法進行計算，計算方法可分為 兩部分：第一部份為非黏滯性項的尤拉方程式採用 Roe 方法計算通量，並且加入 Preconditioning 矩陣，使程式在計算低速可壓縮流可獲得良好之收斂結果，而由於使 用 Preconditioning 法時加入 Artificial time term 破壞了整個統御方程式，因此需使用 Dual time stepping 疊代使其在 Artificial domain 收斂才能進入下一個真實時階；第二 部份為黏滯性項的計算，採用二階中央插分法。在時間項方面則採用 LUSGS 隱式法， 利用 LUSGS 疊代以求出下一時階物理量。出口設非反射性邊界條件避免可壓縮流中 壓力波的干擾。由於本研究為可壓縮流之模擬，不需使用 Boussinesq assumption，適 用於溫差高於 30K 之情形，應用範圍更為廣泛。另外為提升計算速度，本文採用 OpenMP 方法進行平行化運算。 本研究以振幅波長比及瑞利數為主要參數進行模擬。對於三維波形管，流體流 經波峰時類似流經漸縮管，邊界層縮小，熱傳能力較好。而流經波谷時，流體在壁 面處近乎停滯，熱傳效率較差。在所計算之瑞利數(104 ~ 106 )和振幅波長比(0、0.1、 0.2)之間流場並無出現不穩定的情形。三維波形管在自然對流時，平均紐塞數隨振 幅波長比增加而下降，且雖然散熱面積增加，但總熱傳量反而下降。ii
Analysis of Heat Transfer in Three Dimensional Wavy Duct
Student: Chi-Xiu Lo Advisor: Wu-Shung Fu Department of Mechanical Engineering
National Chiao Tung University
Abstract
An investigation of heat transfer in a three-dimensional wavy duct 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. Due to the consideration of flow compressibility , Boussinesq assumption is not used. Therefore, applicable to cases with temperature difference higher than 30K. Besides, to enhance the computing efficiency, the OpenMP method is also used.
In this study, amplitude-wavelength ratio and Rayleigh number is set to be primary variables. The results show that when flow passes through the crests, the boundary layer gets smaller, therefore a better heat transfer rate. And when flow passes through the troughs, fluid is almost stagnant near the walls , which leads to a poor heat transfer rate. Flow field shows no instability for Rayleigh number ranging from 104~106 and amplitude-wavelength ratio 0~0.2. In all cases studied, the averaged Nusset number is decreased as amplitude-wavelength ratio increases, and the total heat flux decreases too even though the cooling area is increased .
iii 誌謝 首先感謝指導教授傅武雄老師，在大學時豐富的課程讓我對熱流領域感興趣， 以及研究所期間對於課業論文及生活上耐心的教導與提點，指引我做研究時正確的 方法及方向。同時也感謝機械系師長在這六年來課業上的指導及良好的學習環境。 另外要特別感謝實驗室的大家，感謝李崇綱、黃玠超、王威翔、黃耘、黃上豪等學 長姊的教學，研究上的問題總是能夠迎刃而解，還有同儕的鄭景木、李世豪和黃崑 榕一同為了畢業而互相打氣努力，總是在精神上支持協助的學弟們蔡承志、范家魁、 吳佩蓉，為我抒解了許多的壓力。另外也要感謝常到實驗室陪伴我們的同學劉岳儒、 陳彥志、顏士傑。最後更要感謝我的父母，給予我衣食無缺的環境，感謝您們對我 的照顧與包容，今日方能順利完成學業。最後僅將此喜悅和所有關心我的人分享。
iv 目錄 中文摘要 ... i 英文摘要 ... ii 目錄 ... iii 表目錄 ... iv 圖目錄 ... v 符號表 ... vii 一、緒論 ... 1 二、物理模式 ... 6 2-1 物理尺寸與分析模式 ... 6 2-2 分析假設及統御方程式 ... 7 2-3 邊界條件 ... 9 三、數值計算模式 ... 13 3-1 統御方程式 ... 13 3-2 橢圓 PDE 網格生成法 ... 17 3-3 Roe scheme ... 26 3-4 MUSCL 法 ... 34 3-5 Preconditioning ... 36 3-6 黏滯性項之差分 ... 40
3-7 Dual time stepping ... 43
3-8 LUSGS 法 ... 45
3-9 反射性邊界 ... 47
四、結果與討論 ... 51
五、結論 ... 90
v
表目錄
表 3-1 精度係數值 ... 35 表 4-1 二維波形板改變振幅波長比之結果 ... 63 表 4-2 三維波形管散熱面積之比較 ... 84
vi 圖目錄 圖 2-1(a.b) 波形管道之物理模式 ... 11 圖 3-1 對邊界加密之網格 ... 20 圖 3-2 對邊界正交化之網格 ... 21 圖 3-3 沿 X 方向波形壁面之管道 ... 22 圖 3-4 管道在 Z=0.5 的 XY 平面網格分布 ... 23 圖 3-5 管道在波谷處之 YZ 平面網格分布 ... 24 圖 3-6 管道在波峰處之 YZ 平面網格分布 ... 25 圖 3-7 黎曼問題特徵值結構圖 ... 33 圖 3-8 差分示意圖 ... 42 圖 3-9 L1、L2、L3、L4與L5於管道兩端的方向性示意圖 ... 50 圖 4-1(a.b) 和 Guzman[26]結果之比對圖 ... 58 圖 4-2 和 Bhavnani[2]實驗對照所設定之物理模式 ... 60 圖 4-3 和 Bhavnani[2]振幅波長比 0.1 之局部熱傳係數比對圖 ... 61 圖 4-4 和 Bhavnani[2]及 Ostrach[28]振幅波長比 0 之局部熱傳係數比對圖 ... 62 圖 4-5 三維波形管網格測試圖 ... 64 圖 4-6(a~c) Ra=105 ,α=0 管道中央XY截面之結果 ... 65 圖 4-7(a.b) Ra=105 ,α=0 管道YZ截面之速度場 ... 66 圖 4-8(a~c) Ra=105 ,α=0 管道中央XY截面之結果 ... 67 圖 4-9(a.b) Ra=105 ,α=0.1 管道YZ截面之速度場 ... 68 圖 4-10(a~c) Ra=105 ,α=0.2 管道中央XY截面之結果 ... 69 圖 4-11(a.b) Ra=105 ,α=0.2 管道YZ截面之速度場 ... 70 圖 4-12(a.b) a=105沿X方向紐塞數比較圖 ... 71
vii 圖 4-14(a~c) Ra=5.5*104流線圖 ... 73 圖 4-15(a~c) Ra=5.5*104速度場圖 ... 74 圖 4-16(a~c) Ra=5.5*104等溫線圖 ... 75 圖 4-17(a~c) Ra=2.5*104流線圖 ... 76 圖 4-18(a~c) Ra=2.5*104速度場圖 ... 77 圖 4-19(a~c) Ra=2.5*104等溫線圖 ... 78 圖 4-20 振幅波長比 0 時沿 x 方向管道中央壁面局部紐塞數圖 ... 79 圖 4-21 振幅波長比 0.1 時沿 x 方向管道中央壁面局部紐塞數圖 ... 80 圖 4-22 振幅波長比 0.2 時沿 x 方向管道中央壁面局部紐塞數圖 ... 81 圖 4-23 平均紐塞數對瑞利數之趨勢圖 ... 82 圖 4-24 總熱傳量對瑞利數之趨勢圖 ... 83 圖 4-25(a~c) 振幅波長比 0.1 之流線圖 ... 85 圖 4-26(a~c) 振幅波長比 0.1 之等溫線圖 ... 86 圖 4-27(a~c) 振幅波長比 0.2 之流線圖 ... 87 圖 4-28(a~c) 振幅波長比 0.2 之等溫線圖 ... 88 圖 4-29 平均紐塞數對瑞利數之趨勢圖 ... 89
viii
Nomenclature
a wavy amplitude [m]
p
C constant-pressure specific heat [J kg⋅ −1⋅k−1]
v
C constant-volume specific heat [J kg⋅ −1⋅k−1]
d width of the vertical tube [ m ]
e internal energy [J kg⋅ −1] g acceleration of gravity [m s⋅ −2]
h Enthalpy [J ]
x
h local heat transfer coefficient [Wm-2K-1]
h average heat transfer coefficient [Wm-2K-1] k thermal conductivity [Wm-1k-1]
l length of the vertical tube [ m ]
Nu
Nusselt number defined in Eq.(4-6)x x h d Nu k =
Nu
Paverage Nusselt number 1 A Nu Nudxdy A =
∫∫
Pressure [Pa]q total heat flux (watt)
R gas constant [J kg⋅ −1⋅k−1]
ix 2 3 0 2 ( ) Pr Pr ( ) w g T T d Ra Gr T T r µ − = ⋅ = T Kelvin temperature[K] 0 T Surrounding temperature[K] w
T temperature of heat surface[K]
u velocity component in x -direction[m s⋅ −1]
v velocity component in y-direction[m s⋅ −1]
w velocity component in z -direction[ 1
m s⋅ − ] x , y , z Cartesian coordinate (m)
x Greek symbols α amplitude-wavelength ratio a a λ = r density[kg m⋅ −3] 0 r surrounding density[kg m⋅ −3] υ kinematics viscosity[ 2 1 m ⋅s− ] µ absolute viscosity[kg.m−1.s−1] γ λ
specific heat ratio wavelength (m) w
τ friction force acting per unit area on the surface
, ,
1
第一章、緒論
近年來，科技產業蓬勃發展，電子設備及產品為提高競爭力，不斷的在效能 上衝刺，加上半導體製程技術之進步，促使消費性電子產品的使用與生活型態做 結合，因而許多電子產品的設計走向高效能、微小化的趨勢，致使電子元件舉凡 電腦 CPU、顯示卡、筆記型手機等其封裝元件之發熱密度越來越高，其單位熱通 量亦相對地不斷增加，對性能及可靠度等方面均造成不容忽視之影響。因此在未 來產品的發展趨勢走向更輕薄短小之際，散熱技術對於其性能及可靠度的提升將 扮演重要的角色。一般正常電腦中央處理器在執行程式下溫度大約為攝氏 50~60 度之間，超過攝氏八十度則電腦基於保護狀態，自動關閉系統，長期下來內部零 組件壽命將會減短。因此為了維持元件於額定溫度下運作，必須將此密集的熱量 能有效散逸於系統外之環境。傳統的散熱片多使用鋁合金製造，其熱傳導性只屬 於中等程度，對於今日元件發熱功率越來越高的情況，逐漸有捉襟見肘的情況發 生，目前電子元件的發熱量達到每平方公分數十瓦的等級，且接點可承受之溫度 約在攝氏 150 度以下。因此如何在有限之條件下提升熱傳效率為當前重要之研究 課題之一。 一般提高熱傳效率，增加主動式散熱元件功率為最直接的選擇，但如此往往 伴隨高噪音高耗能或體積增加等之缺點，因此在工程上應用時，通常會藉由改變 散熱區域之外型以達到同功率輸入下之最佳熱傳，其重要之設計原理除了增加散 熱面積、降低熱阻抗、減少流動所需之壓差外，設計之外形能否有效干擾熱交換 面邊界層形成使壁面高溫和外界低溫流體混合也是重要參考之因素之一。而其中 ㄧ常見之外型是波形壁面，波形壁面在工業上應用廣泛，在許多產品中都可見到 其身影，如板形熱交換器(plate heat exchanger)、半透膜加氧器(membrane oxygenators)、電透析器(electrodialyzer)等，由於波形壁面在波谷處容易引 起渦流，使得流體混合效率提升並破壞熱邊界層，因此可有效增加熱傳效率。 探討波形壁面之文獻不少，其中自然對流的部分，Yao[1]利用座標轉換的方式模2 擬二維垂直波形壁面於開放空間之自然對流現象，其採用 Boussinesq 假設，壁 面振幅對波長比 0~0.3，結果顯示局部紐塞數對流動方向呈正弦函數震盪，波長 為壁面波長之一半，振幅則隨流動方向遞減。而平均紐塞數部分，所測試之案例 皆略低於平板。Bhavnani 和 Bergles[2]以實驗方法觀測等溫波形板之熱傳，其 研究結果和 Yao[1]相異，局部紐塞數震盪之週期和壁面週期相同，他們推測此 差異可能是由於 Yao 在統御方程式做了過多之近似。隨後 Moulic 及 Yao[3]針對 波形壁面之混合對流進行模擬，試圖解釋和實驗結果差異的原因，他們指出自然 對流的浮力效應會造成類似混合對流的現象，而模擬結果顯示混合對流在平均紐 塞數圖出現兩個諧波(two harmonics)震盪的現象，分別和自然對流及強制對流 相關。Ashjaee 等人[4]使用干涉計測量雷射束通過流體時產生之相位差(phase shift)得到波形壁面自然對流之溫度場，結果在壁面振幅波長比為 0.2、 Ra=5.59*10^5 的案例中，局部熱傳系數具有兩個諧波震盪的現象，但在壁面振 幅波長比低於 0.2 的案例中第二諧波並不明顯，在這篇研究中也利用 SIMPLE 法 模擬，其結果和實驗結果相當接近，並且將結果之局部紐塞數對流動方向座標、 瑞利數和壁面振幅波長比製作成經驗方程式。 而在波形壁面強制對流的部分， Wang和Vanka[5]利用橢圓PDE生成網格，將 二維管道映射至正交網格計算，流動方向採週期性邊界，探討其暫態及穩態熱傳， 結果顯示當雷諾數低於 180 時，其流場為穩態層流，高於此雷諾數，則流場會出 現不穩定的情形。而隨著雷諾數上升，非線性之對流項影響逐漸提高，因此流場 由上下對稱趨於不對稱，另外在上壁面波鋒處（下壁面波谷處）可觀察到渦流之 現象，當雷諾數高於 180 時，流場不穩定的現象使壁面渦流週期性生成消散，破 壞熱邊界層，因此熱傳效率較平板高。Rush等人[6]利用可視法觀察波形管道內 過渡流場發生位置與雷諾數、振幅、波長、上下壁面相位差之關係，實驗結果發 現即使入口速度為穩定且均勻，當雷諾數高於一過渡雷諾數時，流場即會產生不 穩定的現象，此現象發生之位置隨雷諾數增加往入口方向移動，且和壁面相位差、
3 相對振幅有關，而在測試案例中相對波長對此位置影響不大。Pham 等人[7]利用 大渦流模擬(LES)配合沉浸式邊界計算波形壁面之流場及熱傳現象，其物理模式 在流動方向設定為八個波長，藉以觀測流場由層流轉變為紊流之過程。測試之雷 諾數範圍為 750 到 4500，當雷諾數在 1500 以下，流場雖可觀測不穩定的現象但 仍維持層流，雷諾數 3000 以上則出現明顯紊流的特徵。而壓降及熱傳分析顯示， 間隙比(channel spacing ratio 振幅除以上下壁面間距)越低則速度場越均勻且 紊流強度也越小，Fanning factor隨間隙比減少而下降，Colburn factor和間隙 比較無關，最佳j/f比發生於流場產生不穩定現象但未進入紊流時。由結果可得 縮小間隙比對流場分離及渦流產生有抑制之效果，能在相同之熱傳量下以較低之 壓降進行。Sun等人[8]利用 直接數值模擬法(DNS)配合座標轉換模擬波形壁面之 紊流現象，繪出鄰近壁面之平均速度及紊流強度。Alawadhi和Bourisli[9]使用 套裝軟體FLUENT模擬波形壁面管道之流場，其數值方法採用SIMPLE法，入口之邊 界條件為等速流外加一正弦振盪之擾動，出口為定壓，結果顯示當速度擾動和渦 流生成之頻率一致時，熱傳效率最好，在特定條件時熱傳效率比沒有擾動提高 60%。SAMAD等人[10]使用套裝軟體ANSYS-CFX 11.0 計算交錯排列凹槽壁面管道之流場， 數值方法採用RANS，以無因次化參數凹槽深除以凹槽半徑、管寬除以凹槽半徑、凹槽
半徑除以凹槽間距為變數，利用polynomial RSA surrogate model進行最佳化，結果得
到凹槽深除以凹槽半徑越大時，其所定義的熱傳表現越好。Haitham[11]探討穩態完 全發展流流經二維波形及弧形管道之流場，在特定條件下，熱傳能提升 80%，但 大部分時熱傳表現都低於直管。Mohammed[12]等人 模 擬三維波形微渠道散熱器， 結果得到振幅波長比低於 0.25 的例子時，振幅越大平均熱傳係數也越好，但振 幅波長比 0.25 則反而下降且低於直管微渠道。綜觀前述之研究，波形壁面之熱 傳效率和相關參數影響甚大，某些參數設定下，波形壁面熱傳和平板熱傳相比並 無顯著改善，但如果當參數調整合宜使流場出現不穩定的現象時，熱傳效率則能 有效提升，此外臨界雷諾數也會因相關參數影響，從穩定層流至不穩定層流的臨
4 界雷諾數約在 200~800 之間，而進入紊流的臨界雷諾數應大於 1500。 目前波形壁面之模擬文獻，多需借助許多假設，如流體為不可壓縮，使用 Boussinesq assumption 以模擬密度差所造成之浮力等，這些假設限制了程式使 用之範圍，在工程上應用時容易造成不可忽視之誤差。此外為節省網格使用量， 早期之研究多假設為二維流場，因此實際三維流場之流動現象無法仔細探討，而 在邊界條件上，也常見在不適合處使用完全發展流邊界條件。 本研究針對三維之流場作模擬，工作流體為可壓縮流，密度受溫度及壓力而 改變，不需使用 Boussinesq assumption，因此適用於高溫差及高速之案例，對 於工業之應用較為廣泛及實際。在波形壁面變形部分採用 Thompson J.F.等[13] 及 Hoffmann[14]之座標轉換式，配合橢圓 PDE 網格生成法生成，由於本研究為 三維管道，在流動方向網格一般採等間距分布即可，只需對垂直流向之平面調整， 因此使用二維橢圓 PDE 網格法生成並針對邊界處加密。程式主體參考李[15]發展 出之黏性流場之全域速度場數值解法，此種方法最大的困難處在於計算低速流場 時，由於可壓縮流必須遵守 Courant-Friedrichs-Lewy (CFL)條件，因此在低速 流體時，受限於流體變化傳遞速度(約等於聲速)，時階將會極小。在此種情況造 成計算過程將耗費極大的時間與整體過低的計算效率。為了改善此缺點，本研究 在計算可壓縮流時加入 Preconditioning 法，藉此讓流體即使在低速時，也可 有較高的效率與良好的收斂性。在計算此種低馬赫數流體的方面，目前有密度基 底法(Density-based method)與壓力基底法(Pressure-based method)。本研究 以密度基底法為主。而在密度基底法中又以 Turkel[16]提出 Preconditioning 法最為廣泛應用，不僅可同時應用在可壓縮流與不可壓縮流中，更可以讓程式的 收斂性增加。在邊界設定最棘手的問題是流體的流動速度和壓力波的速度相差過 大，導致在出口反射回彈的壓力波會干擾流體流動，故在入口及出口處設置非反 射性(non-reflecting)邊界條件[17]，避免干擾發生，若設定為完全展開流，則 會與實際流場有所差異。
5
故本文需要求解完整的 Navier-Stokes 方程式與理想氣體方程式(ideal gas equation)以得此可壓縮流中密度的變化，以期能同時考慮密度與壓力變化之效 應符合實際物理現象。採用中央插分法處理黏性項部分；而非黏性項則必須採用 Roe 法[18]，以解決偏微分方程式的不連續問題；本研究在作數值計算時，主要 是 利 用 網 格 之 間 的 物 理 量 ， 因 此 使 用 MUSCL(Monotone Upstream-centered Schemes for Conservation Laws)法[19]來計算 ROE 法中網格與網格間的物理量； 為了處理在低速可壓縮流時不易收斂的困擾，計算時加入了 Preconditioning 法[20]，可有效提高收斂效率，使其在高速與低速可壓縮流場均可適用，但卻破 壞原統御方程式，為了彌補此一缺點，故需要加入 Artificial time term 來修 正 ； 並 使 用 Dual time stepping[21] 來 計 算 暫 態 的 物 理 量 ； 最 後 使 用 LUSGS(implicit lower-upper symmetric Gauss-Seidel algorithm)法[22]做迭 代運算。 此外，在處理複雜的流體力學問題時需要大量的計算過程，利用多核心處理器來 提升運算速度已為目前的發展主流。多核心處理器對於單一執行緒在平行運算方 面，並無法提升計算速度；若利用多執行緒的程式架構，可透過不同核心來同時 計算，達到提升計算效率、節省計算時間之目的；但是多執行緒的程式在撰寫、 編輯上，也都比單一執行緒的程式架構複雜。常見的平行運算方法有 MPI 和 OpenMP 兩種方法；本文採用 OpenMP 方法提升運算速度。MPI (Message Passing Interface)是一種分散式記憶體(distributing memory)的觀念，而 OpenMP (Open Multi-Processing)為共享式記憶體(sharing memory)的觀念，而兩種方 法各有優缺點。MPI 在撰寫程式上面相較於 OpenMP 較困難，且計算速度較慢且 會受限於網路效率，在設備擴充方面較昂貴；反之 OpenMP 在程式撰寫上較簡單， 計算速度快且不會受限於網路效率，在設備擴充方面需較少的經費。因此本程式 利用 OpenMP 來進行平行化運算，效率為原程式的四倍，減少計算時間及成本。
6
第二章 物理模式
2-1、物理尺寸與分析模式： 本文之物理模式如圖 2-1 所示，為一正方形截面之三維管道，由上而下可分 為三個段落：上直管、波形管、下直管。上下直管截面之邊長為 d，長度上直管 為λ、下直管 λ，邊界設定為絕熱壁面、不可滑移條件。波形管道段落長度為3λ， 由振幅 a、波長λ
之正弦函數三個週期組成，邊界設定為等溫壁面、不可滑移條 件。外界環境及流體初始條件設定為壓力 0P
=101300Pa、溫度T
0=298.0592K，受 到波形壁面高溫 wT
影響，密度改變造成浮力效應驅動流場，使得流體進出上下 管道邊界．以往大多將此截面設為速度完全發展流條件，但此方式在可壓縮流中 壓力波於出口處易反射回計算區域而影響收斂性，較不適用於低速可壓縮流流 場；且為使熱管道內的熱流場更容易成為完全發展流，管道的長寬將有所限制， 因此常使得網格需求增加，計算時間較長，因此本研究在出口條件部分設為非反 射性(non-reflecting)邊界條件[14]，此時可視其出口為極遠處的邊界條件，大 幅減少使用網格數量。 此外定義無因次參數振幅波長比a
a
λ
=
。另外如圖 2-1a，以管道中央為基 準，定義壁面之波峰波谷。7 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 22
xx xy xz xx xy xzu
u
P
uv
F
uw
V
T
e
u
Pu
k
u
v
w
x
ρ
ρ
τ
ρ
τ
ρ
τ
ρ
τ
τ
τ
+ −
−
=
−
∂
+
+
−
−
−
−
∂
8 − − − ∂ ∂ − + + − − + − = 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。 1 2 2 2 ( ) 2 v e=C T+ u +v +w ，Cv為等容比熱。 理想氣體狀態方程式： 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 ρ ρ τ ρ τ ρ τ ρ τ τ τ − − = + − ∂ + + − − − − ∂ 9 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：非反射性邊界 入口速度 v：非反射性邊界 入口速度 w：非反射性邊界 入口壓力 p：非反射性邊界 入口溫度 T：非反射性邊界 2-3.3 上出口邊界條件： 出口速度 u：非反射性邊界 出口速度 v：非反射性邊界 出口速度 w：非反射性邊界 出口壓力 p：非反射性邊界 入口溫度 T：非反射性邊界 2-3.4 波形壁面邊界： 邊界速度： 不可滑移條件，u =v=w=0m/s 邊界溫度： 加熱壁面溫度 w
T
邊界壓力：在垂直壁面方向，梯度為零， p 0 n ∂ = ∂10 2-3.5 上直管壁面邊界： 邊界速度： 不可滑移條件，u =v=w=0m/s 邊界溫度： 絕熱壁面，垂直壁面方向溫度梯度為零， T 0 n ∂ = ∂ 邊界壓力：垂直壁面方向，梯度為零， p 0 n ∂ = ∂ 2-3.6 下直管壁面邊界： 邊界速度： 不可滑移條件，u =v=w=0m/s 邊界溫度： 絕熱壁面，垂直壁面方向溫度梯度為零， T 0 n ∂ = ∂ 邊界壓力：垂直壁面方向，梯度為零， p 0 n ∂ = ∂
11
12
13
第三章 數值計算模式
本章主旨在說明本論文的數值計算所使用的所有模式。第一節介紹統御方程式 Navier-Stokes 方程式。參考 Fu 和 Li 等人[23]之數值方法，將方程式拆解為非黏滯項與 黏滯項。而由於研究主題為曲線管道，為簡化計算過程將物理空間映射至計算空間 求解，因而統御方程式也須經過轉換。第二節介紹橢圓 PDE 網格生成法[14]及利用差 分法求得轉換導數。第三節介紹的為黎曼解中的 ROE 法(Roe Scheme)[18]，利用 ROE 法來求出非黏滯項的通量。接著第四節介紹 MUSCL 法(Monotone Upstream-centered Schemes for Conservation Laws)[19]，此法是為了要解出 ROE 法中使用的網格之間的物 理量，然後為了防止在高階插分時產生震盪現象，在 MUSCL 法插分的結果方程式中 加入 Minmod limiter 以確保程式不會發散。第五節為介紹 Preconditioning 法[20]，因為 當計算低速可壓縮流時，因速度和音速的數量級(order)上差距過大，在數值分析時不好 計算，所以為彌補此一缺點須使用 Preconditioning 法。第六節介紹計算黏滯性項則採 用二階精度的中央差分法。第七節為 Dual time stepping，因使用 Preconditioning 法時， 加入 Artificial time term 而破壞了完整的統御方程式，為了使計算暫態結果較準確，因 此需使用 Dual time stepping 疊代，使其在 Artificial domain 收斂時才能進入下一個真實 時 階 ， 更 提 高 程 式 的 效 率 。 第 八 節 為 LUSGS 法 (implicit lower-upper symmetric Gauss-Seidel algorithm)[22]，程式因為在使用 Preconditioning 法，故而在統御方程式中 修改了計算時階，故需要加入 Artificial time term，待時階項收斂時才能達到真實穩態。 第九節為完全非反射邊界之介紹，為了使自然對流之邊界可以更接近真實狀態，在程 式的出入口邊界使用完全非反射邊界。綜合上述，本論文在數值上的計算過程為，利 用 MUSCL 法算出 ROE 法所需要的網格間物理量，求出非黏滯的通量，並且在計算通 量時加入 Preconditioning 法，以拉近與音速的數量級(order)。接下來使用二階中央插分 法對黏滯項做插分進而求出黏滯項；然後再與 ROE 法求出的通量結合得到真正的物理 通量。最後使用 Dual time stepping 及 LUSGS 法疊代以求出正確時階的物理量。14 3-1、統御方程式： 本研究在計算流場的方面其統御方程式分可為兩大部分，第一部份為非黏滯性項 的尤拉方程式，第二部份為黏滯性項。 S z H y G x F t U = ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ (3-1) 其中
(
)
T U = ρ ρu ρv ρw ρe (3-2) 2 22
xx xy xz xx xy xzu
u
P
uv
F
uw
V
T
e
u
Pu
k
u
v
w
x
ρ
ρ
τ
ρ
τ
ρ
τ
ρ
τ
τ
τ
+ −
−
=
−
∂
+
+
−
−
−
−
∂
−
−
−
∂
∂
−
+
+
−
−
+
−
=
yz yy yx yz yy yxw
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 22
zx zy zz zx zy zzw
wu
wv
H
w
P
V
T
e
w
Pw k
u
v
w
z
ρ
ρ
τ
ρ
τ
ρ
τ
ρ
τ
τ
τ
−
−
=
+ −
∂
+
+
−
−
−
−
∂
(3-3) (3-4) (3-5)15
ρ
為密度，p為壓力。u、v、w分別為x、y、z方向的速度。k為 thermal diffusivity。 2 2 2 1 ( ) 2 v e=C T+ u +v +w ，C 為等容比熱。 v 上式可拆解為黏滯性項與非黏滯性項： + + − ∂ ∂ − + + + = + = 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 因此統御方程式可改寫如下inv inv inv v v v
E
F
G
E
F
G
Q
S
t
x
y
z
x
y
z
∂
∂
∂
∂
∂
∂
∂
+
+
+
=
+
+
+
∂
∂
∂
∂
∂
∂
∂
左式由非黏滯項組成的部分稱為尤拉方程式。 為了簡化曲線模型易造成的複雜計算過程，將物理空間( , , )x y z 映射至計算空間 ( , , )ξ η ζ 求解，利用下列關係式 (3-6) 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) (3-10)16 t t t x x x y y y z z z
t
x
y
z
ξ
η
ζ
τ
ξ
η
ζ
ξ
η
ζ
ξ
η
ζ
ξ
η
ζ
ξ
η
ζ
ξ
η
ζ
ξ
η
ζ
∂
∂
∂
∂
∂
=
+
+
+
∂
∂
∂
∂
∂
∂
=
∂
+
∂
+
∂
∂
∂
∂
∂
∂
=
∂
+
∂
+
∂
∂
∂
∂
∂
∂
=
∂
+
∂
+
∂
∂
∂
∂
∂
將(3-11)代入(3-10)改寫，得 v v v E F G Q E F G Sτ
ξ
η
ζ
ξ
η
ζ
∂ ∂ ∂ ∂ +∂ +∂ +∂ = + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ 其中 1 ( ) 1 ( ) 1 ( ) 1 ( ) ( ) ( ) t x y z t x y z t x y z v x v y v z v v x v y v z v v x v y v z v Q Q J E Q E F G J F Q E F G J G Q E F G J E E F G J Q F E F G J Q G E F G J ξ ξ ξ ξ η η η η ζ ζ ζ ζ ξ ξ ξ η η η ζ ζ ζ = = + + + = + + + = + + + = + + = + + = + + 式子(3-12)即為統御方程式在計算空間之形式。 (3-11) (3-12)17 3-2、橢圓 PDE 網格生成法： 3-1 小節介紹了流場的統御方程式，但由於其解析解並不容易求得，因此對於較複 雜的問題，常藉由數值方法估算，將偏微分導數改寫為有限差分的形式進行求解。而 在求有限差分時，需要將計算區域分割為許多離散的點，藉由這些離散的點之值估算 偏微分導數，因此標示這些點的位置是在求解前必要的準備工作。 通常為了方便計算，希望這些格點是在一方型區域內沿著正交的網格線等間距分 布，但是對於曲線系統，這樣的網格邊界處理不易且容易造成誤差。為了克服這樣的 問題，一般會將其分為物理領域和計算領域，藉由轉換的方法將物理領域非正交的網 格轉換到一正交且等間距的計算領域來求解，只要設定好計算領域和物理領域網格一 對一對應關係，計算領域網格所求得的數值即可視作為物理領域相對應位置的數值。 但這樣的做法會產生一重要的議題，如何求得物理領域的網格點？也就是如何求得每 個計算領域的網格其在物理領域所對應的座標，一個好的對應關係必須具備幾個特 徵，一對應關係必須是一對一，也就是同族的網格線不能互相交錯，二網格線最好是 平滑且偏斜越少越好，另外網格最好易於調整，像在數值高變化處能增加網格的密度 等等。本研究中所採用的生成法為橢圓 PDE 生成法[14]，和代數法相比，橢圓 PDE 網 格生成法具有許多優點，像是網格較為平滑、可由邊界網格生成內部網格、易於調整 網格之正交性等等，因此橢圓 PDE 網格生成法應用之範圍廣泛。其原理為解 Laplace 方程式或 Poisson 方程式時，從等高圖可描繪出等高線，如果利用兩組不同軸向之等高 線則可交織出一完整之 2D 網格。其公式如下：
0
0
yy zz yy zzη
η
ζ
ζ
+
=
+
=
將獨立變數和因變數互換，則可得2
0
2
0
ay
by
cy
az
bz
cz
ηη ηζ ζζ ηη ηζ ζζ−
+
=
−
+
=
其中 (3-13) (3-14)18 2 2 2 2 a y z b y y z z c y z η η ξ η ξ η ξ ξ = + = + = + 在解上列方程組時，必須要輸入一初始條件，利用下列方程式生成 3 3 3 3 2 3 3 3 3 2 * 2 ( ) 2 ( )
* *sin[ ] 2* * *sin[ ] for l l l * for otherwise
2 ( ) 2 ( )
* *sin[ ] 2* * *sin[ ] for l l l * for otherwise x x l x l y d a a x y d x l x l z d a a x z d l ξ π π η η l l η π π ζ ζ l l ζ = − − = + − < < + = − − = + − < < + = 而轉換導數(transform derivatives)則利用有限差分法得到，內部點使用中央差分，在邊 界處之點則採用前進／後退差分求得。進行疊代至收斂即可求得網格。而對於需要邊 界加密或需要網格正交的問題則可將生成網格之統御方程式改寫為下式
( , )
( , )
yy zz yy zzP
Q
η
η
η ζ
ζ
ζ
η ζ
+
=
+
=
其中邊界及局部加密之Ｐ和Ｑ函數如下 1 2 2 2 1 1 1 2 2 2 1 1 exp( | |) exp{ [(| |) (| |) ] } | | | | exp( | |) exp{ [(| |) (| |) ] } | | | | M N m n m m m n n n n m m n n M N m n m m m n n n n m m n n P a c b d Q a c b d η η η η η η η η ζ ζ η η η η ζ ζ ζ ζ ζ ζ η η ζ ζ ζ ζ ζ ζ = = = = − − = − − − − − − + − − − − − = − − − − − − + − − −∑
∑
∑
∑
其中 a、b、c、d 為用以調整之常數。 經過加密後的網格範例如圖 3-1，在接近壁面處由於邊界層的影響，在數值上會有較大 的變化，因此對其加密。而網格正交化可分為 Neumann orthogonality 及 Dirichlet orthogonality，使用 Neumann orthogonality 調整的網格範例如圖 3-2，根據 Thompson 等人的研究，網格偏斜容易造 成截斷誤差(truncation error)，尤其是在邊界處影響特別大，因此針對邊界處調整，使
19 網格線垂直。 一般而言，對於管道流，流動方向之網格多以平均分布即可，因此適合使用 2D 之 橢圓 PDE 網格對於垂直流向的每個平面調整。以 (圖 3-3)之管道為例，將ｘ方向切為 150 個 YZ 截面，輸入 YZ 截面邊界之座標進行疊代，則可得到如(圖 3-4~3-6)之網格分 布。而對於需要使用到 3D 橢圓 PDE 網格生成法之問題則可參考 Upender[24]。
20
21
22
23
24
25
26 3-3、Roe scheme： 在雙曲線的守恆形式方程式中，若其初始條件包含有不連續的片段連續(piecewise) 常數，此類型的問題通稱為黎曼(Riemann)問題。因為其包含有不連續解，因此在流體 計算上有著相當廣泛的應用。一維線性黎曼方程式如下： 0 u u a t x ∂ + ∂ = ∂ ∂ (3-18) 其中a為一常數 Jacobian 矩陣。 初始條件為
( )
, 0 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− ，其中Λ為特徵值矩陣： 10
0
m l l Λ = 。 (1) ( ) , , m T K = K K 為特徵向量，故 ( )i ( )i i AK =lK 。 接著定義特徵變數W〈characteristic variables〉，其定義如下： ( , ) W =W t x ， ( 1) W =K − 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-19)方程式(3-19)稱為 canonical form 或 characteristic form。 將以上的結果簡單整理如下：
27 0 i i i W W t l x ∂ + ∂ = ∂ ∂ ，或 1 1 1 2 2 3 3 ... 0 0 ... 0 0 : : : : : 0 ... m t x w w w w w w l l + = (3-20) 上式可由特徵曲線法求得其解為： (0) 0 ( , ) ( ) 0 i i i i i i i x t w x t w x t x t α l l β l − < = − = − > (3-21) 其中，αi與βi為初始值的特徵變數。由於U=KW，可以得到 (0) ( ) ( , ) ( ) m i i i i u x t =
∑
w x−lt K 參照圖(3-7)，可以進一步推導出 ( ) ( ) 1 1 ( , ) p m i i i i i i p U x t α K βK = = + =∑
+∑
(3-22) 除此之外，還可決定出U x t( , )中的 jump ∆U： ( ) 1 m i R L i i U U U α K = ∆ = − =∑
(3-23) 其中αi =β αi− i。 在一維線性黎曼問題中，雖然有 exact solution，但在非線性問題裡需利用疊代等方 法，這些動作將耗費大量的時間，因此在實際應用上並不廣泛。為了解決此問題，一 般皆求解近似黎曼問題〈approximation Riemann problem〉解而不直接求其 exact solution。在求解近似黎曼問題中最被廣泛應用的方法為 Roe 所提出，亦即為 Roe scheme，其內容如下： 假設一維尤拉方程式： 0 U F t x ∂ +∂ = ∂ ∂ (3-24) 根據 chain rule，可將方程式(3-24)改寫如下：28 0 U F U t U x ∂ + ∂ ∂ = ∂ ∂ ∂ 再令A U( ) F U ∂ = ∂ ，於是方程式(3-24)可以表示成： ( ) 0 U U A U t x ∂ + ∂ = ∂ ∂ (3-25) 其中，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-26) 於是前述方法可以得到(3-26)的近似解。由以上的原理可得知，在近似黎曼問題 上，Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性轉變成線 性，但是初始條件並沒有改變，因此可以得到方程式(3-26)的近似解。為了要求得合理 的常數 Jacobian 矩陣，須合乎 Roe 所提出的四項條件： 1. U 與F 之間，存在著線性轉換的關係。 2. 當UR −UL → ，則U A U U( L, R)→A U( )，此處A F U ∂ = ∂ 。 3. A U( L−UR)=FL−FR。 4. 矩陣A的特徵向量必須線性獨立。 這四項條件都是雙曲線方程式所需具備的，這同時也說明了 Roe 所提出的常數 Jacobian 矩陣必須有實數特徵值，其所對應的特徵向量必須線性獨立。除此之外，條件 3.則是為了符合守恆定律（conservation law）與 Rankine-Hugoniot 條件。
線性黎曼問題的解析解，可以直接從(3-21)與(3-23)式得到， 1 2 ( / ) i U x t + 的解可以利 用下面的方程式計算：
29 ( ) 1 0 2 ( / ) i i L i i U x t U K l α + = +
∑
< (3-27) 或 ( ) 1 0 2 ( / ) i i R i i U x t U K l α + = −∑
> (3-28) 其中 1 2 i+ 表示網格與網格之間的交界面(face)。 而黎曼問題的近似解，則須從解近似黎曼問題著手： ( ) 0 U F U t x ∂ +∂ = ∂ ∂ ，根據(3-26)式可得知F = AU 為了符合守恆的條件，因此下式必須成立： ( R) ( L) ( R) ( L) F U −F U =F U −F U (3-29) 接著在固定體積的條件下，積分近似解 1 2 (0) i U + ，可得到通量(flux)的數值公式： 1 1 2 2 ( (0))+ ( R) ( R) i i F F U F U F U + = + − (3-30) 再從F = AU 的關係中可進一步求得： 1 1 2 2 (0)+ ( R) R i i F AU F U AU + = + − (3-31) 再根據(3-27)式與(3-28)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i R i i i i F F U A K F U K l α l α+ + = − ∑
> = −∑
= (3-32) 或 ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i L i i i i F F U A K F U K l α l α− + = + ∑
> = +∑
= (3-33) (3-32)與(3-33)所指的l
i−與l
i+分別是代表負的特徵值與正的特徵值，接著再利用平均的 方法將 1 2 i F + 更進一步表示成： ( ) 1 1 2 1 ( ) ( ) 2 m i R L i i i i F F U F U l α K + = = + − ∑
(3-34) 再由(3-22)式可再次改變 1 2 i F + 的形式如下：30 1 2 1 ( ) ( ) 2 R L i F F U F U A U + = + − ∆ (3-35) 其中∆ =U UR −UL、 1 A = A+−A−=K Λ K− ， 1
0
0
m l l Λ = 。 接下來需找出 A中所需的物理量，必須利用下列方法： 現考慮一維等溫尤拉方程式： ( ) 0 t x U +F U = (3-36) 其中 1 2 u U u u ρ ρ = = ； 1 2 2 2 f u F f u a p ρ ρ = = + ，a為聲速 方程式(3-36)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示： 2 2 0 1 ( ) 2 F A U a u u U ∂ = = − ∂ (3-37) 特徵值：l1 = − ，u a l = + 2 u a 特徵向量： (1) 1 K u a = − ， (2) 1 K u a = + 接著選定 parameter vector Q 1 2 q U Q q u ρ ρ ρ = = (3-38) 再將F與U 利用Q表示： 2 1 1 1 2 1 2 u q U q Q u q q = = = (3-39) 1 1 2 2 2 2 2 2 1 f q q F f q a q = = + (3-40) 為了表示出∆U與∆F 需在定義 averaged vectorQ：31 1 2 1 1 ( ) 2 2 L R L R L L R R q Q Q Q q u u ρ ρ ρ ρ + = = + = + (3-41) 再找出B =B Q( ) 與C=C Q ( )使得 U B Q ∆ = ∆ ；∆ = ∆F C Q (3-42) 將(3-42)結合可得 1 ( ) F CB− U ∆ = ∆ (3-43) 再根據上述條件 3 求出近似 Jcaobian 矩陣 1 A =CB − (3-44) 為了滿足(3-42)，可以求得 1 2 1 2q 0 B q q = ； 2 1 2 2 1 2 2 q q C a q q = (3-45) 再帶入(3-44)可得 2 2 0 1 2 A a u u = − (3-46)
u為 Roe averaged velocity L L R R L R u u u ρ ρ ρ ρ + = + (3-47) 因此可以用同樣方法得到以下物理量： L L R R L R v v v ρ ρ ρ ρ + = + (3-48) L L R R L R w w w ρ ρ ρ ρ + = + (3-49) L L R R L R H H H ρ ρ ρ ρ + = + (3-50) 1/ 2 [( 1)( 1/ 2 )] a= γ − H − V (3-51)
32
其中u、v、w 分別代表x方向、y方向、z方向的速度。 H 、a則分別為焓和音速。
33
y
x<0 x=0 x>0 x
圖 3-7 黎曼問題特徵值結構圖
34
3-4、Monotonic Upstream-Centered Scheme for Conservation Laws(MUSCL)：
本論文使用的是採用I. Abalakin等[19]中所使用的插分法。其方程式如下： 1/ 2 1/ 2 1/ 2 L L i i i u+ = +u ∆u+ (3-52) 1/ 2 1/ 2 1/ 2 R R i i i u+ = −u ∆u+ (3-53) 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-54) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i i u+ β u+ u β u+ u+ ∆ = − − + − +θc(−ui−1+3ui−3ui+1+ui+2) 1 2 3 ( 3 3 ) d i i i i u u u u θ + + + + − + − + (3-55) 其中(3-54)、(3-55)式中的β 、θ 、c θ 值可由表(3-1)中查得。代入不同的值可以得到d 不同的精度。本論文則是使用三階精度，以減少數值計算的消散性。 在程式中，高次項的插分法在不連續的情況下，容易使震盪變大，為了降低震盪，本 研究在 MUSCL 法插分出來的方程式中加入 minmod limiter，用來確保程式不會發散。 因此(3-52)與(3-53)式需改寫如下： 1/ 2 1/ 2 min mod( 1/ 2) L L i i i u+ = +u ∆u+ (3-56) 1/ 2 1/ 2 min mod( 1/ 2) R R i i i u+ = −u ∆u+ (3-57)
35 表 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 2 3 4 5
36 3-5、Preconditioning 法： 為了提高程式的應用範圍，在 Navier-Stokes 方程式中加入 preconditioning 法，讓程式 不論在高速或低速流內計算可壓縮流，皆可獲得精確的結果。 S U F G H t x y z ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-58) 上式為原始方程式，接著將保守形式（conserved variables）轉變成主要變數形式 （primitive variables），其形式如下： S p U F G H M t x y z ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-59) 其中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-60) 其中 p p ρ ρ =∂ ∂ ； T T ρ ρ =∂ ∂ 接著將(3-59)式的方程式乘上矩陣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-61) 再將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-62)
37 將(3-62)式帶入(3-59)式，連續方程式： ( ) S p p u v w t x y z ρ ρ ρ ρ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-63) 在理想氣體中可將(3-55)再表示成 2( ) S p u v w C t x y z γ ∂ +∂ρ +∂ρ +∂ρ = ∂ ∂ ∂ ∂ (3-64) 其中C為聲速 由(3-64)式可以看出，在等密度條件下，由於ρ 為零，(3-55)式將變成 p S u v w x y z ρ ρ ρ ∂ +∂ +∂ = ∂ ∂ ∂ (3-65) 上式即為不可壓流的連續方程式。 綜上所述，可以得知只要改變(3-62)式中的ρ 項，利用當地流場速度（local velocity）p 的倒數取代，即可轉換系統中的特徵值，藉此改變低速情況下流場的聲速，使聲速與 流場速度冪次級數(order)相同，系統不再受到 CFL(Courant-Friedrichs-Lewy Condition) 條件的限制，提高程式的計算效率。 利用θ取代ρ 項： p 2 1 1 ( ) r p U TC θ = − (3-66) max if if if r U u C U u C u C C u C ε ε ε × < × = × < < > (3-67) 其中
ε
為一極小的值，約等於 5 10− ，其主要是用來防止停滯點（stagnation point）在 計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，U 必須大於流體r的當地擴散速度（local diffusion velocity），因此U 還需加入下列限制： r
max( , ) r r U U x ν = ∆ 將θ帶入(3-62)式後，可得到一新矩陣Γ nc
38 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T nc p C θ ρ ρ ρ ρ ρ Γ = − (3-68) 經過上述推導之後，方程式從(3-59)式轉變如下： ( ) S p nc U F G H K t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-69) 為了讓(3-69)式中的通量項再度轉換成保守形式，在乘上 1 K− 1 (K nc) Up F G H S t x y z −Γ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-70) 根據(3-70)式，定義 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-71)
其中
Γ
為 preconditioning 矩陣，Up為 primitive form[ , , , , ]TP u v w T 由於方程式在時間項經過改變，因此必須重新推導 Roe 所提出的近似黎曼解。在(3-35) 式中，可以觀察到 1 2 i F + 項，是由 1 ( ( ) ( )) 2 F UR +F UL 的中央差分法加上為了解決不連續面
問題的 artificial viscosity term 1
39
artificial viscosity term 做改變即可，其推導如下：
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-72) 其中 p U M U ∂ = ∂
所以 artificial viscosity terms 改寫如下：
1 1 2 1 1 ( ) 2 R L 2 P i F F F − AM U + = + − Γ ∆ (3-73) 其中 1 1 AM KA DA KA − − Γ = × ×
40 3-6、黏滯性項之差分: 黏滯性項；在黏滯性項方面，採用二階中央差分法。由於在尤拉方程式中計算的範圍 皆為網格與網格之間的通量項，因此在黏滯性項方面，所需要得到的速度梯度項也必 須是網格之間的通量項。下列以三維的 X 方向為例，圖 3-8 為其示意圖。 圖 3-8 中各編號所代表的位置分別為： 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-74) X V V X V X V ∆ − = ∆ ∆ = ∂ ∂ (9) (7) (3-75) X W W X W X W ∆ − = ∆ ∆ = ∂ ∂ (9) (7) (3-76) Y U U Y U Y U ∆ − = ∆ ∆ = ∂ ∂ 2 ) 14 ( ) 2 ( 2 (3-77) 其中 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-78) 同理 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-79)
41 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-80) Z U U Z U Z U ∆ − = ∆ ∆ = ∂ ∂ 2 ) 5 ( ) 11 ( 2 (3-81) 其中 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-82) 同理 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-83) 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-84)
42
![圖 4-1(a) 和 Guzman[26]結果之比對圖 左:Guzman 右:本研究之結果](https://thumb-ap.123doks.com/thumbv2/9libinfo/8036830.161691/70.892.135.715.219.880/圖41a和Guzman26結果之比對圖左Guzman右本研究之結果.webp)
![圖 4-2 和 Bhavnani[2]實驗對照所設定之物理模式](https://thumb-ap.123doks.com/thumbv2/9libinfo/8036830.161691/72.892.147.773.318.821/圖42和Bhavnani2實驗對照所設定之物理模式.webp)
