往復運動薄板對管道內加熱壁面熱傳影響之研究

國 立 交 通 大 學

機械工程學系

碩士論文

往復運動薄板對管道內加熱壁面熱傳影響之研究

A Study of a Thin Plate Reciprocated Moving on a Heat

Surface in a Duct

研 究 生：陳又旭

指導教授：傅武雄 博士

往復運動薄板對管道內加熱壁面熱傳影響之研究

A Study of a Thin Plate Reciprocated Moving on a

Heat Surface in a Duct

研 究 生： 陳又旭 Student: Yu-Hsu Chen

指導教授： 傅武雄 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

July 2008

往復運動薄板對管道內加熱壁面熱傳影響之研究

研究生：陳又旭 指導教授：傅武雄 國立交通大學機械工程學系 摘要 本研究利用數值方法分析往復運動的薄板對管道內加熱壁面熱傳的影響；且 近一步討論壁面熱傳的增益。數值計算方面，採用密度基底法解可壓縮流，並且 在方程式中加入 Preconditioning 法，程式中所求解的是網格交界面的通量，方程 式中的差分法皆為顯性差分。首先探討往復運動薄板提升管道內壁面熱傳效率的 機制，進而比較不同薄板移動速度時，對壁面熱傳效率的影響。 由數值計算的結果得知，壁面上往復運動薄板推擠前方及上方流體而破壞其 邊界層。同時牽引後方流體填補往復運動所產生的空間，因而生成新的邊界層， 並且引導流體流向高溫壁面。邊界層的生成與低溫流體流向高溫壁面為提高熱傳 效率的主要機制。此外隨著薄板移動度的增加，邊界層的破壞與生成較為頻繁， 使的熱傳量大幅提升。

II

A Study of a Thin Plate Reciprocated Moving on a Heat

Surface in a Duct

Student : Yu-Hsu Chen Advisor : Wu-Shung Fu

Department of Mechanical Engineering National Chao Tung University

Abstract

The study investigated the heat transfer enhancement by a reciprocated moving plate on the heat surface inside a duct. The characteristics of the flow and thermal fields are analyzed numerically. In the numerical analysis, taking density based method in to consideration, also adding preconditioning method in equations. All finite difference methods are explicit method. Observing the mechanisms of the heat transfer rate then discussing how the different plate moving speed effect the heat transfer rate.

The numerical result shows that the motion of the plate destroys the boundary layer on the heat surface before the moving plate, and the boundary layer reforms behind the moving plate immediately. The boundary layer reformation and cold fluid induced by the moving plate flowing toward the heat surface are the major heat transfer enhancement mechanisms. Besides when the plate moves faster, the boundary layer reforms frequently. As a result, the heat transfer enhancement becomes huge.

誌謝 我由衷的感謝指導老師傅武雄教授在這兩年來給予課業和論文上的指導，以及在 生活各方面上的關心與照顧，在此謹致最高的謝忱與敬意。同時也感謝機械系諸 師長在課業方面的指導。另外要特別感謝博士班學長李崇綱以及連信宏在數值模 擬計算以及軟體使用上的協助以及指導。還有實驗室同窗王俊傑、楊宗霖、許友 維、林敬倫、王威翔、學弟學長們在精神上的鼓勵讓我可以順利突破難關。 更要感謝父母含辛茹苦的養育之恩，與家人和親友的支持與關懷，今日方能順利 完學業。最後僅將此喜悅與所有關心我的人分享。

IV 目錄 中文摘要···Ⅰ 英文摘要···Ⅱ 誌謝···Ⅲ 目錄···Ⅳ 表目錄···Ⅴ 圖目錄···Ⅵ 符號表···Ⅷ 一、緒論···1 二、物理模式···5 2-1 物理尺寸與模型···5 2-2 分析假設及統御方程式···6 2-3 邊界條件及初始狀態···7 三、數值模式···13 3-1 統御方程式 ···14 3-2 Roe scheme···15 3-3 MUSCL 法···22 3-4 Preconditioning··· 24 3-5 Runge-Kutta method··· 30 3-6 座標轉換···32 四、結果與討論···35 五、結論···68 參考文獻···69

表目錄

表 3-1 精度係數值···23 表 4-1 網格與薄板速度···42 表 4-2 Nu number 增加百分比···45

VI 圖目錄 圖 1-1 各種熱傳模式與冷卻流體之熱傳量與可達成溫度差關係圖···4 圖 2-1 管道及薄板之物理模式圖···10 圖 2-2 YZ 平面物理模式圖···11 圖 2-3 XY 平面物理模式圖···12 圖 3-1 黎曼問題特徵值結構圖··· 21 圖 3-2 黏性項差分示意圖··· 29 圖 4.1 XY 平面網格加密圖···43 圖 4.2 YZ 平面網格加密圖···43 圖 4.3 出口速度確認圖··· ·44 圖 4.4.1(a~f) XZ 平面速度向量圖、5~6 週期、雷諾數 225、薄板速度：u* =0.5···46 圖 4.5.1(a~f) YZ 平面速度向量圖、5~6 週期、雷諾數 225、薄板速度：u* =0.5···48 圖 4.6.1(a~f) YZ 平面等溫線圖、5~6 週期、薄板速度：0.005m/s、T_h =310K···50 圖 4.7.1(a) 紐瑟數在網格(11,2,15)、一個週期、薄板速度：u* =0.5···53 圖 4.7.1(b) 紐瑟數在網格(15,2,20)、一個週期、薄板速度：u* =0.5···54 圖 4.7.1(c) 紐瑟數在網格(30,2,25)、一個週期、薄板速度：u* =0.5···55 圖 4.7.1(d) 平均紐瑟數分佈、0~6 週期、薄板速度：u* =0.5···56 圖 4.7.1(e) 總平均紐瑟數增加率、薄板速度：u* =0.5···57 圖 4.7.2(a) 紐瑟數在網格(11,2,15)、一個週期、薄板速度：u* =0.25···58 圖 4.7.2(b) 紐瑟數在網格(15,2,20)、一個週期、薄板速度：u* =0.25···59 圖 4.7.2(c) 紐瑟數在網格(30,2,25)、一個週期、薄板速度：u* =0.25···60 圖 4.7.2(d) 平均紐瑟數分佈、0~4 週期、薄板速度：u* =0.25···61 圖 4.7.2(e) 總平均紐瑟數增加率、薄板速度：u* =0.25···62 圖 4.7.3(a) 紐瑟數在網格(11,2,15)、一個週期、薄板速度：u* =0.75···63 圖 4.7.3(b) 紐瑟數在網格(15,2,20)、一個週期、薄板速度：u* =0.75···64 圖 4.7.3(c) 紐瑟數在網格(30,2,25)、一個週期、薄板速度：u* =0.75···65

