國立交通大學
機械工程學系
博士論文
管道熱傳增益之動態分析研究
Investigation of dynamic analysis of the heat transfer
enhancement in the channel flow
研 究 生：黃玠超
指導教授：傅武雄 博士
管道熱傳增益之動態分析研究 研究生：黃玠超 指導教授：傅武雄 博士 國立交通大學機械工程學系
摘要
隨著工程上發展的需要，高溫壁面熱傳效應對機械的影響越益受到重視，而 元件其性能及可靠度的提升將取決於其散熱技術，為了要達到最佳的熱傳增益， 瞭解並分析高溫壁面其對流熱傳的問題，勢必為必要的研究課題。因此本計畫首 先研究關於一具有高溫壁面的三維管道可壓縮流，在其高溫壁面上加裝一移動薄 板，探討薄板往復運動於強制對流及混合對流下的現象，造成高溫壁面上溫度與 速度邊界層不斷的被移除與重新生成，因此可大幅的提升熱傳效率。此外，並研 究關於三維垂直管道自然對流的問題，探討在高雷利數時從層流轉變成不穩定之 現象，以及對於高溫壁面熱傳效率的影響，之後藉此基礎分析研究三維管道純自 然對流誘導衍生成紊流之問題，以期能提升熱傳效率，提高對實際工業應用之範 圍。Investigation of dynamic analysis of the heat transfer enhancement in the channel flow
Student：Jieh-Chau Huang Advisor：Wu-Shung Fu Department of Mechanical Engineering
National Chiao Tung University Abstract
With the improvement of the technology, the effect of heat transfer of a heat surface is getting more seriously concern in industry. The promotion of reliability and performance is dependent on the technology of heat transfer enhancement. So the problem of heat transfer of a heat surface in a convection flow is a very important issue. First, this study investigates the problem of three-dimensional compressible convection flow with insertion of a moving block numerically. The phenomena of interaction between the moving block and channel flow are observed. The enhancement of heat transfer is remarkable due to destroying the boundary layers by the moving block. Besides, the varied processes of natural convection in three-dimensional vertical channels of single heated wall with different aspect ratios and Rayleigh numbers from stable to unstable situation are investigated numerically. The phenomena from laminar to unstable transition and the relationship between Rayleigh number and heat transfer rate are observed. The results are expected to be useful for the future study about turbulence flow induced by pure natural convection. Therefore, this study is a foundation potentially for the heat transfer enhancement of practical industry.
誌謝 從大學、碩士班一直到博士班，在交大的十一年間，度過了非常充實且多采 多姿的生活。由衷地感謝指導老師傅武雄教授在這些年來給予課業和論文上的指 導，以及在生活各方面上的關心與照顧，在此謹致最高的謝忱與敬意。同時也感 謝機械系諸位師長在課業方面的指導。另外要特別感謝博士班學長連信宏與李崇 綱，在研究計算上的協助以及指導，還有實驗室學弟妹們在各方面的協助和精神 上的鼓勵讓我可以順利完成論文。更要感謝父母親在人生道路上一路的支持與幫 助，無論在心靈或是物質上的投資，從不間斷地給予耐心與支持，讓我毫無後顧 之憂專心在學業上。最後感謝朋友、同學以及親友們的鼓勵與關心，今日才能順 利完成學業。
目錄 摘要... i 英文摘要... ii 誌謝... iii 目錄... iv 表目錄... vi 圖目錄... vii 符號表... x 一、緒論... 1 二、物理模式... 7 2-1 三維水平管道加裝移動平板強制對流物理模式 ... 7 2-2 三維垂直管道加裝移動平板混合對流物理模式 ... 11 2-3 三維垂直管道自然對流物理模式 ... 15 三、數值方法... 18 3-1 統御方程式 ... 19 3-2 Roe scheme ... 22 3-3 MUSCL 法 ... 29 3-4 Preconditioning 法 ... 31 3-5 LUSGS 法... 39 3-6 非反射性邊界 ... 41 3-7 座標轉換 ... 45 3-8 沉浸邊界法 ... 48 四、結果與討論... 50 4-1 三維管道加裝移動平板強制對流 ... 50 4-2 三維垂直管道加裝移動平板混合對流 ... 82
4-3 三維垂直管道自然對流 ... 102 五、結論... 119 參考文獻... 121
表目錄 表 3-1 精度係數值 ... 30 表 4-1 本研究模擬之兩個不同雷諾數下各兩種薄板速度的案例 ... 56 表 4-2 不同雷諾數與薄板移動之熱傳增益比較 ... 81 表 4-3 本研究在 T 12K的模擬案例 ... 86 表 4-4 本研究在 T 100K的模擬案例 ... 86 表 4-5 不同浮慣比與薄板移動之熱傳增益比較 En ( T 100K) ... 101 表 4-6 不同浮慣比與薄板移動之熱傳增益比較 En ( T 12K) ... 101 表 4-7 三維垂直管道自然對流案例 ... 106
圖目錄 圖 2-1 三維水平管道加裝移動平板物理模式 ... 9 圖 2-2 鄰近邊界網格示意 ... 10 圖 2-3 三維垂直管道加裝移動平板物理模式 ... 14 圖 2-4 三維垂直管道自然對流物理模式 ... 17 圖 3-1 黎曼問題特徵值結構圖 ... 28 圖 3-2 差分示意圖 ... 36 圖 3-3L 、1 L 、2 L 、3 L 與4 L 於管道兩端的方向性示意圖 ... 44 5 圖 4-1 網格測試圖 ... 57 圖 4-2 解析解與數值解之比較 ... 58 圖 4-3 雷諾數 200 平均紐塞數隨時間周期變化示意圖 ... 59 圖 4-4 雷諾數 600 平均紐塞數隨時間周期變化示意圖 ... 60 圖 4-5 x2x3平面流場示意圖.(雷諾數=200, V =1/4) ... 61 b 圖 4-6 x2x3平面溫度場示意圖.(雷諾數=200, V =1/4) ... 62 b 圖 4-7 局部紐塞數示意圖.(雷諾數=200, V =1/4) ... 63 b 圖 4-8 x2x3平面流場示意圖.(雷諾數=200, V =1/4) ... 64 b 圖 4-9 x2x3平面溫度場示意圖.(雷諾數=200, V =1/4) ... 65 b 圖 4-10 局部紐塞數示意圖.(雷諾數=200, V =1/4) ... 66 b 圖 4-11 x2x3平面流場示意圖.(雷諾數=200, V =1/4) ... 67 b 圖 4-12 x2x3平面溫度場示意圖.(雷諾數=200, V =1/4) ... 68 b 圖 4-13 局部紐塞數示意圖.(雷諾數=200, V =1/4) ... 69 b 圖 4-14 x2x3平面流場示意圖.(雷諾數=600, V =1/2) ... 70 b
圖 4-15 x2x3平面溫度場示意圖.(雷諾數=600, V =1/2) ... 71 b 圖 4-16 局部紐塞數示意圖.(雷諾數=600, V =1/2) ... 72 b 圖 4-17 x2x3平面流場示意圖.(雷諾數=600, V =1/2) ... 73 b 圖 4-18 x2x3平面溫度場示意圖.(雷諾數=600, V =1/2) ... 74 b 圖 4-19 局部紐塞數示意圖.(雷諾數=600, V =1/2) ... 75 b 圖 4-20 x2x3平面流場示意圖.(雷諾數=600, V =1/2) ... 76 b 圖 4-21 x2x3平面溫度場示意圖.(雷諾數=600, V =1/2) ... 77 b 圖 4-22 局部紐塞數示意圖.(雷諾數=600, V =1/2) ... 78 b 圖 4-23 雷諾數 200 下流動方向加熱壁面中央處之平均局部紐塞數比較圖 ... 79 圖 4-24 雷諾數 600 下流動方向加熱壁面中央處之平均局部紐塞數比較圖 ... 80 圖 4-25 不同浮慣比下的管道中央 x1x2平面溫度場比較圖 ... 87 圖 4-26 x2x3平面流場示意圖(Gr/ Re2 0.11 , T = 12K ) ... 88 圖 4-27 x2x3平面溫度場示意圖(Gr/ Re2 0.11 , T = 12K ) ... 89 圖 4-28 局部紐塞數示意圖( 2 / Re 0.11 , T = 12K Gr ) ... 90 圖 4-29 x2x3平面流場示意圖(Gr/ Re2 0.11 , T = 12K ) ... 91 圖 4-30 x2x3平面溫度場示意圖(Gr/ Re2 0.11 , T = 12K ) ... 92 圖 4-31 局部紐塞數示意圖( 2 / Re 0.11 , T = 12K Gr ) ... 93 圖 4-32 x2x3平面流場示意圖(Gr/ Re2 11.5 , T = 12K ) ... 94 圖 4-33 x2x3平面溫度場示意圖(Gr/ Re2 11.5, T = 12K ) ... 95 圖 4-34 局部紐塞數示意圖( 2 / Re 11.5 , T = 12K Gr ) ... 96 圖 4-35 x2x3平面流場示意圖(Gr/ Re2 11.5 , T = 12K ) ... 97 圖 4-36 x2x3平面溫度場示意圖(Gr/ Re2 11.5, T = 12K ) ... 98 圖 4-37 局部紐塞數示意圖( 2 / Re 11.5 , T = 12K Gr ) ... 99
圖 4-38 x1流動方向加熱壁面中央處之平均局部紐塞數比較圖 ... 100 圖 4-39 R=3.75、 8 1.35 10 x Ra 、t=6.21 秒時流線圖與等溫線分布圖 ... 107 圖 4-40 R=3.75、 8 1.35 10 x Ra 、t=6.27 秒時流線圖與等溫線分布圖 ... 108 圖 4-41 R=3.75、 8 1.35 10 x Ra 、t=6.35 秒時流線圖與等溫線分布圖 ... 109 圖 4-42 R=3.75、 8 1.35 10 x Ra 、t=6.4 秒時流線圖與等溫線分布圖 ... 110 圖 4-43 R=3.75、 8 1.35 10 x Ra 、t=6.45 秒時流線圖與等溫線分布圖 ... 111 圖 4-44 R=3.75、 8 1.35 10 x Ra 加熱壁面上 ABC 三點隨時間變化之局部紐塞數 ... 112 圖 4-45 R=2.75、 7 5.35 10 x Ra 、t=6.15 秒時流線圖與等溫線分布圖... 113 圖 4-46 R=2.75、 7 5.35 10 x Ra 、t=6.2 秒時流線圖與等溫線分布圖... 114 圖 4-47 R=2.75、 7 5.35 10 x Ra 、t=6.25 秒時流線圖與等溫線分布圖... 115 圖 4-48 R=2.75、 7 5.35 10 x Ra 、t=6.3 秒時流線圖與等溫線分布圖... 116 圖 4-49 R=2.75、 7 5.35 10 x Ra 加熱壁面上 ABC 三點隨時間變化之局部紐塞數 ... 117 圖 4-50 R=2.75、 7 1.34 10 x Ra 流線圖與等溫線分布圖 ... 118 圖 4-51 R=2.75、 7 1.34 10 x Ra 加熱壁面上 ABC 三點隨時間變化之局部紐塞數 ... 119 圖 4-52 三個模擬案例隨時間變化平均紐塞數比較圖 ... 120 圖 4-53 等效雷利數 Ra*與平均紐塞數之關係圖 ... 121
Nomenclature
1
w dimensional width of the channel (m )
2
w dimensional height of the channel (m )
3
w dimensional width of the slender block (m )
1
l dimensional length of the channel (m )
2
l dimensional length of the slender block (m )
3
l the distance between the inlet and the heat surface (m )
4
l the distance between the outlet and the heat surface (m ) h dimensional height of the slender block (m )
g acceleration of gravity ( 2
/
m s )
k thermal conductivity (W mK/ )
M local Mach number
P pressure (Pa) 0 P surrounding pressure (Pa) Re Reynold number R Aspect Ration Gr Grashof number R gas constant (J kg K ) / / Ra Rayleigh number
Rax Rayleigh number defined in Eq. (4-16)
Ra* Rayleigh number defined in Eq. (4-17)
t time ( s )
t* Dimensionless time
0
T temperature of surroundings ( K )
h
T temperature of heat surface ( K ) , ,
x y z Cartesian coordinates (m )
, ,
X Y Z dimensionless Cartesian coordinates , ,
u v w velocities in
1
x , x and 2 x directions (3 m s/ )
b
v moving velocity of the slender block (m s/ ) , ,
U V W dimensionless velocities in X , Y and Z directions
b
V dimensionless velocity of the slender block H dimensionless height of the slender block
X
Nu local Nusselt number defined in Eq. (4-17) 1
X
Nu local Nusselt number defined in Eq. (4-1)
Nu average Nusselt number defined in Eq. (4-3) 1
X
Nu average Nusselt number defined in Eq. (4-4) (Nu ) average Nusselt number defined in Eq. (4-5) (Nu)t average Nusselt number defined in Eq. (4-12)
Greek symbols density ( 3 / kg m ) 0 surrounding density (kg m ) / 3 viscosity ( 2 / N s m ) 0 Surrounding viscosity (N s m / 2)
specific heat ratio
dimensionless temperature
, ,
第一章 緒論
近年來為因應許多工程上的需要，機械元件所需承受的熱負載不斷增加，造 成元件損害與故障，故高溫壁面的熱傳效應問題及對機構的影響也一再被研究。 此外隨著消費性電子產品蓬勃發展，高功率電子元件的發展與半導體製程技術的 進步，目前許多電子產品皆走向高效能、微小化的趨勢，致使電子元件舉凡電腦 CPU，顯示卡，筆記型電腦或智慧型手機等其封裝元件之發熱密度愈來愈高，其 單位熱通量亦相對地不斷增加，對性能及可靠度等方面均造成不容忽視之影響。 因此在未來電子產品的發展趨勢走向更輕薄短小之際，其性能及可靠度的提升將 取決於其散熱技術。而為了要達到最佳的熱傳增益，瞭解並分析高溫壁面其對流 熱傳的問題，勢必為必要的研究課題。一般正常電腦 CPU 在執行程式下溫度大 約為攝氏 50~60 度之間，超過攝氏八十度則電腦會基於保護狀態，自動關閉系統， 長期下來內部零組件壽命將會減短。因此為了維持元件於額定溫度下運作，必須 將此密集的熱量能有效散逸於系統外之環境。傳統的散熱片大多數是由鋁合金所 製造，其熱傳導性只屬於中等程度，對於目前元件發熱功率越來越高的情況，逐 漸有捉襟見肘的情形發生，目前電子元件的發熱量達到每平方公分數十瓦的等級， 且接點可承受的溫度約在攝氏 150 度以下。因此如何改良空氣對流的方式以及 增加對流所能散逸之熱量成為研究的主要課題之一。 混合對流的研究種類以浮力方向及強制對流流場之速度方向亦分為三種形 式。一種為浮力方向與流場之速度方向垂直的交錯流(Cross flow)，Najam 等人[1] 以數值方法研究二維水平管道混合對流，在管道中佈上加熱塊體，藉以探討不同 雷利數、雷諾數及加熱體高度對流場之影響。Kitamura 等人[2]以實驗探討加熱 圓柱誘導對流之現象，以及紐塞數和分離點位置與加熱體熱通量的關係。而另一 種浮力方向與流場之速度方向相同稱為助流(Aiding flow)，Gau 等人[3]以實驗研 究垂直管道中的浮力助流與熱傳現象，觀察到迴流現象及浮力影響參數等關係。 Desrayaud 等人[4]研究層流混合對流在浮慣比大於 1 的現象，發現隨著雷諾數的上升其迴流現象會越不明顯。第三種為浮力方向與流場之速度方向相反的逆流 (Opposing flow)，Martinez-Suastegui[5]等人以數值方法研究暫態二維垂直管道層 流混合對流之現象，其結果分析了浮慣比與雷諾數對流場之影響及紐塞數分布。 Joye 等人[6]在一垂直圓管實驗中，觀察浮力參數與熱傳增益之相互關係，以及 雷諾數與紐塞數的關係。 過 去 有 多 篇 文 獻 提 出 各 種 不 同 的 方 法 來 增 加 高 溫 壁 面 的 熱 傳 效 能 。 Bergles[7,8]將常見的增加熱傳效率方式，詳細地討論與整理並將其略分為兩大類： 一為不需另外作功的被動方式，另一類為需外加能量的主動式方法。Hwang 與 Liou[9]以實驗的方法探討在壁面加裝有縫隙的肋狀紊流產生器，對方形管道流的 熱傳與熱阻的影響。結果顯示縫隙狀肋狀紊流產生器可避免壁面產生熱點，同時 對熱傳效能有所增益。Park 等[10]研究不同形狀之熱交換管內插入結構對提高熱 傳效率的影響。結果顯示圓錐形盤管與肋狀的插入結構約可提高 30％的熱傳效 率。Iida 等[11]在液體中加入鋁金屬微粒，當液體發生沸騰現象時，鋁金屬微粒 也被液化，可將熱傳效率提高十倍。Fusegi[12]數值分析有肋狀突起管道中之振 動對向流，當提高振動對黏性的相對強度時，可顯著的提高熱傳效率。Sitter 等 [13]實驗探討高強度音場對重力場與微重力場下池沸騰的影響，發現音場有助於 提高熱傳效率。Liu[14] 在通道內設凸起物作往復運動時混合對流之熱傳研究， 可得知於往復設凸起物的垂直通道中往復運動可使流場的熱傳效率提昇至靜止 通道的 2~3 倍。Wang[15] 在其研究中探討了往復運動的薄塊對噴流中壁面熱傳 影響之研究，可發現往復運動的薄塊對於壁面邊界層不斷的破壞與生成，達成提 升熱傳效率的目的。Wu[16]以數值分析研究二維含加熱方塊之管道中，置入薄 板最多可提昇 39.5%的紐塞數。Dogan [17]利用實驗探討混合對流中，在水平矩 形管道內設置散熱片，研究在不同散熱片高度及散熱片間隔對於熱傳效果之影響， 可得知在散熱片高度越高，可增加熱傳面積，且不易使散熱片內較高溫流體從尖 端流進另一片散熱片，達到增加熱傳之效果。楊[18]探討在二維、不可壓縮流中， 混合對流對傾斜平板管道內設置橫向散熱薄板熱傳之研究，由研究可得知，在最
理想的長寬比下，熱傳量伴隨著雷諾數之增加。Hamouche[19]探討二維、不可壓 縮混合對流中，在水平管道內設置突出的熱源方塊，探討空氣冷卻之效果，由研 究可得知，增加分離區之距離可增加熱傳量。Dogan[20]探討在混合對流中，水 平管道內設置分離的多個熱源在管道上方及下方之熱傳研究，可得知將電子元件 設置在底部最前及最後列可得到最大的熱傳效率。Youssef[21]討論在層流混合對 流中，浮力向上或向下對熱傳及質傳的影響，由研究中可得知，在重力與流動方 向相反的流動中，存在壁面與進口流體的高溫差，會產生流體不穩定性及紊流產 生。 比較上述各種增加熱傳量的方法，不論是被動式或是主動式方法增加熱傳效 率似乎都會受到限制。主要的原因是流場在熱傳面上生成的速度邊界層與溫度邊 界層阻礙了熱量傳出。當流體流經高溫壁面時，壁面上的邊界層有速度邊界層與 溫度邊界層，所謂速度邊界層，是指當連續體之黏性流體以一接近速度流經固體 物體時，在雷諾數較高的情況下，流體在固體表面附近受到黏性的影響，形成速 度變化非常激烈且厚度極薄的區域，在此區域內，流體速度從壁面的零漸次變化 到此區域的邊緣速度。而所謂溫度邊界層則是流體受壁面高溫影響，在壁面產生 熱傳，其靠近壁面的溫度逐漸變化至與壁面溫度相同的區域。在溫度邊界層中， 流體的溫度變化和緩，溫度梯度較小。根據熱傳導之傅利葉定律，在溫度梯度較 小的情況下，熱傳量也較小。因此，溫度邊界層的存在將限制壁面高溫所能傳出 之熱傳量。而為了大幅提高高溫壁面的熱傳效率，就必須移除壁面上的溫度邊界 層，使得高溫壁面與低溫流體直接接觸，進而增加高溫壁面與低溫流體的溫度梯 度，達成增加壁面熱傳量的目的。 本研究的工作流體將考慮可壓縮流，因此密度隨著壓力與溫度的變化而改變， 對於實際工業界中的應用較為廣泛與實際。而在薄版移動的同時，流體受到薄版 的牽引，使溫度與速度場隨薄版移動而產生變化，流體與移動物體相互影響的問 題，屬於動態的移動邊界問題，又因為在處理實際散熱問題時，工作流體多為可 壓縮流，因此在本研究中的工作流體不論高速或低速皆視為可壓縮流，以增加應
用範圍。Fu[22]發展出黏性流場之全域速度場數值解法，此種方法最大的困難處 在於計算低速流場時，由於可壓縮流必須遵守 Courant-Friedrichs-Lewy (CFL)條 件，因此在低速流體時，受限於流體變化傳遞速度(約等於聲速)，時階將會極小。 在此種情況造成計算過程將耗費極大的時間與整體過低的計算效率。為了改善此 缺點，本研究在計算可壓縮流時加入 Preconditioning 法，藉此讓流體即使在低 速時，也可有較高的效率與良好的收斂性。在計算此種低馬赫數流體的方面，目 前有密度基底法(Density-based method)與壓力基底法(Pressure-based method)。本 研究以密度基底法為主。而在密度基底法中又以 Turkel[23]提出 Preconditioning 法最為廣泛應用，不僅可同時應用在可壓縮流與不可壓縮流中，更可以讓程式的 收斂性增加。而本研究在做數值計算時，主要是利用到網格之間的物理量，因此 用到 MUSCL 法來處理網格間的物理量，在計算流場時空間部分以 ROE 法[24] 來計算出非黏滯性項的通量，黏滯性項部分則採用二階中央插分，重力項則不使 用 Boussinesq approximation，由 Gray[25]可得知，在冷熱表面溫差小於 30K 時才 可適用，但在現實工業應用中，例如：在半導體製程，蒸鍍、乾燥過程中，往往 溫差高達上百度，故較不適用；在時階部份為了能夠使程式能夠加速收斂，因此 採用 LUSGS implicit[26]方法；而當在程式中加入 Preconditioning 法時，同時破 壞了統御方程式，為了彌補此一缺點，必須加入 Artificial time term 來修正方程 式；最後利用 Dual time stepping[27]來計算暫態的物理量。在靠近壁面的地方， 為了增加程式的收斂效果以及更清楚的觀察壁面效應，對於壁面網格做加密的處 理；而在處理移動邊界的問題時，本研究採用 Immersed boundary[28]中的移動邊 界方法，此方法在處理複雜外型及移動邊界的問題十分的有效率。 此外，在處理複雜的流體力學問題時需要大量的計算過程，如何利用多核心 處理器來提升運算速度已為目前的發展主流。多核心處理器對於單一執行緒在平 行運算方面，並無法提升計算速度；若利用多執行緒的程式架構，可透過不同核 心來同時計算，達到提升計算效率、節省計算時間之目的；但是多執行緒的程式 在撰寫、編輯上，也都比單一執行緒的程式架構要複雜。本文採用 OpenMP 方
法來提升運算速度。常見的平行運算方法有 MPI 和 OpenMP 兩種方法；MPI (Message Passing Interface)是一種分散式記憶體(distributing memory)的觀念，而 OpenMP (Open Multi-Processing)為共享式記憶體(sharing memory)的觀念，而兩種 方法各有優缺點。MPI 在撰寫程式上面相較於 OpenMP 較容易，但計算速度較 慢且會受限於網路效率，在設備擴充方面較便宜；反之 OpenMP 在程式撰寫上 較困難，但計算速度快且不會受限於網路效率，在設備擴充方面需較大的經費。 因此本程式利用 OpenMP 來進行平行化運算，效率為原程式的四倍，減少計算 時間及成本。 故本研究首先研究關於一具有高溫壁面的三維管道流，在其高溫壁面上加裝 一薄平板，此薄平板直接與高溫壁面接觸，並且在高溫壁面上做與主流方向垂直 之往復運動。受到薄平板移動的影響，薄平板前進方向前的流體受到推擠而往薄 板移動方向移動，使的薄板前方的溫度與速度邊界層被破壞。而在薄板後方的流 體，則受到薄板的牽引而向薄板流動，在填補因薄版移動所產生的空間的同時， 低溫的流體將與高溫的壁面直接接觸，因而產生最大的溫度梯度，並且同時在壁 面上生成新的溫度與速度邊界層，熱傳量因此增加。由於薄版在高溫壁面上不斷 的做往復運動，高溫壁面上的溫度與速度邊界層不斷的被移除與重新生成，因此 可大幅的提升熱傳效率，相對於過往文獻中的增加熱傳方法大都是增加噴流上的 熱傳量，此種利用往復薄板增加熱傳的方法更能增加下游區的熱傳。 此外，本研究將考慮有限長度管道自然對流衍生不穩定之熱傳問題，此類研 究在現今文獻中少有論述，過去大都為討論強制對流下的紊流模擬，由於為強制 對流，故能一開始得知流量而給予特定擾動進行循環計算，其結果只能類似於穩 態紊流流場，然而對於自然對流慢慢轉變不穩定進而成紊流之過程，無法預先得 知由於氣體受熱其密度與壓力變化所導致的流量，亦即無法確定誘導成紊流的驅 動力，即無法使用前述方法計算。過去常以擾動流場利用紊流現象增加熱傳效益， 而大多的數值模擬研究皆是以周期性邊界加以些許紊流強度反覆計算得到一個 紊流流場，並非真實物理尺寸下的紊流模擬，若考慮自然對流狀態下的物理模型，
無法預先得知因流場變化而產生的進口流量，此種方法即無法模擬出自然對流從 層流轉變成紊流之過程，故研究此類自然對流誘導紊流之研究亦為重要方向之一。 為了提高實際應用範圍，為了瞭解上述自然對流誘導不穩定之熱傳問題，將研究 模擬垂直管道自然對流誘導紊流的問題，並對其從層流轉變成不穩定過程中的機 制加以分析，以期日後能了解紊流流場對熱傳的影響並加以改善。
第二章 物理模式
2-1、三維水平管道加裝移動平板強制對流物理模式： 圖2-1為一水平管道在底面部份加熱管道中加裝一往復移動薄板之物理模式 圖。一長度為l1，寬度為w1，高度為w2的方形管道，溫度T0的冷卻空氣以等速 度u0由管道入口噴入，在管道的底部ABDC部份區域為高溫壁面，對管道內的空 氣加熱，進行熱交換；在高溫壁面上有一長度為l2，寬度為w3，高度為h的絕熱 薄板。當時間為零時，薄板靜止在加熱底面的3/4位置(pr)，也就是薄板做往復 運動的起點，待流場穩定後，薄板開始以小於進口速度的移動速度vb往加熱底面 1/4的位置(pl)，也就是薄板移動的終點移動。當薄板移動又回到最初的起始位置 時即完成一週期的運動。本研究即當薄板在管道內做往復運動時，對管道內的熱 場及流場進行分析。x1方向為streamwise 的方向、x2為垂直壁面方向，x3為 spanwise 的方向。u1、u2、u3分別為其所對應的速度。而為了有效的增進熱傳 效率以及讓管道內的流場更容易成為完全發展流，在管道的長度以及寬度上需要 一定的比例，進口速度與薄板移動速度也有一定的關係，本研究選定管道長度為 寬度的15倍長，以確保其出口流體為完全發展流形式；而薄板的移動速度則為進 口速度的1/4、1/2，而薄板的移動速度不可大於入口流體速度，這是為了要使入 口進來之流體能順利往出口流動，並且達到完全發展流的條件。 本研究選擇層流流場作為模擬流場，流場作以下假設： 1. 可壓縮流，空氣密度會隨溫度與壓力而改變。 2. 工作流體為空氣，假設為理想氣體。流體性質為牛頓流體(Newtonian fluid)， 黏滯係數為等方向性。 3. 忽略重力效應影響。 4. 考慮溫度變化對流體所造成的影響。 此外除了加熱壁面外，其餘壁面皆考慮無滑移條件(No-slip condition)，其表示如 下:1 1 2 2 3 3 ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) P i k P i k u i k u i k u i k u i k u i k u i k T i k T i k (2-1) 而加熱壁面的邊界條件如下列所示: 1 1 2 2 3 3 ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) 2 h ( ,1, ) P i k P i k u i k u i k u i k u i k u i k u i k T i k T T i k (2-2) 其中Th為加熱壁面的溫度。如圖 2-2 所示，0 指的是 ghost cell，1 為最接近壁面 的網格。 進口等速及出口完全發展流的各表示如下: 1 0 2 3 0 (0, , ) 2 (1, , ) (2, , ) (0, , ) (0, , ) 0 (0, , ) 0 (0, , ) P j k P j k P j k u j k u u j k u j k T j k T (2-3) 1 1 2 2 3 3 ( 1, , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) atm P nx j k P u nx j k u nx j k u nx j k u nx j k u nx j k u nx j k T nx j k T nx j k (2-4)
圖 2-1 三維水平管道加裝移動平板物理模式 B A D C E F 1
w
2w
1l
2l
3l
4l
rp
lp
0u
fully-developed
flow
h 3w
0T
hT
2x
1x
3x
B A D C E F 1w
2w
1l
2l
3l
4l
rp
lp
0u
fully-developed
flow
h 3w
0T
hT
2x
1x
3x
圖 2-2. 鄰近邊界網格示意 1
0
2-2、三維垂直管道加裝移動平板混合對流物理模式： 圖 2-3 為一垂直管道在底面部份加熱管道中加裝一往復移動薄板之物理模式 圖。在管道的底部 ABDC 部份區域為高溫壁面，對管道內的空氣加熱，進行熱 交換；在高溫壁面上有一長度為l2，寬度為w3，高度為h的絕熱薄板。當時間為 零時，薄板靜止在加熱底面的 3/4 位置(pr)，也就是薄板做往復運動的起點，待 流場穩定後，薄板開始以小於進口速度的移動速度vb往加熱底面 1/4 的位置(pl)， 也就是薄板移動的終點移動。當薄板移動又回到最初的起始位置時即完成一週期 的運動。本研究即當薄板在管道內做往復運動時，對管道內的熱場及流場進行分 析。x1方向為 streamwise 的方向、x2為垂直壁面方向，x3為 spanwise 的方向。 1 u 、u2、u3分別為其所對應的速度。而為了有效的增進熱傳效率以及讓管道內 的流場更容易成為完全發展流，在管道的長度以及寬度上需要一定的比例，進口 速度與薄板移動速度也有一定的關係，本研究選定管道長度為寬度的 15 倍長， 以確保其出口流體為完全發展流形式；而薄板的移動速度則為進口速度的 1/4、 1/2，而薄板的移動速度不可大於入口流體速度，這是為了要使入口進來之流體 能順利往出口流動，並且達到完全發展流的條件。 圖 2-3 為一垂直管道在底面部份加熱管道中加裝一往復移動薄板之物理模式 圖。如圖 2-3 以及圖 2-4 所示，一長度為l1，寬度為w1，高度為w2的方形管道， 溫度T0的冷卻空氣以等速度u0由管道入口噴入，在管道的底部則部份區域為高 溫壁面 Th，對管道內的空氣加熱，進行熱交換；在高溫壁面上有一長度為l2，寬 度為 w3，高度為 h 的絕熱薄板。當時間為零時，薄板靜止在加熱底面的 3/4 位置 ( Ps )，也就是薄板做往復運動的起點，待流場穩定後，薄板開始以小於進口速度 的移動速度 vb往加熱底面 1/4 的位置(Pe )，也就是薄板移動的終點移動。當薄板 移動又回到最初的起始位置時即完成一週期的運動。本研究即當薄板在管道內做 往復運動時，對管道內的熱場及流場進行分析。x 方向為 streamwise 的方向、y 為 垂直壁面方向最後 z 為 spanwise 的方向。u 、v 、w 分別為其所對應的速度。 而為了有效的增進熱傳效率以及要熱管道內的熱流場更容易成為完全發展流，在
管道的長度以及寬度上需要一定的比例，進口速度與薄板移動速度也有一定的關 係，本研究選定管道長度為寬度的 15 倍長，以確保其出口流體為完全發展流形 式；而薄板的移動速度則為進口速度的 1/4、1/2，要注意薄板的移動速度不可大 於入口流體速度，這是為了要使入口進來之流體能順利往出口流動，並且達到完 全發展流的條件。物理尺寸比例如下： 2 1 3 1 1 1 1 2 1 3 1 4 / 1 / 0.05 / 0.5 / 15 / / 2.25 10.5 w w w w h w l w l w l w l (2-5) 本研究選擇層流流場作為模擬流場，流場作以下假設： 1. 可壓縮流，空氣密度會隨溫度與壓力而改變。 2. 工作流體為空氣，假設為理想氣體。流體性質為牛頓流體(Newtonian fluid)， 黏滯係數為等方向性。 3. 考慮溫度變化對流體所造成的影響。 此外除了加熱壁面外，其餘壁面皆考慮無滑移條件(No-slip condition)，其表示如 下: 1 1 2 2 3 3 ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) P i k P i k u i k u i k u i k u i k u i k u i k T i k T i k (2-6) 而加熱壁面的邊界條件如下列所示: 1 1 2 2 3 3 ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) ( ,1, ) ( , 0, ) 2 h ( ,1, ) P i k P i k u i k u i k u i k u i k u i k u i k T i k T T i k (2-7) 其中Th為加熱壁面的溫度。如
進口等速及出口完全發展流的各表示如下: 1 0 2 3 0 (0, , ) 2 (1, , ) (2, , ) (0, , ) (0, , ) 0 (0, , ) 0 (0, , ) P j k P j k P j k u j k u u j k u j k T j k T (2-8) 1 1 2 2 3 3 ( 1, , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) ( 1, , ) ( , , ) atm P nx j k P u nx j k u nx j k u nx j k u nx j k u nx j k u nx j k T nx j k T nx j k
2-3、三維垂直管道自然對流物理模式： 本研究為有限長度三維管道自然對流衍生不穩定之熱傳問題，共有兩種不同 長寬比之物理模型觀察其機制變化。三維垂直管道所採用的物理模式示意圖如圖 2-4 所示，主要為一個高h與寬w的三維垂直管道，z 方向採用循環性邊界。本研 究一共有兩種長寬比 R(Rh w/ )，分別為 R=2.75 和 R=3.75，其中高溫壁面位於 左側整個壁面，該一高溫壁面的溫度為Th 400K及 700K 兩種，在垂直管道中 除了高溫壁面外，另一側的壁面為絕熱狀態。而外面的環境溫度與壓力分別為 300 c T K與P0 101300Pa。 在計算時所採取的基本假設如下： 1. 流場為三維流場。 2. 工作流體為空氣，假設其為理想氣體。 3. 進出口條件皆為完全非反射條件。 求解的方程式為考慮黏滯性、流體可壓縮性與重力影響的三維 N-S 方程式與理 想氣體方程式： U F G S t x y (2-9) 與 PRT (2-10) 其中，U 、 F 、G與S的值如下所示： u U v E 2 xx xy xx xy u u P F uv T Eu Pu k u v x
2 yy yx yy v vu G v P T Ev Pv k u v y 0 0 0 ( ) 0 ( ) g S gu 此處 1 2 2 ( ) ( 1) 2 p E u v 。 空氣的黏滯係數與熱傳導係數則根據 Sutherland’s law 求得： 2 0 3 0 0 110 ( ) ( ) 110 T T T T T (2-11) ( ) ( ) ( 1) Pr T R k T (2-12) 此處的 3 0 1.1842kg m/ 、 2 9.81 / g m s、 5 2 0 1.85 10 N s m/ ，T0 298.0592K， 1.4 ，R287 /J kg K/ 與Pr0.72。
第三章 數值方法
本章主旨在說明本論文的數值計算所使用的所有模式。 第一節整理所求解 的 Navier-Stokes 方程式。將 Navier-Stokes 方程式拆解為非黏滯項與黏滯項。第 二節介紹的為黎曼解中的 ROE 法，利用 ROE 法來求出非黏滯項的通量。 接著 第三節介紹 MUSCL 法，此法是為了要解出 ROE 法中使用的網格之間的物理量， 然後為了防止在高階插分時產生震盪現象，在 MUSCL 法插分的結果方程式中 加入 Minmod limiter 以確保程式不會發散。第四節為介紹 Preconditioning 法， 因為當計算低速可壓縮流時，因速度和音速的數量級上差距過大，在數值分析時 造成計算的困難，所以為彌補此一缺點須使用 Preconditioning 法。 最後第五節 為 LUSGS Scheme，為了加快收斂速度並且避免能量耗散問題，因此利用 LUSGS Scheme，而程式因為在使用 Preconditioning 時，加入 Artificial time term 時，已 破壞了整個統御方程式，因此需使用 Dual time stepping 疊代使其在 Artificial domain 收斂時才能進入下一個真實時階，將在此小節做詳細的解說。第六節為 座標的轉換，本研究為了要觀察邊界的地方流場及熱傳的變化及邊際效應的影響， 在靠近壁面的網格部份做了加密的動作，當程式作計算時使用曲線座標系做計算， 但是當插分時是在直角座標系做插分的動作，因此對於方程式座標的轉換在此節 做說明。第七節則說明處理移動邊界的沉浸邊界法，此方法除了適當的處理流場 外型複雜度，同時於網格節點保留計算結果之準確性，因此在卡式座標中可處理 移動或變形之複雜外型流場問題，卻不需要隨時間產生網格。綜合上述，本論文 在數值上的計算過程為，首先將 Navier-Stokes 方程式拆解為非黏滯項與黏滯項。 利用 MUSCL 法算出 ROE 法所需要的網格間物理量並搭配 Minmod limiter 以 確保程式不會發散，求解出非黏滯性項的通量，並且在計算通量時加入Preconditioning 法，以拉近與音速的數量級。 接下來利用二階中央插分法對黏 滯項做插分進而求出黏滯性項；然後再與 ROE 法求出的非黏滯性項通量做結合 得到真正的物理通量。最後使用 LUSGS Scheme 疊代以求出下一時階的物理量。
3-1、統御方程式： 本研究在計算流場的方面其統御方程式分可為兩大部分，第一部份為非黏滯 性項的尤拉方程式， 第二部份為黏滯性項。 U F G H S t x y z PRT 其中在水平管道中 E w v u U xz xy xx xz xy xx w v u x T k Pu Eu uw uv P u u F 2 yz yy yx yz yy yx w v u y T k Pv Ev vw P v u v G 2 zz zy zx zz zy zx w v u z T k Pw Ew P w v u w H 2 gv g S ) ( 0 ) ( 0 0 0 0 (3-1) (3-2) (3-3) (3-4) (3-5) (3-6) (3-7)
在垂直管道中 E w v u U xz xy xx xz xy xx w v u x T k Pu Eu uw uv P u u F 2 yz yy yx yz yy yx w v u y T k Pv Ev vw P v u v G 2 zz zy zx zz zy zx w v u z T k Pw Ew P w v u w H 2 gu g S ) ( 0 0 ) ( 0 0 0 為密度，P 為壓力。u、v、w 分別為 x、y、z 方向的速度。 2 2 2 1 ( ) ( 1) 2 P E u v w ， 為密度，P 為壓力。u、v、w 分別為 x、y、z 方向的速度。，黏滯係數與熱傳 導係數 k 利用 Sutherlands’s law 110 110 ) ( ) ( 3 0 2 0 0 T T T T T (3-8) (3-9) (3-10) (3-11) (3-12) (3-13) (3-14)
Pr ) 1 ( ) ( ) ( T R T k 其中 0=1.1842kg/m3，g=9.81m/s2， 0=1.85×10-5Ns/m2，T0=298.0592K， =1.4，R=287J/kg/K ，Pr=0.72。 上式可拆解為黏滯性項與非黏滯性項： 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 左式由非黏滯項Finviscid組成的方程式即稱為尤拉方程式。 (3-15) (3-16)
3-2、Roe scheme： 在雙曲線的守恆形式方程式中，若其初始條件包含有不連續的片段連續 (piecewise)常數，此類型的問題通稱為黎曼(Riemann)問題。因為其包含有不連續 解，因此在流體計算上有著相當廣泛的應用。一維線性黎曼方程式如下： 0 U U A t x ，其中 A 為一常數 Jacobian 矩陣。 (3-17) 初始條件為 (0) (0) (0) 1 ( , , m )T U u u 。左上方括號代表時間 t 為 0。 求出 A 之特徵值矩陣以及特徵向量。 1 A K K ，其中為特徵值矩陣： 1
0
0
m 。 (1) ( ) , , m T K K K 為特徵向量，故 ( )i ( )i i AK K 。 接著定義特徵變數W〈characteristic variables〉，其定義如下： ( , ) W W t x ， ( 1) WK 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-18) 方程式(3-5)稱為 canonical form 或 characteristic form。將以上的結果簡單整理如下： 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-19) (3-11)可由特徵曲線法求得其解為： (0) ( , ) ( ) i i i w x t w xt i xit0 (0) ( , ) ( ) i i i w x t w xt i xit0 (3-20)
其中，i與i為初始值的特徵變數。由於U KW，可以得到 (0) ( ) ( , ) ( ) m i i i i U x t
w xt K 參照圖(3-1)，可以進一步推導出 ( ) ( ) 1 1 ( , ) p m i i i i i i p U x t K K
(3-21) 除此之外，還可決定出 ( , )U x t 中的 jump U： ( ) 1 m i R L i i U U U K
(3-22) 其中i i i。 在一維線性黎曼問題中，雖然有 exact solution，但在非線性問題裡需利用疊 代等方法，這些動作將耗費大量的時間，因此在實際應用上並不廣泛。為了解決 此問題，一般皆求解近似黎曼問題〈approximation Riemann problem〉解而不直 接求其 exact solution。在求解近似黎曼問題中最被廣泛應用的方法為 Roe 所提出， 亦即為 Roe scheme，其內容如下： 假設一維尤拉方程式： 0 U F t x (3-23) 根據 chain rule，可將方程式(3-23)改寫如下： 0 U F U t U x 再令 ( ) F A U U ，於是方程式(3-23)可以表示成： ( ) 0 U U A U t x (3-24) 其中， ( )A U 就稱為 Jacobian 矩陣。而 Roe scheme 將原本的 Jacobian 矩陣 ( )A U 用一常數 Jacobian 矩陣A U U( L, R)代 替，因此本來的黎曼問題可以改寫成近似黎曼問題： ( ) 0 U U A U t x ( , 0) U x UL x0 (3-25) ( , 0) U x UR x0
於是前述方法可以得到(3-23)的近似解。由以上的原理可得知，在近似黎曼 問題上，Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性 轉變成線性，但是初始條件並沒有改變，因此可以得到方程式(3-18)的近似解。 為了要求得合理的常數 Jacobian 矩陣，須合乎 Roe 所提出的四項條件： 1. U與 F 之間，存在著線性轉換的關係。 2. 當URUL U，則A U U( L, R)A U( )，此處A F U 。 3. A U( LUR)FLFR。 4. 矩陣A 的特徵向量必須線性獨立。 這四項條件都是雙曲線方程式所需具備的，這同時也說明了 Roe 所提出的 常數 Jacobian 矩陣必須有實數特徵值，其所對應的特徵向量必須線性獨立。除此 之外，條件 3.則是為了符合守恆定律（conservation law）與 Rankine-Hugoniot 條 件。 線性黎曼問題的解析解，可以直接從(3-20)與(3-21)式得到， 1 2 ( / ) i U x t 的解 可以利用下面的方程式計算： ( ) 1 0 2 ( / ) i i L i i U x t U K
(3-26) 或 ( ) 1 0 2 ( / ) i i R i i U x t U K
(3-27) 其中 1 2 i 表示網格與網格之間的交界面(face)。 而黎曼問題的近似解，則須從解近似黎曼問題著手： ( ) 0 U F U t x ，根據(3-17)式可得知F AU 為了符合守恆的條件，因此下式必須成立： ( R) ( L) ( R) ( L) F U F U F U F U (3-28) 接著在固定體積的條件下，積分近似解 1 2 (0) i U ，可得到通量(flux)的數值公式：1 1 2 2 ( (0)) ( R) ( R) i i F F U F U F U (3-29) 再從F AU的關係中可進一步求得： 1 1 2 2 (0) ( R) R i i F AU F U AU (3-30) 再根據(3-18)式與(3-19)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i R i i i i F F U A K F U K
(3-31) 或 ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i L i i i i F F U A K F U K
(3-32) (3-31)與(3-32)所指的i與i分別是代表負的特徵值與正的特徵值，接著再利用 平均的方法將 1 2 i F 更進一步表示成： ( ) 1 1 2 1 ( ) ( ) 2 m i R L i i i i F F U F U K
(3-33) 再由(3-8)式可再次改變 1 2 i F 的形式如下： 1 2 1 ( ) ( ) 2 R L i F F U F U A U (3-34) 其中 U URUL、 1 A AAK K ， 10
0
m 。 接下來需找出 A 中所需的物理量，必須利用下列方法： 現考慮一維等溫尤拉方程式： ( ) 0 t x U F U (3-35) 其中 1 2 u U u u ； 1 2 2 2 f u F f u a p ，a為聲速 方程式(3-30)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示：2 2 0 1 ( ) 2 F A U a u u U (3-36) 特徵值：1 u a，2 u a 特徵向量： (1) 1 K u a ， (2) 1 K u a 接著選定 parameter vector Q 1 2 q U Q q u (3-37) 再將 F 與U利用 Q 表示： 2 1 1 1 2 1 2 u q U q Q u q q (3-38) 1 1 2 2 2 2 2 2 1 f q q F f q a q (3-39) 為了表示出U與 F 需在定義 averaged vector Q ： 1 2 1 1 ( ) 2 2 L R L R L L R R q Q Q Q q u u (3-40) 再找出BB Q( )與CC Q( )使得 U B Q ； F C Q (3-41) 將(3-36)結合可得 1 ( ) F CB U (3-42) 再根據上述條件 3 求出近似 Jcaobian 矩陣 1 ACB (3-43) 為了滿足(3-36)，可以求得
1 2 1 2q 0 B q q ； 2 1 2 2 1 2 2 q q C a q q (3-44) 再帶入(3-30)可得 2 2 0 1 2 A a u u (3-45)
u為 Roe averaged velocity
L L R R L R u u u (3-46) 因此可以用同樣方法得到以下物理量： L L R R L R v v v (3-47) L L R R L R w w w (3-48) L L R R L R H H H (3-49) 1/ 2 [( 1)( 1/ 2 )] a H V (3-50) 其中u、v、w分別代表x方向、 y 方向、 z 方向的速度。 H 、a則分別為焓和 音速。(3-46)～(3-50)式中的UL以及UR則是利用 MUSCL 法求出。
y x0 x0 x0 x 圖 3-1 黎曼問題特徵值結構圖 p i 2 1 1 p i 1 m m
3-3、Monotonic Upstream-Centered Scheme for Conservation Laws(MUSCL)： 本論文使用的是採用I. Abalakin中所使用的插分法。其方程式如下： 1/ 2 1/ 2 1/ 2 L L i i i u u u (3-51) 1/ 2 1/ 2 1/ 2 R R i i i u u u (3-52) 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-53) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i i u u u u u c(ui13ui3ui1ui2) 1 2 3 ( 3 3 ) d i i i i u u u u (3-54) 其中(3-45)、(3-46)式中的、c、d 值可由表(3-1)中查得。代入不同的值可以 得到不同的精度。本論文則是使用三階精度，以減少數值計算的消散性。 在程式中，高次項的插分法在不連續的情況下，容易使震盪變大，為了降低震盪， 本研究在 MUSCL 法插分出來的方程式中加入 minmod limiter，用來確保程式不 會發散。 因此(3-51)與(3-52)式需改寫如下： 1/ 2 1/ 2 min mod( 1/ 2) L L i i i u u u (3-55) 1/ 2 1/ 2 min mod( 1/ 2) R R i i i u u u (3-56) min mod( , )x y Sgn x Max( ) {0,Min x[ ],ySgn x( )} (3-57)
表 3-1：精度係數值 c d Order 1/3 1/3 1/3 1/3 1/3 0 -1/6 0 -1/10 -1/10 0 0 -1/6 -1/15 -1/15 2 3 4 5 6
3-4、Preconditioning 法：
為了增加 N-S 方程式於低馬赫速可壓縮流的準確度與效率，因此於方程式 中增加 preconditioning 法。本程式採用 Weiss and Smith 的 preconditioning method [29]，撰寫於三維曲線座標，其方程式如下： U F G H S t x y z (3-58) 上式為原始方程式，接著將保守形式（conserved variables）轉變成主要變數形式 （primitive variables），其形式如下： p U F G H M S 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-60) 再將 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-61)
將(3-61)式帶入(3-59)式，連續方程式： ( ) p p u v w S t x y z (3-62) 在理想氣體中可將(3-49)再表示成 2 ( ) p u v w S C t x y z (3-63) 其中C為聲速 由(3-63)式可以看出，在等密度條件下，由於p為零，(3-62)式將變成 u v w S x y z (3-64) 上式即為不可壓流的連續方程式。 綜上所述，可以得知只要改變(3-61)式中的p項，利用當地流場速度（local velocity）的倒數取代，即可轉換系統中的特徵值，藉此改變低速情況下流場的 聲速，使聲速與流場速度冪次級數(order)相同，系統不再受到 CFL 條件的限制， 提高程式的計算效率。 利用 取代p項： 2 1 1 ( ) r p U TC (3-65) max U if u C r U u if C u C (3-66) C if u C 其中 為一極小的值，約等於 5 10，其主要是用來防止停滯點（stagnation point） 在計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，Ur必須
大於流體的當地擴散速度（local diffusion velocity），因此Ur還需加入下列限制：
max( , ) r r U U x 將 帶入(3-61)式後，可得到一新矩陣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-67) 經過上述推導之後，方程式從(3-62)式轉變如下： ( ) p nc U F G H K S t x y z (3-68) 為了讓(3-62)式中的通量項再度轉換成保守形式，在乘上 1 K 1 1 (K nc) Up F G H SK t x y z (3-69) 根據(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 (3-70) 最後方程式簡化成如下形式： p U F G H S t x y z (3-71) p U 為 primitive form[ , , , , ] /T P u v w T J ，由於方程式在時間項經過改變，因此必須 重新推導 Roe 所提出的近似黎曼解。在(3-34)式中，可以觀察到 1 2 i F 項，是由 1 ( ( ) ( )) 2 F UR F UL 的中央差分法加上為了解決不連續面問題的 artificial viscosity term 1
2 A U 所組成。加入 preconditioning 的方程式只需在 artificial viscosity term 做改變即可，其推導如下：
1 1 1 1 1 1 ( ) ( ) ( ) p p p p p p p U F G H S t x y z U F G H S t x y z U U U U A B C S t x y z U U U U AM BM CM S t x y z 其中 p U M U
所以 artificial viscosity terms 改寫如下：
1 1 2 1 1 ( ) 2 R L 2 P i F F F AM U (3-72) 其中 1 1 AM KA DA KA 解完非黏滯性項之後，接著要解的是黏滯性項；在黏滯性項方面，採二階中 央差分法。由於在尤拉方程式中計算的範圍皆為網格與網格之間的通量項，因此 在黏滯性項方面，所需要得到的速度梯度項也必須是網格之間的通量項。下列以 三維的 X 方向為例，圖 3-2 為其示意圖。 圖 3-2 中各編號所代表的位置分別為： 1(i, j1,k)；2 , , ) 2 1 (i j k ；3(i1, j1,k)； 4(i, j,k1)；5 , , 1) 2 1 (i j k ；6(i1,j,k1； 7(i, j,k)；8 , , ) 2 1 (i j k ；9(i1,j,k)； 10(i, j,k1)；11 , , 1) 2 1 (i j k ；12(i1,j,k1)； 13(i, j1,k)；14 , 1, ) 2 1 (i j k ；15(i1,j1,k)； 其各速度梯度差分分別如下表示： X U U X U X U (9) (7) (3-73) X V V X V X V (9) (7) (3-74)
X W W X W X W (9) (7) (3-75) Y U U Y U Y U 2 ) 14 ( ) 2 ( 2 (3-76) 其中 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-77) 同理 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-78) 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-79) Z U U Z U Z U 2 ) 5 ( ) 11 ( 2 (3-80) 其中 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-81) 同理 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-82) 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-83)