圖 4.7.3(d) 平均紐瑟數分佈、0~4 週期、薄板速度：u* =0.75···66 圖 4.7.3(e) 總平均紐瑟數增加率、薄板速度：u* =0.75···67

VIII 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

ρ 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[m s⋅ −1]

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

w velocity component in z -direction[m s⋅ −1]

β pressure gradient[N m⋅ −3] Ac area[m ] 2 t Δ time difference[s] a sound speed[m s⋅ −1] R gas constant[J kg⋅ −1⋅k−1] H enthalpy[ 1 1 J kg⋅ − ⋅k− ]

w

τ friction force acting per unit area on the surface

u_τ friction velocity Greek symbols

v kinematics viscosity[m s⋅ −1]

μ absolute viscosity[kg.m−1.s−1] δ the channel half-width subscripts

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

1

第一章 緒論

在工業界中，高溫壁面熱傳的問題十分常見，邁入 3G 時代，高功率電子 元件的發展與半導體製程技術的進步，促使消費性電子產品的使用與生活型態做 結合，因而目前許多電子產品皆走向高性能、微小化的趨勢，致使電子元件舉凡 電腦 CPU，VGA 或通訊元件 PA 等其封裝元件之發熱密度愈來愈高，其單位熱 通量亦相對地不斷增加，對性能及可靠度等方面均造成不容忽視之影響。因此在 未來電子裝備的發展趨勢走向輕薄短小之際，其性能及可靠度的提升將取決於其 冷卻技術。依據美國航太相關統計資料顯示：電子零件發生故障時，55％主要來 自於溫度因素。因此為了維持元件於額定溫度下運作，必須將此密集的熱量能有 效散逸於系統外之環境。傳統的散熱片大多數是由鋁合金所製造，其熱傳導性只 屬於中等程度，對於目前元件發熱功率越來越高的情況，逐漸有捉襟見肘的情形 發生，目前電子元件的發熱量達到每平方公分數十瓦的等級，且接點可承受的溫 度約在攝氏 150 度以下。相較於圖 1-1 之各種冷卻高溫壁面方法對壁面熱傳量與 溫度差示意圖，此等級的熱傳量已達到空氣強制對流的極限，因此如何提升空氣 強制對流的熱傳效能成為研究的主要課題之一。 過去有多篇文獻提出各種不同的方法來增加高溫壁面的熱傳效能。Hwang 與 Liou[1]以實驗的方法探討在壁面加裝有縫隙的肋狀紊流產生器，對方形管道 流的熱傳與熱阻的影響。結果顯示縫隙狀肋狀紊流產生器可避免壁面產生熱點， 同時對熱傳效能有所增益。Park 等[2]研究不同形狀之熱交換管內插入結構對提 高熱傳效率的影響。結果顯示圓錐形盤管與肋狀的插入結構約可提高 30％的熱 傳效率。Iida 等[3]在液體中加入鋁金屬微粒，當液體發生沸騰現象時，鋁金屬微 粒也被液化，可將熱傳效率提高十倍。Fusegi[4]數值分析有肋狀突起管道中之振 動對向流，當提高振動對黏性的相對強度時，可顯著的提高熱傳效率。Sitter 等 [5]實驗探討高強度音場對重力場與微重力場下池沸騰的影響，發現音場有助於 提高熱傳效率。

比較上述各種增加熱傳量的方法，不論是被動式或是主動式方法增加熱傳效 率似乎都會受到限制。主要的原因是流場在熱傳面上生成的速度邊界層與溫度邊 界層阻礙了熱量傳出。當流體流經高溫壁面時，壁面上的邊界層有速度邊界層與 溫度邊界層，所謂速度邊界層，是指當連續體之黏性流體以一接近速度流經固體 物體時，在雷諾數較高的情況下，流體在固體表面附近受到黏性的影響，形成速 度變化非常激烈且厚度極薄的區域，在此區域內，流體速度從壁面的零漸次變化 到此區域的邊緣速度。而所謂溫度邊界層則是流體受壁面高溫影響，在壁面產生 熱傳，其靠近壁面的溫度逐漸變化至與壁面溫度相同的區域。在溫度邊界層中， 流體的溫度變化和緩，溫度梯度較小。根據熱傳導之傅利葉定律，在溫度梯度較 小的情況下，熱傳量也較小。因此，溫度邊界層的存在將限制壁面高溫所能傳出 之熱傳量。而為了大幅提高高溫壁面的熱傳效率，就必須移除壁面上的溫度邊界 層，使的高溫壁面與低溫流體直接接觸，進而增加高溫壁面與低溫流體的溫度梯 度，達成增加壁面熱傳量的目的。王[6]在其研究中探討了往復運動的薄塊對噴 流中壁面熱傳影響之研究，由王的研究可發現往復運動的薄塊對於壁面邊界層不 斷的破壞與生成，達成提升熱傳效率的目的；但美中不足的是王的物理模式為二 維的流場假設，且其工作流體為不可壓縮流，對於實際工業的應用範圍有限；因 此本研究針對三維的流場作模擬，觀察沿縱軸方向的熱流場變化，再者，本研究 的工作流體為可壓縮流，因此密度隨著壓力與溫度的變化而改變，對於實際工業 界中的應用較為廣泛與實際。 本研究為在一具有一高溫壁面的三維管道流內，在其高溫壁面上加裝一薄平 板，此薄平板直接與高溫壁面接觸，並且在高溫壁面上往復運動。受到薄平板移 動的影響，薄平板前進方向前的流體受到推擠而往薄板移動方向移動，使的薄板 前方的溫度與速度邊界層被破壞。而在薄板後方的流體，則受到薄板的牽引而向 薄板流動，在填補因薄版移動所產生的空間的同時，低溫的流體將與高溫的壁面 直接接觸，因而產生最大的溫度梯度，並且同時在壁面上生成新的溫度與速度邊

3 的溫度與速度邊界層不斷的被移除與重新生成，因此可大幅的提升熱傳效率。而 在薄版移動的同時，流體受到薄版的牽引，使溫度與速度場隨薄版移動而產生變 化，流體與移動物體相互影響的問題，屬於動態的移動邊界問題，又因為在處理 實際散熱問題時，工作流體多為可壓縮流，因此在本研究中的工作流體不論高速 或低速皆視為可壓縮流，以增加應用範圍。此種方法最大的困難處在於計算低速 流場時，由於可壓縮流必須遵守 Courant-Friedrichs-Lewy (CFL)條件，因此在低 速流體時，受限於流體變化傳遞速度(約等於聲速)，時階將會極小。在此種情況 造成計算過程將耗費極大的時間與整體過低的計算效率。為了改善此缺點，本研 究在計算可壓縮流時加入 Preconditioning 法，藉此讓流體即使在低速時，也可有 較高的效率與良好的收斂性。在計算此種低馬赫數流體的方面，目前有密度基底 法(Density-based method)與壓力基底法(Pressure-based method)本研究以密度基底 法為主。而在密度基底法中又以 Turkel[7]提出的 Preconditioning 法最為廣泛應 用，不僅可同時應用在可壓縮流與不可壓縮流中，更可以讓程式的收斂性增加。 而本研究在做數值計算時，主要是利用到網格之間的物理量，因此用到 MUSCL 法來處理網格間的物理量，在計算流場時空間部分以 ROE 法[8]來計算出非黏滯 性項的通量，黏滯性項部分則採用二階中央插分；在時階部份為了能夠使程式能 夠作平型化運算，因此採用三階 Runge-Kutta 插分方法；而當在程式中加入 Preconditioning 法時，同時破壞了統御方程式，為了彌補此一缺點，必須加入 Artificial time term 來修正方程式；最後利用 Dual time stepping[9]來計算暫態的物 理量。在靠近壁面的地方，為了增加程式的收斂效果以及更清楚的觀察壁面效 應，對於壁面網格做加密的處理；而在處理移動邊界的問題時，本研究採用 Immersed boundary[10]中的移動邊界方法，此方法在處理複雜外型及移動邊界的 問題十分的有效率。

5

第二章 物理模式

2.1 物理尺寸與分析模式： 圖 2-1 為一在底面部份加熱管道中加裝一往復移動薄板之物理模式圖。如圖 2.2 以及圖 2.3 所示，一長度l₁為，寬度為w₁，高度為w₂的方形管道，其中l₁ =5m、 m w₁ =0.35 、w₂ =0.35m。溫度T₀ =298K 的冷卻空氣以等速度u=0.01m/s由管 道入口噴入，在管道的底部則部份區域為高溫壁面T ，_h T_h =310K對管道內的空 氣加溫，進行熱交換；在高溫壁面上有一長度為l₂，寬度為w ，高度為₃ h的絕熱 薄板，其中l₂ =0.725m、w₃ =0.0174m、h=0.175m。當時間為零時，薄板靜止 在加熱底面的1 4位置(Ps)，也就是薄板做往復運動的起點，待流場穩定後，薄板 開始以小於進口速度的移動速度往加熱底面的3 4位置(Pe)，也就是薄板移動的終 點移動。當薄板移動又回到最初的起始位置時即完成一週期的運動。本研究即當 薄板在管道內做往復運動時，對管道內的熱場及流場進行分析。x方向為 streamwise 的方向、 y 為垂直壁面方向最後 z 為 spanwise 的方向。u、v、w分 別為其所對應的速度。而為了有效的增進熱傳效率以及要熱管道內的熱流場更容 易成為完全發展流，在管道的長度以及寬度上需要一定的比例，進口速度與薄板 移動速度也有一定的關係，本研究選定管道長度為寬度的 15 倍長，以確保其出 口流體為完全發展流形式；而薄板的移動速度則為進口速度的1 4、 1 2、 3 4，要注 意薄板的移動速度不可大於入口流體速度，這是為了要使入口進來之流體能順利 往出口流動，並且達到完全發展流的條件。

2.2 分析假設及統御方程式： 本研究選擇層流流場作為模擬流場，流場作以下假設： 1. 可壓縮流，空氣密度會隨溫度與壓力而改變。 2. 工作流體為空氣，假設為理想氣體。流體性質為牛頓流體(Newtonian fluid)， 黏滯係數為等方向性。 3. 忽略重力效應影響。 4. 考慮溫度變化對流體所造成的影響 5. 所有壁面均為不可滑移(No slip condition)

統御方程式分別為連續方程式(Continuity equation)、動量方程式(Momentum equations)與能量方程式(Energy equations)，壓力方面則假設流體為理想氣體，利 用理想氣體狀態方程式定義。 連續方程式： ( _j) 0 j u t x ρ _ρ ∂ ₊ ∂ ₌ ∂ ∂ (2-1) 動量方程式： ( ) i j i ij j i j u p u u t x x x ρ _{ρ σ} ∂ ∂ ∂ ∂ + = − + ∂ ∂ ∂ ∂ 其中 2 [ ( )] 3 j i k ij jf ij j i k u u u x x x σ =μ_∈ ∂ +∂ − δ ∂ ∂ ∂ ∂ (2-2) 能量方程式： ( _{j j}) _j ( _{i ij}) j j j E u E pu q u t x x x ρ _{ρ σ} ∂ ₊ ∂ _{+ = −} ∂ ₊ ∂ ∂ ∂ ∂ ∂ (2-3) 理想氣體狀態方程式： p=ρRT (2-4)

7

2.3 邊界條件：

本研究所採用的統御方程式為可壓縮 Navier-stokes 方程式，因此需要給定的 邊界條件有：初始狀態(Initial condition)、入口條件(Inlet condition)、出口條件 (Outlet condition)、壁面邊界(Boundary condition)；另外由於本研究屬於移動網格 問題，因此在方程式中加入沉浸邊界法來解決移動網格的問題。 2.3.1 初始狀態：初始速度、初始壓力、初始密度 初始速度 u： 0.01 /m s 初始速度 v： 0 /m s 初始速度 w： 0 /m s 初始壓力 p： 一大氣壓力(101300Pa) 初始密度ρ： 空氣密度( 3 1.28kg m ) / 2.3.2 入口條件： 入口速度 u： 0.01 /m s 入口速度 v： 0 /m s 入口速度 w： 0 /m s 2.3.3 出口條件： 出口速度 u： 完全發展流 出口溫度 T： 完全發展溫度分布 出口壓力 p： 一大氣壓力(101300Pa) 2.3.4 壁面邊界： 邊界速度： 不可滑移條件，u=v=w=0m/s 邊界溫度： 加熱壁面溫度T_h =310K 邊界密度： 在垂直壁面方向，梯度為零， 0 n ρ ∂ ₌ ∂ 邊界壓力： Neumann 條件（利用 Poisson 方程式來求得邊界壓力）。 根據高斯散度定理，可將 Poisson 方程式寫成如下形式

Poisson 方程式： 2 2 A C p p dy dl y n ∂ ₌ ∂ ∂ ∂

∫∫

∫

(2-5) 所以 2 1 1 2 2 2 y y y p p p p p n y y + − + − ∂ ₌ ∂ ₌ ∂ ∂ Δ 因為 p 0 n ∂ ₌ ∂ ，由此可得到 py−1= py 所以 2 1 1 2 2 2 y y y p p p p y y + − + − ∂ ₌ ∂ Δ 可以改寫成 1 2 y y p p y + − Δ 因此在邊界的壓力可以設定成p_y+1= p_y 2.3.5 沉浸邊界法(Immersed boundary) 處理複雜形狀的能力是數值模擬的動大問題之ㄧ。許多實際之移動或變形邊 界之非線性流體系統問題模擬，不可避免的需要於每個時間項網格重建，因而造 成計算困難度。過去於複雜計算領域中，貼體法(Body-fitting)即沉浸邊界法 (Immersed boundary method)成功的模擬外形複雜流場。Peskin[10]發展沉浸邊界 法以作為廣義座標變換的替換，於控制方程式中增加外力項以代表流場內任何形 狀的物體，以避免轉換過程。增加流體內特定節點之外力項可模擬物體於流場之 影響，故可在卡式座標中模擬任意形狀物體。沉浸邊界法的開發，除了適當的處 理流場外型複雜度，同時於網格節點保留計算結果之準確性。因此在卡式座標中 可處理移動或變形之複雜外型流場問題，卻不需要隨時間產生網格。沉浸邊界法 的特色在於運動方程式中加入源項。為了產生邊界效應所提出的外立方程式通常 有兩種方式。分別為連續外力逼近與離散外力逼近。第一種方法主要是在離散為 分方成之前加入外力方程式。連續外力法的應用將流經沉浸彈性邊界的物理機制 耦合於方程式中。離散外力法則是在方程式被離散後才引入外力項。因此相較於 連續外力逼近，離散外力逼近較不直接，而加入外力方式與離散方式有直接關 係，可使沉浸邊界有較準確表現。根據沉浸邊界上不同描述邊界條件的方法，離 散外力法可分類為間接跟直接使用邊界條件。而本研究採用的是離散外力法，藉

9 界鄰近點之動量方程式可表示為： ( ) i j i ij j i j u p u u f t x x x ρ _{ρ σ} ∂ ∂ ∂ ∂ + = − + + ∂ ∂ ∂ ∂ (2-6) 式中外力項 f 可由下式之離散動量方程式獲得 1 ( ) ( ) n n j i ij j i j u u p u u f t x x x ρ + _{ρ σ} ∂ − ₊ ∂ _{= −}∂ ₊ ∂ ₊ Δ ∂ ∂ ∂ (2-7) 此處外加力 f 僅作用於物體邊界或物體內，經由每一次時間步伐計算獲得以滿足 沿任意沉浸邊界Ω 上之速度使其等於其給定之速度函數V_Ω。施加之外力項可表 示為 n V u f RHS t Ω− = − + Δ (2-8) 上式中Δt為時間增量，V_Ω為沿邊界上給定之速度，RHS 為動量方程式中之壓力 梯度、對流及擴散向組成。如果對於靜止物體，給定之無滑動邊界條件V_Ω = 。 0 經移項之後可以得到 i i n i n i f RHS t u u + = Δ − + ) ( 1 ρ (2-9) 再將外力項 f 帶入方程式中可以得到最後結果 1 1 + + ₌ n i n i V u (2-10)

11

13

第三章 數值計算模式

本章主旨在說明本論文的數值計算所使用的所有模式。 第一節首先要知道所求 解的方程式為 Navier-Stokes 方程式。本研究將方程式拆解為非黏滯項與黏滯項。 第二節介紹的為黎曼解中的 ROE 法，利用 ROE 法來求出非黏滯項的通量。 接 著第三節介紹 MUSCL 法，此法是為了要解出 ROE 法中使用的網格之間的物理 量，然後為了防止在高階插分時產生震盪現象，在 MUSCL 法插分的結果方程式 中加入 Minmod limiter 以確保程式不會發散。 第四節為介紹 Preconditioning 法， 因為當計算低速可壓縮流時，因速度和音速的數量及上差距過大，在數值分析時 不好計算，所以為彌補此一缺點須使用 Preconditioning 法。 最後第五節為三階 的 Runge-Kutta Scheme，在時間項的插分時為了配合程式的平行化運算，因而使 用 Runge-Kutta 法作為計算的分法，而程式因為在使用 Preconditioning 時，加入 Artificial time term 時，已破壞了整個統御方程式，因此需使用 Dual time stepping 疊代使其在 Artificial domain 收斂時才能進入下一個真實時階，將在此小節做詳 細的解說。 第六節為座標的轉換，本研究為了要觀察邊界的地方流場及熱場的 變化及邊際效應的影響，在靠近壁面的網格部份做了加密的動作，因此當程式作 計算時使用曲線座標系做計算，但是當插分時是在直角座標系做插分的動作，因 此對於方程式座標的轉換在此節做說明。綜合上述，本論文在數值上的計算過程 為，利用 MUSCL 法算出 ROE 法所需要的網格間物理量，求出非黏滯的通量， 並且在計算通量時加入 Preconditioning 法，以拉近與音速的 Order。 接下來使用 2 階中央插分法對黏滯項做插分進而求出黏滯項；然後再與 ROE 法求出的通量 結合得到真正的物理通量。最後使用 Runge-Kutta 法疊代以求出下一時階的物理 量。

3-1、統御方程式： 本研究在計算流場的方面其統御方程式分可為兩大部分，第一部份為非黏滯性項 的尤拉方程式， 第二部份為黏滯性項。 3 1 2 0 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 mj j m u u u P F u u P u u P T e P u u x ρ ρ δ τ ρ δ τ ρ δ τ ρ λ τ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − ⎟ ⎜ ⎟ =_⎜ + − _⎟ ⎜ _{+ −} ⎟ ⎜ ⎟ ∂ ⎜ _{+ − −} ⎟ ⎜ _∂ ⎟ ⎝ ⎠ ∀ =m 1, 2,3。 (3-2) ρ為密度， P 為壓力。u 、₁ u 、₂ u 分別為₃ x、 y 、 z 方向的速度。λ ρ= C k_p ，k 為 thermal diffusivity。 2 2 2 1 2 3 1 ( ) 2 v e=C T+ u +u +u ，C 為等容比熱。_v τ_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組成的方程式即稱為尤拉方程式。

15 3-2、Roe scheme： 在此節首先介紹一維線性黎曼方程式，並導出其 exact solution。其次將其推 廣至非線性的黎曼方程式，同時也利用 Roe 法求出其近似解。 首先考慮一維線性黎曼方程式如下： 0 U U A t x ∂ ₊ ∂ ₌ ∂ ∂ ，其中 A 為一常數 Jacobian 矩陣。 (3-4) 初始條件為 (0) (0) (0) 1 ( , , _m )T U = u _L u 。左上方括號代表時間 t 為 0。 求出 A 之特徵值矩陣以及特徵向量。 1 A= ΛK K− ，其中Λ 為特徵值矩陣： 1

0

0

m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ L M O M L 。 (1) ( ) , , m T K = ⎣⎡K _L K ⎤_⎦ 為特徵向量，故 ( )i ( )i i AK =λK 。 接著定義特徵變數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 ∂ ₊ ∂ ₌ ∂ ∂ ，或 1 1 1 2 2 3 3 ... 0 0 ... 0 0 : : : : : 0 ... m t x w w w w w w λ λ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ₊⎢ ⎥ ⎢ ⎥ ₌ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (3-6) (3-6)可由特徵曲線法求得其解為： (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 i i p U x t αK βK = = + =

∑

+

∑

(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[15]所 提出，亦即為 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 就稱為 joacobian 矩陣。

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

Jacobian matrix〉A U U%( _L, _R)代替，因此本來的黎曼問題可以改寫成近似黎曼問題： ( ) 0 U U A U t x ∂ ∂ + = ∂ % ∂ ( , 0) U x =U _L x<0 (3-12)

17

於是前述方法可以得到(3-10)的近似解。由以上的原理可得知，在近似黎曼 問題上，Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性 轉變成線性，但是初始條件並沒有改變，因此可以得到方程式(3-10)的近似解。 為了要求得合理的常數 Jacobian 矩陣，須合乎 Roe 所提出的四項條件： 1. U與 F 之間，存在著線性轉換的關係。 2. 當U_R −U_L→ ，則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) 接著在固定體積的條件下，積分近似解 ₁ 2 (0) i U + ，可得到通量(flux)的數值公式：

1 1 2 2 ( (0)) ( _R) ( _R) i i F F U F U F U + = % + − − % (3-16) 再從F% = %AU的關係中可進一步求得： 1 1 2 2 (0) ( _R) _R i i F AU F U AU + = % + − − % (3-17) 再根據(3-13)式與(3-14)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i R i i i _i F F U A K F U K λ α λ α+ + = − %

∑

_> % = −

∑

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

∑

_> % = +

∑

₌ % % % (3-19) (3-18)與(3-19)所指的λ%_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-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

0

m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ L % _{M O M} L 。 接下來需找出 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為聲速 方程式(3-22)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示： 2 2 0 1 ( ) 2 F A U a u u U ⎡ ⎤ ∂ = _{= ⎢ ⎥} − ∂ _{⎣ ⎦} (3-23)

19 特徵向量： (1) 1 K u a ⎡ ⎤ = ⎢ ₋ ⎥ ⎣ ⎦， (2) 1 K u a ⎡ ⎤ = ⎢ ₊ ⎥ ⎣ ⎦ 接著選定 parameter vector Q 1 2 q U Q q _u ρ ρ ρ ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥}= ⎢ ⎥ ⎣ ⎦ ⎢_⎣ ⎥_⎦ (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 1 f q q F f q a q ⎡ ⎤ ⎡ ⎤ =_{⎢ ⎥ ⎢}= _⎥ + ⎣ ⎦ ⎣ ⎦ (3-26) 為了表示出Δ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-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)，可以求得 1 2 1 2q 0 B q q ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ % % % % ； 2 1 2 2 1 2 2 q q C a q q ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ % % % % % (3-31) 再帶入(3-30)可得 2 2 0 1 2 A a u u ⎡ ⎤ = ⎢ ₋ ⎥ ⎣ ⎦ % % % (3-32)

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 則是利用 MUSCL 法求出。 _R

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

3-3、Monotonic Upstream-Centered Scheme for Conservation Laws(MUSCL)：

本論文使用的是採用I. Abalakin etl.[14]中所使用的插分法。其方程式如下：

1/ 2 1/ 2 1/ 2 L L i i i u₊ = +u Δu₊ (3-38) 1/ 2 1/ 2 1/ 2 R R i i i u₊ = −u Δu₊ (3-39) 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-40) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i i u₊ β u₊ u β u₊ u₊ Δ = − − + − +θc(−u_i−1+3u_i−3u_i+1+u_i+2) 1 2 3 ( 3 3 ) d i i i i u u u u θ _{+ + +} + − + − + (3-41) 其中(3-40)、(3-41)式中的β、θc 、θd 值可由表(3-1)中查得。代入不同的值可以 得到不同的精度。本論文則是使用三階精度，以減少數值計算的消散性。 在程式中，高次項的插分法在不連續的情況下，容易使震盪變大，為了降低震盪， 本研究在 MUSCL 法插分出來的方程式中加入 minmod limiter，用來確保程式不 會發散。 因此(3-38)與(3-39)式需改寫如下： 1/ 2 1/ 2 min mod( 1/ 2) L L i i i u₊ = +u Δu₊ (3-42) 1/ 2 1/ 2 min mod( 1/ 2) R R i i i u₊ = −u Δu₊ (3-43)

23 表 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 法：

為了提高程式的應用範圍，加入 preconditioning 法，讓程式不論在高速或低速流 體內，皆可獲得精確的計算結果。本程式採用 Weiss and Smith 的 preconditioning method [17]，方程式改變如下： 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 T p p T p T p u u U v v M U w w H u v w H C ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ∂ _{⎢ ⎥} = = ∂ ⎢ ⎥ ⎢ ⎥ ⎢ _{− +} ⎥ ⎣ ⎦ (3-46) 其中 _p p ρ ρ =∂ ∂ ； T T ρ ρ = ∂ ∂ 接著將(3-45)式的方程式乘上矩陣 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-47) 再將 K 與 M 相乘 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 p T KM C ρ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢₋ ⎥ (3-48)

25 將(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 條件的限制， 提高程式的計算效率。 利用θ 取代ρ_p項： 2 1 1 ( ) r p U TC θ = − (3-52) max U ε× if u < ×ε C r U = u if ε× <C u <C (3-53) C if u >C 其中ε為一極小的值，約等於 5 10−，其主要是用來防止停滯點（stagnation point） 在計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，U 必須_r

大於流體的當地擴散速度（local diffusion velocity），因此U 還需加入下列限制： _r

max( , ) r r U U x ν = Δ 將θ帶入(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)式轉變如下： ( ) 0 p nc U F G H K t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-55) 為了讓(3-55)式中的通量項再度轉換成保守形式，在乘上 1 K− 1 (K _nc) Up F G H 0 t x y z −_Γ ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (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 ρ θ ρ θ ρ ρ θ ρ ρ θ ρ ρ θ ρ ρ ρ ρ − − ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ₋ ⎥ ⎢ ⎥ Γ = Γ = ⎢ ⎥ ⎢ ₋ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ − ⎥ − + ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (3-57) 最後方程式簡化成如下形式： 0 p U _{F G H} t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-58) 由於方程式在時間項經過改變，因此必須重新推導 Roe 所提出的近似黎曼解。 在(3-21)式中，可以觀察到 ₁ 2 i F + 項，是由 1 ( ( ) ( )) 2 F UR +F UL 的中央差分法加上為 了解決不連續面問題的 artificial viscosity term 1

2 A U% Δ 所組成。加入

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

所以 artificial viscosity terms 改寫如下：

1 1 2 1 1 ( ) 2 R L 2 P i F F F − AM U + = + − Γ Δ (3-59) 其中 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) 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) W W W W(10)+ (12) (4)+ (6)

29

3-5、Runge-Kutta method：

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

首先，先在原始 Navier-Stokes 方程式加入一虛擬時間項，稱為 artificial time term。 方程式改變如下： 0 U U F G H t x y z τ ∂ ₊∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ ∂ (3-71) 其中τ 即為 artificial time， t 為 physical time

接著在 artificial time term 加入 preconditioning 的方法：

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

最後對 artificial time term 採一階的有限差分離散，對 physical time term 採二階的 後項差分離散，可得 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 ( ) ( ) ( ) 0 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 t x y z τ + + − + + + + + + − + − + − + − − + Γ + + − + − + − = Δ Δ Δ Δ Δ (3-73) 為了求得方程式中的 k 1 p U + ，本研究採用顯示三階 Runge-Kutta 法，選用此法主要 是因為此法不僅程式易於撰寫，且為了程式計算的速度，此法易於平行化，加速 程式的計算時間。其方程式如下所示： 1 1 3 2 1 1 1 3 3 3 2 2 1 3 1 3 3 1 1 4 4 4 1 2 2 3 3 3 k k k p p k k k k p p p k k k k p p p U U R U U U R U U U R + ₋ + + ₋ + + + + − = + Γ = + + Γ = + + Γ (3-74) 其中 _{1 1 1 1 1 1} , , , , , , , , , , , , 2 2 2 2 2 2 1 1 1 [ ( ) ( ) ( )] i j k i j k i j k i j k i j k i j k R F F G G G G S x y z τ Φ Φ Φ Φ Φ Φ Φ Φ + − + − + − = −Δ − + − − + Δ Δ Δ

31 上述方程式，當 artificial time term 收斂至

1 0 k k p p U U τ + ₋ Γ = Δ 時，方程式即會回復到 原始的 Navier-Stokes 方程式，並且包含著時間項，故程式可以計算暫態結果。 接著整理(3-73)式，先將其線性化 1 3( ) 4 ( ) ( ) ( ) 0 2 k n n p p k k k k k k x p p y p p z p p U U M U U U F A U G B U H C U t δ δ δ τ − Δ + Δ − + Γ + + + Δ + + Δ + + Δ = Δ Δ (3-75) 其中 k 1 k p p p U U + U Δ = − 、 _{1 1} , , , , 2 2 1 ( ) k k k x i j k i j k F F F x δ + − = − Δ 、Ap =AM 、 1 1 , , , , 2 2 1 ( ) ( ) k k k x P P P i j k i j k A A A x δ + − ⎡ ⎤ = _⎢ − _⎥ Δ ⎣ ⎦。 再將ΔU_p項，放置於等號左邊，其餘則於等號右邊： 1 3 1 [ ( )] 2 k k k k x p y p z p p I M A B C U R t δ δ δ τ − − + Γ + Γ + + Δ = Δ Δ (3-76) 此處 1 3 4 ( ) ( ) 2 k n n k k k k x y z U U U R F G H t δ δ δ − − + = − − + + Δ ， I 為單位矩陣。 1 1 3 1 [ ( )] 2 k k k p p k k k x p y p z p R U U I M A B C t δ δ δ τ + − − = + + Γ + Γ + + Δ Δ (3-77) Eq. (3-77)同樣採用顯示三階 Runge-Kutta 法求解。

3-6、座標轉換：

本研究為了要使程式的應用範圍更廣，在程式中加入了座標轉換，本節即是推導 座標轉換後的方程式變化。 首先，原座標系統，也就是卡式座標系統，我們稱 作 Physical space ( , , )x y z ，轉換過後的曲線座標系統則稱為 Computational

space ( , , )ξ η ς 。其關係式如下： ( , , , ) ( , , , ) ( , , , ) t t x y z t x y z t x y z τ ξ ξ η η ς ς = = = = (3-78) 經由 Chain rule 對偏微分作計算，可以知道兩座標間的轉換公式如下： t t t x x x y y y z z z t x y z ξ η ς τ ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ∂ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∂ ∂ ₌ ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∂ (3-79)

由上式可以知道在 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-80) 本研究為在靠近壁面的地方做網個加密的動作，而此加密函數的轉換矩陣，轉換 公式如下： ( ) 1 ( ) 1 1 (2 )( ) 2 1 1 (2 1)(1 ( ) ) 1 Y H η α α η α α β α β α β β β α β − − − − + + + − − = + + + − (3-81) (2 1) 2 [ ] (2 1) 2 (1 ) 1 y H Log y H α β α α β α η α α _β + + − + − + = + − + (3-82)

33 當程式在做計算時，需要的是偏微分η_y及 y_η的值，其偏微分的結果如下： ( ) 1 ( ) ( ) ₂ 1 1 1 1 (2 )( ) [ ] 1 1 1 1 ( 1 )(2 1)(( ) ( ) ) 1 1 H Log Y η α α η α η α α β β β β β β β α α β β + − − − + + − − = − + + − + + + − − (3-83) ( 1 ) ( ) ( ) ₂ 1 1 1 1 1 1 ( 1 )(2 1)( ) (( ) ( ) ) 1 1 1 2 Y HLog α η α η α α α β β β α α β β β η β − − + − − + − + + + + − + + + − − − = − (3-84) 接著考慮統御方程式： 0 U F G H t x y z ∂ ₊∂ ₊∂ ₊∂ ₌ ∂ ∂ ∂ ∂ (3-85) 本研究將黏滯項與非黏滯項拆開，方程式變為： F G H U F G H t x y z x y z ν ν ν ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + + = + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ (3-86) LHS = t t t x x x y y y z z z U U U U F F F G G G H H H ξ η ς τ ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ξ η ς ∂ ∂ ∂ ∂ + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + + ∂ ∂ ∂ ∂ ∂ ∂ + + + ∂ ∂ ∂ ∂ ∂ ∂ + + + ∂ ∂ ∂ (3-87) 同理對 RHS 也做一樣的計算，最後對 LHS 以及 RHS 同除以轉換矩陣 J，統御方 程式最後結果如下： LHS = 1 ( ) [ ( )] 1 [ ( )] 1 [ ( )] t x y z t x y z t x y z U U F G H J J U F G H J U F G H J ξ ξ ξ ξ τ ξ η η η η η ς ς ς ς ς ∂ ∂ + + + + ∂ ∂ ∂ + + + + ∂ ∂ + + + + ∂ (3-88)

RHS = 1 [ ( )] 1 [ ( )] 1 [ ( )] x y z x y z x y z F G H J F G H J F G H J ν ν ν ν ν ν ν ν ν ξ ξ ξ ξ η η η η ς ς ς ς ∂ + + ∂ ∂ + + + ∂ ∂ + + + ∂ (3-88) 上式即為本研究在做計算時所使用的統御方程式，在計算時一樣使用 Roe scheme 來解非黏滯性項的答案，並且加入 Preconditioning 法以及最後使用 Runge-Kutta 法來做 Dualtime steping 的暫態計。

35

第四章 結果與討論

本研究探討往復運動的薄板對高溫壁面附近溫度與速度場的影響，並說明高 溫壁面上的平均紐瑟數隨著薄板運動而產生的變化，進而比較不同的薄板運動速 度對流場與溫度場的影響與對熱傳效率的提升。本研究總共計算三個不同薄板速 度的案例，如表格 4.1 所示。接著要說明的是，本研究在程式計算時為有因次的 計算，但是輸出的數據在做圖時轉為無因次的表示法，主要是為了方便討論與使 用無因次之參數。以下為無因次化結果： * u u U_∞ = 、 * v v U_∞ = 、 * w w U_∞ = 、U_∞ =0.01 * * 2 * 2 * 2 ( ) ( ) ( ) Vmag = u + v + w c h c T T T T T − − = * ，Tc =298K，Th =310K * 2 1 2 p p U ρ_{∞ ∞} = 再次確認我們的物理模式圖如圖 2.2 所示，可以知道加熱底板為等溫，薄板 為絕熱，薄板延著 z 方向做往復運動如圖 2.3 所示。 接著要驗證的是本研究使用的座標轉換公式是否正確，如第三章所說明，本 研究使用的方程式為： ( ) 1 ( ) ( ) ₂ 1 1 1 1 (2 )( ) [ ] 1 1 1 1 ( 1 )(2 1)(( ) ( ) ) 1 1 H Log Y η α α η α η α α β β β β β β β α α β β + − − − + + − − = − + + − + + + − − (3-72) ( 1 ) ( ) ( ) ₂ 1 1 1 1 1 1 ( 1 )(2 1)( ) (( ) ( ) ) 1 1 1 2 Y HLog α η α η α α α β β β α α β β β η β − − + − − + − + + + + − + + + − − − = − (3-73) 其中β =1.2為其加密係數、α =0.5是決定其加密以 y 軸中央為對稱。網格加密 的結果如圖 4.1、及圖 4.2 所示，可以很清楚的看出來網格加密只在 y 方向加密， x 和 z 方向的網格為等間距，與本研究的預期相符合。 再來要確認速度場的正確性，本研究初使條件為方形管道完全發展流的出口

速度解析解，其出口速度解析解如下： ( 1) 2 2 3 3 1,3,5,... cosh( ) cos( ) 16 _{2 2} ( , ) ( ) ( 1) [1 ] cosh( ) 2 i i i z i y a dp _{a a} u y z i b dx i a π π π μπ − ∞ = = −

∑

− − × (4-1) 3 5 5 1,3,5,... tanh( ) 4 192 ₂ Q ( )[1 ] 3 i i b ba dp a _a dx b i π μ π ∞ = = − −

∑

(4-2) Q 為流量，其結果如圖 4.3 所示，可以看出來本研究的數值解和解析解十分的吻 合，出口最大速度大約為進口速度的兩倍，但是因為壁面效應消耗了些許能量， 因此比 0.02m/s 略小一點，由此結果可以知道本研究的初始條件正確。 當確認過了網格與初始流場的正確性之後，接著要確認的是流場與熱場模擬 的結果，圖 4.4.1(a~f)至圖 4.6.1(a~f)分別為雷諾數為 225、無因次薄板移動速度 為u* =0.5的情況，隨著薄板往復運動而在高溫壁面上到達不同的位置時，其 xz 平面速度向量圖、yz 平面速度向量圖、等溫線圖與高溫壁面上的特定點的紐瑟 數分布圖。由於本研究探討三維的空間，因此除了 x 方向之外還探討了 z 方向的 速度變化，以看清楚管道內因薄板移動所造成之流場。 第一個 case 是當薄板移動速度為 u* =0.5時的結果。 圖 4.4.1(a)至圖 4.4.1(c)、與相對應的 4.5.1(a)~4.5.1(c)為薄板從網格 z 位置 15 移動至位置 20 再移 動至位置 25，薄板由左向右移動，往復運動前半週期的階段，其相對於在 z 軸 上的物理位置分別為 0.1312m、0.175m、0.2188m。如圖 4.4.1(a)所示，薄板在網 格 15 的位置時，薄板已由起始位置向右移動了五個網格的距離，由於圖 4.4.1(a) 為 xz 平面上視圖的結果，所以薄板在圖中移動方向為上，在圖 4.5.1(a)中可以看 出來薄板右方的流體受到擠壓而向右流動。較接近薄板位置的流體受到薄板的推 擠較為劇烈，其流動速度較遠處的流體為快，因此接進薄板前方的流體受擠壓而 向上方流動。而由於薄板向右移動，在薄板左方的流體為了填補薄板移動所產生 的空間而亦向右流動，同時在薄板後方生成新的速度邊界層。此外，薄板向右運

37 小的迴流產生。 當薄板繼續移動至網格 20 與網格 25 的位置時，如圖 4.5.1(b)， 薄板前方的流體持續被推擠，薄板後方的流體也不斷的向右填補，使的薄板上方 與前方的迴流加劇，而隨著薄板向右移動到達右端的終點後，薄版右端的流體衝 突現象逐漸消失。 圖 4.4.1(d)至圖 4.4.1(f)與 4.5.1(d)至 4.5.1(f)則為薄板到達右端 終點後向左移 5 格、10 格、與 15 格的距離，為往復運動後半週期的階段。由圖 4.5.1(d)中可以看出來，當薄板往左移動時，同樣的薄板左方的流體被推擠而堆 積，最後往薄板上方及前方流動，而薄板後方的流體不斷的補進來，因而產生迴 流區，同樣的，在薄板走過之後產生新的邊界層。從圖 4.4.1(a)至圖 4.4.1(f)為止， 為第六次週期往復運動的結果；而因為薄板的移動速度為u* =0.5，為入口速度 的 1/2，對於程式的收斂性有相當的影響，因為薄板在第一次週期之後對管道內 的流體流場已經造成了震盪，因此當薄板進行第二次到第六次週期的往復運動 時，可以知道在 x 方向以及 z 方向的擾動會較第一次週期為大，但也因此使的熱 傳的效果增加。 第二個 case 為薄板移動速度為u* =0.25時的結果，由於薄板移動速度較慢， 第一次週期後的結果和薄板移動第二次到第四次週期之後的結果十分的相似，由 此可以知道管道內部流場在程式計算時都有達到收斂，誤差很小，原因是此時的 薄板移動速度為u* =0.25為入口速度的 1/4，因為薄板移動較慢，對管道內流場 的影響較小，因此流場較容易達到收歛。 第三個 case 為薄板移動速度為u* =0.75時的結果，因為薄板的移動速度過 快，所以在薄板移動了半個週期之後，其內部流場的擾動十分的大，隨著薄板移 動從一個週期至四個週期，可以知道在薄板的前方延 x 方向以及 z 方向有三維擾 動現象產生，並且隨著週期越來越大，擾動延伸的範圍越廣泛，因此可以知道薄 板移動的速度有其極限值，而其極限值主要和其入口速度相關，當薄板移動速度 越接近入口速度時，其管道內部的流體受到薄板推擠的能量越大，因此內部流場 的變化與擾動越大，當其能量大過入口速度給流體的能量時，從入口進來之流體 將無法留至出口，程式因此會無法收歛，所以薄板移動速度的極限值就是不能大

於入口速度。 接著要討論熱場的變化，圖 4.6.1(a~f)為u* =0.5且與圖 4.4.1(a~f)以及圖 4.5.1(a~f)相同位置及條件下薄板與高溫壁面間的等溫線圖。由於薄板在高溫壁面 上往復週期運動，薄板移動方向前方的邊界層被移動中的薄板破壞，且薄板移動 方向前方的流體受推擠而往薄板上方及前方向流動，因此如圖 4.6.1(a~c)所示， 薄板前進方向前方高溫壁面上的等溫線分布較疏，熱傳效率也較差。相反的，在 薄板的移動方向的後方，熱邊界層重新生成在高溫壁面上，同時受迴流區與薄板 移動的牽引，流體朝向高溫壁面流動，使的薄板移動方向後方的高溫壁面上等溫 線分布較為密集，熱傳量亦較大。而隨著薄板的移動速度改變，溫度場的變化也 隨之而產生變化。當薄板的移動速度為u* =0.25時，由於流場速度相對於薄板移 動速度為大，此時管道內流場較為穩定，同樣的管道內熱場變化也較為穩定，因 此熱傳量較差。而當薄板移動速度為u* =0.75s 時由於管道內流場變動較為劇 烈，連帶的影響到管道內的熱場變化較為密集且變化劇烈，因此有較大的熱傳 量。由於u* =0.25和u* =0.75的等溫線分布圖和u* =0.5的情況十分相似，所以 不以圖示說明。

圖 4.7.1(a~c)、 4.7.1(a~c)以及圖 4.7.3(a~c)分別為u* =0.5、u* =0.25和 75 . 0 * = u 時的局部紐瑟數分布圖，其圖 abc 的物理位置分別為位置(11,2,15)、 (25,2,20)和(40,2,25)在相同薄板位置情況下高溫壁面上局部紐瑟數分佈圖，局部 紐瑟數的定義如下式： 0 = ∂ ∂ − = = Y x x Y k W h Nu θ 式中h 為局部熱對流係數；W 為管道的代表長度；k 為冷卻流體的熱傳導係數。 _x 首先討論薄板移動速度為u* =0.5時的情況；圖 4.7.1(a)中的實線(-)代表薄板 由左往右進行前半週期的運動時，在位置(11,2,15)的紐瑟數變化，可以看出來， 受到絕熱薄板的影響，紐瑟數的分布被切為兩部分，前半段紐瑟數下降，是由於