以高階通量限制函數之壓力修正法應用無結構性網格求解全速流流場

國 立 交 通 大 學

機 械 工 程 學 系

博士論文

以高階通量限制函數之壓力修正法

應用無結構性網格求解全速流流場

Pressure-Based Unstructured-Grid Algorithms

Incorporating High-Resolution Schemes for

All-speed Flow Calculations

研 究 生：吳 添 成

指導教授：崔 燕 勇

以高階通量限制函數之壓力修正法應用無結構性網格

求解全速流流場

Pressure-Based Unstructured-Grid Algorithms Incorporating

High-Resolution Schemes for All-speed Flow Calculations

研 究 生：吳添成 Student：Tian-Cherng Wu 指導教授：崔燕勇 Advisor：Yeng-Yung Tsui 國 立 交 通 大 學 機 械 工 程 學 系 博 士 論 文 A Thesis

Submitted to Institute of Mechanical Engineering College of Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Doctor of Philosophy

in

Inctitute of Mechanical Engineering January 2009

Hsinchu, Taiwan, Republic of China

以高階通量限制函數之壓力修正法應用無結構性網格

求解全速流流場

研究生：吳添成 指導教授：崔燕勇 博士 國立交通大學機械工程學系 摘 要 本文以壓力基底法發展可執行低速不可壓縮流到高速超音速可壓縮流 之全速流流場計算方法，運用有限體積法、重置變數及任意邊形之無結構 性網格來離散統御方程式，為了處理震波附近陡峭的梯度變化，採用全變 量消去法(TVD)或正規化變數(NV)法導出之高階對流通量限制函數，通量限 制子則由偵測兩個連續的梯度比來自動調整。在梯度的計算採用二階線性 重置限制修正法，以強化計算過程的穩定性及解的準確度。 本文分別使用原始變數及守恆變數作為求解變數，以原始變數求解 時，壓力修正方程式係由密度變量及速度變量各別與壓力變量之關係式代 入連續方程式所導出，並藉由流場局部馬赫數來自動調整橢圓型式或雙曲 線型式之壓力方程式。以守恆變數求解時，壓力修正方程式只有守恆速度 (ρV )變量與壓力變量之關係式，為了模擬超音速區域之雙曲線型式之流場 特性及穩定計算過程，則使用上風型式之遲滯密度或遲滯壓力作修正。 採用一些策略來增加求解過程的穩定性，諸如(1)對流項使用一階隱式 上風差分混合顯式高階差分之遲緩修正法；(2)將擴散項分解為只包含相鄰 網格之隱式正交項，以及非正交部份之顯式垂直導數修正項並將其置於源 項中；(3)採用不同的局部時步方式，所有控容體之時步由固定的 Courant 數來決定，即網格較小則時步較小；(4)以鬆弛法來進行線性方程的疊代。 上述的方法會導致對角係數值的擴大，因此會使得係數矩陣具有對角佔優。 本文發展之計算解子允許網格為任意多邊形，計算所用之網格可以使 用不同來源的網格產生方式，在本研究中亦發展一套整合式網格介面處理

程式，可以將不同方式產生之不同邊數之區塊網格予以結合，並轉換產出 滿足吾人發展之計算解子所需計算域網格資料。 經由數種流場測試來驗證本文發展的方法，黏性流有(1)流經圓柱之低 速流、(2)低速空穴流、(3)流經 NACA 0012 翼型外流場、(4)雙喉部噴嘴內 流場等。非黏性流則包括(1)漸縮-漸擴噴嘴內流場、(2)流經下壁面圓弧之渠 道流、(3)流經 NACA 0012 翼型外流場、(4)流經圓柱之高速流場、(5)流經 三角柱之高速流場等。由測試結果證明，不論是原始變數或守恆變數求解 方式所建構之流場解子，均能執行低速不可壓縮流到高速可壓縮流之層流 流場計算，均能獲得準確的收歛解且能準確地捕捉高速流場中震波的位置 及強度。

Pressure-Based Unstructured-Grid Algorithms Incorporating

High-Resolution Schemes for All-speed Flow Calculations

Student：Tian-Cherng Wu Adviser：Dr. Yeng-Yung Tsui Institute of Mechanical Engineering

National Chiao Tung University

ABSTRACT

Pressure-based algorithms applicable to all-speed flows, ranging from incompressible to supersonic flows, are developed in this thesis. The finite volume method is employed for discretization. The grids, which can be of arbitrary topology, are arranged in collocated manner. To tackle the abrupt change of gradient in the region of shock, either the total variation diminishing (TVD) scheme or the normalized variable (NV) scheme can be incoporated via the use of flux limiting function. These flux limiters are determined from the ratio of two consecutive gradients. To enhance solution accuracy, the gradients are calculated using a second-order linear reconstruction approach.

In this study, the mathematical formulation is based on either the primitive variables or the conservative variables. In the model using the primitive variables, a pressure-correction equation is obtained from the continuity equation by using the relations between the variations of the velocities and density and that of the pressure. The resulted equation is of mixed type, either elliptic or hyperbolic, depending on the local Mach number. The second model consider the variation of the pressure with the conserved velocities (ρV ). To

account for the hyperbolic character of the supersonic flows, either the density or the pressure is retarded in the upwind direction.

Several strategies are adopted to enhance the stability of the solution iteration procedure as follows: (1) The convective flux is composed of a upwind

part and an anti-diffusion part. The upwind part is treated implicitly and the other part explicitly; (2) The diffusive flux is divided into a part in the direction directed from the considering node to the neighboring node and a part normal to this direction. The former is tackled in an implicit manner while the latter is absorbed into the source term; (3) The time step for each control volume is based on the cell Courant number. With a fixed Courant number for all control volumes, the time steps are different for the control volumes. The smaller the cell volume, the smaller the time step; (4) The difference equations are under-relaxed during iteration. The above methods can enlarged the diagonal coefficients and ,thus, make the coeffient matrix more diagonal dominant.

The algorithm developed allows the control volumes of the meshes to be a polygon of arbitrary geometry. Different sources of grid generator can be adopted to generate computational meshes. An interface is developed to combine the meshes generated in different blocks using different grid generation methods and transfer the grid data into the format required by our computational code.

The methodology is validated via testing on a number of flows. For viscous flows there are (1) low-speed flows over a cylinder, (2) low-speed flows in a cavity, (3) flows over a NACA 0012 airfoil and (4) flows in a double throats. In inviscid flow, test cases include (1) flows in a convergent-divergent nozzle, (2) flows in a channel with a circular arc bump, (3) flows over a NACA 0012 airfoil, (4) high-speed flows over a cylinder, (5) high-speed flows over a triangle. Accurate results can be obtained effectively using the developed methods, regardless of the use of primitive or conservative variables, for the flows ranging from the incompressible to high-speed compressible flows. It is seen that the location and the strength of the shock waves in high-speed flow can be accurately predicted.

誌

謝

衷心感謝指導老師崔燕勇教授，在我博士學涯7 年多來在論文研究上 殷切耐心的指導與鼓勵，每每在山窮水盡疑無路時，總能適時的給予正確 指引而出現柳暗花明的新轉機，讓我能走過面臨重重難題而百Run 不得其 解的恢心低潮期，在此再致上深深謝意。 感謝劉振隆博士提供NACA 翼型網格產生程式及趙茂吉博士提供三邊 形無結構性網格產生程式，以及在CFD 領域上的討論，對我的論文研究助 益斐淺。亦要感謝中山科學研究院一所(航空研究所)准予以公餘進修方式攻 讀博士學位。 特別感謝父母親多年的養育之恩及永遠默默地支持我的愛妻筑萱，在 我博士班的學涯中都由她獨攔家務及照料兩位聰明、活潑可愛的女兒斐蓁 及冠誼，如今她們已是高二及小五了，沒有她的付出，我就無法專心地致 力於研究工作，就沒有今日的我，僅以此小小的成就獻給她們。 最後感謝每一位曾在人生路上伴我成長的師長、同事、朋友、學弟伙 伴們，謝謝您們！ 吳添成 2009.01.14 于台中

目 錄

摘 要 ... i 目 錄 ... vi 表目錄 ... x 圖目錄 ... xi 符號說明...xxviii 第 1 章 緒 論... 1 1.1. 簡介... 1 1.2. 文獻回顧... 6 1.3. 研究方向... 17 1.4. 研究貢獻... 19 1.4.1. 全速流流場解子程式發展... 19 1.4.2. 整合式網格組合程式發展... 19 1.4.3. 論文發表... 20 1.5. 論文綱要... 20 第 2 章 統御方程式... 22 2.1. 積分型式統御方程式... 22 2.2. 散度型式統御方程式... 23 2.3. 張量型式統御方程式... 24 第 3 章 數值方法... 25 3.1. 傳輸方程式之離散... 26 3.1.1. 時間項... 26 3.1.2. 對流項... 27 3.1.3. 擴散項... 29 3.1.4. 源項... 30 3.1.5. 線性代數方程... 33 3.1.5.1. 動量方程式... 34 3.1.5.2. 能量方程式... 34

3.1.6. 鬆弛處理... 35 3.2. 壓力、密度及速度偶合關係式... 36 3.2.1. 面上質量流率的處理... 36 3.2.2. 面上速度的處理... 37 3.2.3. 面上密度的處理... 39 3.2.4. 壓力修正線性代數方程... 40 3.3. 邊界條件... 45 3.4. 疊代殘值計算及收歛條件... 48 3.5. 梯度計算模式... 49 3.6. 求解程序... 53 第 4 章 高階通量限制函數法... 54 4.1. 簡介... 54 4.2. 高階準確通量限制函數之對流項離散... 57 4.3. 以特徵變數求對流通量之高階限制反擴散項... 61 4.3.1. 雙曲線型系統方程式... 62 4.3.2. 特徵變數通量限制子之計算... 63 4.4. 求解程序... 66 第 5 章 守恆變數求解法... 68 5.1. 遲滯密度法... 68 5.2. 遲滯壓力法... 71 5.3. 傳輸方程式的離散... 73 5.3.1. 對流項... 73 5.3.2. 擴散項... 74 5.3.3. 線性代數方程... 75 5.4. 壓力、密度及速度偶合關係式... 76 5.4.1. 面上質量流率的處理... 76 5.4.2. 壓力修正線性代數方程... 77 5.4.3. 邊界條件... 79 5.4.4. 求解程序... 79

第 6 章 網格產生... 81 第 7 章 結果與討論... 83 7.1. 以原始變數法求解... 84 7.1.1. 非黏性流之驗證測試... 84 7.1.1.1. 漸縮-漸擴噴嘴內流場 ... 84 7.1.1.2. 流經下壁面圓弧之渠道內流場... 89 7.1.1.3. NACA 0012 翼型外流場... 105 7.1.2. 黏性流之驗證測試... 107 7.1.2.1. 低速空穴流流場... 107 7.1.2.2. 流經圓柱之低速流外流場... 108 7.1.2.3. 雙喉部噴嘴內流場... 110 7.1.2.4. NACA 0012 翼型外流場... 114 7.2. 以特徵變數通量限制函數法求解... 114 7.2.1. 非黏性流場之驗證測試... 115 7.2.1.1. 漸縮-漸擴噴嘴內流場 ... 115 7.2.1.2. 流經下壁面圓弧之渠道內流場... 116 7.2.1.3. NACA 0012 翼型外流場... 117 7.2.1.4. 流經圓柱之高速流外流場... 118 7.2.1.5. 斜震波流場解析－流經三角柱之高速流場... 118 7.2.2. 黏性流之驗證測試... 120 7.2.2.1. 低速空穴流流場... 120 7.2.2.2. 流經圓柱之低速流外流場... 121 7.2.2.3. 雙喉部噴嘴內流場... 121 7.2.2.4. NACA 0012 翼型外流場... 122 7.3. 以守恆變數法求解... 123 7.3.1. 遲滯壓力法... 123 7.3.1.1. 漸縮-漸擴噴嘴內流場 ... 123 7.3.1.2. 流經下壁面圓弧之渠道內流場... 124 7.3.1.3. NACA 0012 翼型外流場... 125 7.3.2. 以遲滯密度法求解... 126

7.3.2.1. 漸縮-漸擴噴嘴內流場 ... 126 7.3.2.2. 流經下壁面圓弧之渠道內流場... 127 7.3.2.3. NACA 0012 翼型外流場 ... 128 第 8 章 結論... 130 第 9 章 參考文獻... 133 附 表 ... 150 附 圖 ... 156 附錄A：特徵變數限制法公式推導... 310 附錄B：整合式網格轉換產生程式(IGTP23D)發展及使用說明 ... 315 附錄C：二維翼型計算網格產生方法說明 ... 331 附錄D：正規化變數圖(NVD)... 333 附錄E：Sweby 的 TVD 關係圖 ... 338 附錄F：無結構性網格之正規化變數計算法 ... 342 簡 歷： ... 345

表 目 錄

頁次 表4-1：線性 NVD 與 TVD 限制函數關係式 ··· 150 表4-2：非線性 NVD 與 TVD 限制函數關係式 ··· 151 表7-1：圓柱低雷諾數外流場計算數值比較表 ··· 154 表7-2：雙喉部噴嘴內流場計算數值比較表 ··· 155

圖 目 錄

圖1-1：結構性網格示意圖··· 156 圖1-2：四邊形無結構性網格示意圖··· 156 圖1-3：三邊形無結構性網格示意圖··· 157 圖3-1：對流項離散化之標示圖··· 157 圖3-2：擴散項離散化 over-relaxed 法之標示圖··· 158 圖3-3：壁面剪應力計算之示意圖··· 158 圖3-4：對稱面垂直應力計算之示意圖··· 158 圖3-5：雙喉部噴嘴出口流場(a)邊界條件以壁面剪應力方式處理(出口靠壁 面處產生異常迥流)；(b)出口正交網格(出口仍產生小迥流)；(c)邊 界條件直接解動量方程式(無迥流)··· 159 圖3-6：雙喉部噴嘴出口正交網格··· 159 圖3-7：邊界壓力計算之示意圖··· 160 圖3-8：入口或出口邊界之參考點示意圖 ··· 160 圖3-9：入口給定停滯邊界條件處理之示意圖 ··· 161 圖3-10：出口給定靜壓力邊界條件處理之示意圖 ··· 161 圖3-11：典型三角形無結構網格梯度計算之積分路徑標示圖 ··· 162 圖4-1：對流項離散化之標示圖··· 162 圖4-2：(a)流場原始變數示意圖；(b)流場正規化變數示意圖 ··· 163 圖4-3：正規化變數圖(NVD)(a)對流界限準則(CBC)示意圖；(b)TVD 示意圖 ··· 163 圖4-4：Sweby 的 TVD 示意圖··· 163 圖4-5：SMART 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 164 圖4-6：STOIC 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 164 圖4-7：UMIST 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 164 圖4-8：WACEB 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 165 圖4-9：MUSCL 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 165

圖4-10：GAMMA 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 165 圖4-11：SUPERBEE 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 166 圖4-12：MINMOD 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 166 圖4-13：OSHER 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 166 圖4-14：Koren 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 167 圖4-15：CUBISTA 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 167 圖4-16：H-QUICK 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 167 圖4-17：CHARM 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 168 圖 4-18：Van-Leer 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 168 圖4-19：OSPRE 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 168 圖4-20：Hemker 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 169 圖4-21：Van-Albada 高階限制函數(a)NVD 圖；(b)TVD 關係圖 ··· 169 圖4-22：CLAM 高階限制函數(a)NVD 圖；(b)TVD 關係圖··· 169 圖5-1：遲滯壓力或密度處理之示意圖··· 170 圖6-1：自由流流經圓柱的計算域網格··· 170 圖6-2：自由流流經三角柱的計算域網格 ··· 171 圖6-3：NACA 0012 翼型計算網格··· 171 圖7-1：漸縮-漸擴噴嘴之幾何外形示意圖··· 172 圖7-2：二維漸縮-漸擴噴嘴無結構性網格圖(a)200×20CV；(b)100×10 CV172 圖7-3：一維 C-D 噴嘴之馬赫數分佈圖(CDS/UDS 混合法)··· 172 圖7-4：二維 C-D 噴嘴之馬赫數分佈圖(CDS/UDS 混合法)··· 173 圖7-5：二維 C-D 噴嘴之馬赫數分佈圖(SMART) ··· 173 圖7-6：二維 C-D 噴嘴之馬赫數分佈圖(STOIC) ··· 174 圖7-7：二維 C-D 噴嘴之馬赫數分佈圖(UMIST) ··· 174 圖7-8：二維 C-D 噴嘴之馬赫數分佈圖(WACEB)··· 175 圖7-9：二維 C-D 噴嘴之馬赫數分佈圖(MUSCL)··· 175 圖7-10：二維 C-D 噴嘴之馬赫數分佈圖(GAMMA,bm=0.25)··· 176 圖7-11：二維 C-D 噴嘴之馬赫數分佈圖(GAMMA,bm=0.5)··· 176 圖7-12：二維 C-D 噴嘴之馬赫數分佈圖(GAMMA,bm=0.75)··· 177

圖7-13：二維 C-D 噴嘴之馬赫數分佈圖(SUPERBEE)··· 177 圖7-14：二維 C-D 噴嘴之馬赫數分佈圖(MINMOD) ··· 178 圖7-15：二維 C-D 噴嘴之馬赫數分佈圖(OSHER) ··· 178 圖7-16：二維 C-D 噴嘴之馬赫數分佈圖(Koren)··· 179 圖7-17：二維 C-D 噴嘴之馬赫數分佈圖(CUBISTA)··· 179 圖7-18：二維 C-D 噴嘴之馬赫數分佈圖(H-QUICK)··· 180 圖7-19：二維 C-D 噴嘴之馬赫數分佈圖(CHARM)··· 180 圖7-20：二維 C-D 噴嘴之馬赫數分佈圖(Van Leer) ··· 181 圖7-21：二維 C-D 噴嘴之馬赫數分佈圖(OSPRE) ··· 181 圖7-22：二維 C-D 噴嘴之馬赫數分佈圖(Hemker)··· 182 圖7-23：二維 C-D 噴嘴之馬赫數分佈圖(Van Albada) ··· 182 圖7-24：二維 C-D 噴嘴之馬赫數分佈圖(CLAM) ··· 183 圖7-25：二維下壁面圓弧渠道(a)幾何外形示意圖；(b)次/穿音速流場計算網 格90×30 CV；(C)超音速流場計算網格 240×80 CV ··· 184 圖7-26：Min=0.5 次音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (UDS/CDS 混合法) ··· 185 圖7-27：Min=0.5 次音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SUPERBEE) ··· 185 圖7-28：Min=0.5 次音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Albada)··· 186 圖7-29：Min=0.675 穿音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈 圖(CDS/UDS 混合法)··· 186 圖7-30：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(SMART) ··· 187 圖7-31：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(STOIC) ··· 187 圖7-32：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(UMIST) ··· 187 圖7-33：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(WACEB)··· 187 圖7-34：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(MUSCL) ··· 188

圖7-35：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(GAMMA) (a) bm=0.25；(b) bm=0.5··· 188

圖7-36：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(SUPERBEE) (a)dm=2；(b)dm=0.5 ··· 189 圖7-37：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(MINMOD)··· 189 圖7-38：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(OSHER)··· 189 圖7-39：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(Koren)··· 189 圖7-40：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(CUBISTA) ··· 190 圖7-41：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(H-QUICK) ··· 190 圖7-42：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(CHARM) ··· 190

圖7-43：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(Van Leer) ··· 190

圖7-44：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(OSPRE) ··· 191

圖7-45：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(Hemker & Koren)191 圖7-46：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(Van Albada) ··· 191

圖7-47：Min=0.675 穿音速渠道流場壁面馬赫數分佈圖(CLAM) ··· 191

圖7-48：(a)接觸(Attached)及脫離(Detached)斜震波之示意圖；(b)正規反射 波(Regular reflection)之示意圖；(c)馬赫反射波(Mach reflection)之 示意圖··· 192 圖7-49：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (UDS/CDS 混合法) ··· 193 圖7-50：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SMART)··· 193 圖7-51：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (STOIC)··· 194 圖7-52：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (UMIST)··· 194 圖7-53：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (WACEB) ··· 195 圖7-54：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MUSCL) ··· 195 圖7-55：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖

(GAMMA，bm=0.5) ··· 196 圖7-56：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SUPERBEE) ··· 196 圖7-57：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MINMOD) ··· 197 圖7-58：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSHER) ··· 197 圖7-59：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Koren) ··· 198 圖7-60：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CUBISTA) ··· 198 圖7-61：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (H-QUICK) ··· 199 圖 7-62：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CHARM)··· 199 圖 7-63：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Leer)··· 200 圖 7-64：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSPRE)··· 200 圖7-65：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖

(Hemker & Koren) ··· 201

圖7-66：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Albada)··· 201 圖7-67：Min=1.65 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CLAM)··· 202 圖7-68：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CDS/UDS 混合法) ··· 202 圖7-69：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SMART)··· 203

圖7-70：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (STOIC)··· 203 圖7-71：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (UMIST)··· 204 圖7-72：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (WACEB) ··· 204 圖7-73：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MUSCL) ··· 205

圖7-74：Min=1.4 馬赫數等值圖及壁面馬赫數分佈圖(GAMMA) (a)

bm=0.25；(b) bm=0.5··· 206 圖7-75：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SUPERBEE，dm=0.4) ··· 206 圖7-76：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MINMOD) ··· 207 圖7-77：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSHER) ··· 207 圖7-78：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Koren) ··· 208 圖7-79：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CUBISTA) ··· 208 圖7-80：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (H-QUICK) ··· 209 圖7-81：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CHARM)··· 209 圖7-82：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Leer)··· 210 圖7-83：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSPRE)··· 210 圖7-84：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖

(Hemker & Koren) ··· 211 圖7-85：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Albada)··· 211 圖7-86：Min=1.4 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CLAM)··· 212 圖7-87：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CDS/UDS 混合法) ··· 212 圖7-88：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SMART)··· 213 圖7-89：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (STOIC)··· 213 圖7-90：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (UMIST)··· 214 圖7-91：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (WACEB) ··· 214 圖7-92：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MUSCL) ··· 215 圖7-93：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (GAMMA) ··· 215 圖7-94：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (SUPERBEE) ··· 216 圖7-95：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (MINMOD) ··· 216 圖7-96：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSHER) ··· 217 圖7-97：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Koren) ··· 217 圖7-98：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CUBISTA) ··· 218

圖7-99：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (H-QUICK) ··· 218 圖7-100：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CHARM)··· 219 圖7-101：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Leer)··· 219 圖7-102：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (OSPRE)··· 220 圖7-103：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Hemker & Koren) ··· 220

圖7-104：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (Van Albada)··· 221 圖7-105：Min=1.6 超音速渠道流場(a)馬赫數等值圖；(b)壁面馬赫數分佈圖 (CLAM)··· 221 圖7-106：渠道內流場三邊形網格(a)厚度比 t/c=10% (15234 CV)；(b)厚度比 t/c=4% (24562 CV) ··· 222 圖7-107：Min=0.5 次音速流場(a)壁面馬赫數分佈比較圖；(b)u 速度通量限 制子分佈圖··· 223 圖7-108：Min=0.675 穿音速流場(a)壁面馬赫數分佈比較圖；(b)馬赫數等值 圖；(c)u 速度通量限制子分佈圖··· 224 圖7-109：Min=1.4 超音速流場(a)壁面馬赫數分佈比較圖；(b)馬赫數等值圖； (c)u 速度通量限制子分佈圖··· 225 圖7-110：Min=1.65 超音速流場(a)壁面馬赫數分佈比較圖；(b)馬赫數等值 圖；(c)u 速度通量限制子分佈圖··· 226 圖7-111：NACA 0012 翼型示意圖 ··· 226 圖7-112：NACA 0012 翼型 C-型四邊形無結構性網格 ··· 226 圖7-113：NACA 0012 翼型 O-型四邊形無結構性網格 ··· 227 圖7-114：NACA 0012 翼型 O-型四邊形與三邊形之混合網格 ··· 227 圖7-115：NACA 0012 次音速流場 M∞=0.63，α=2°，dm=1.0(Van Albada)(a)

馬赫數等值圖；(b)翼表面壓力分佈··· 227 圖7-116：NACA 0012 穿音速流場 M∞=0.80，α=1.25°，dm=0.6 (Van Albada)

(a)馬赫數等值圖；(b)翼表面壓力分佈 ··· 228 圖7-117：NACA 0012 穿音速流場 M∞=0.85，α=1.0°，dm=0.6 (Van Albada)

(a)馬赫數等值圖；(b)翼表面壓力分佈 ··· 228 圖7-118：NACA 0012 超音速流場 M∞=1.20，α=0.0°，dm=0.6 (Van Albada)

(a)馬赫數等值圖；(b)翼表面壓力分佈 ··· 228 圖7-119：穴流(a)四邊形計算網格(QUSG); (b)三邊形計算網格(TUSG)··· 229 圖7-120：穴流，Re=1000，原始變數法 SUPERBEE 限制函數(a)流線圖；(b) 速度向量圖；(c)x, y 中心截面速度分佈；(d)殘值收歛情形(QUSG) ··· 229 圖7-121：穴流，Re=1000，原始變數法 Van Albada 通量限制函數(a)流線圖；

(b)速度向量圖；(c)x, y 中心截面速度分佈；(d)殘值收歛情形(QUSG) ··· 230

圖7-122：穴流，Re=1000，原始變數法 SUPERBEE 通量限制函數(a)流線圖；

(b)速度向量圖；(c)x, y 中心截面速度分佈；(d)殘值收歛情形 (TUSG)··· 231 圖7-123：穴流，Re=1000，原始變數法 Van Albada 通量限制函數(a)流線圖；

(b)速度向量圖；(c)x, y 中心截面速度分佈；(d)殘值收歛情形 (TUSG)··· 232

圖7-124：圓柱外流場計算網格(a)混合網格(HUSG)；(b)四邊形網格(QUSG)

··· 233 圖7-125：圓柱邊界流場分離角定義示意圖 ··· 233 圖7-126：圓柱低速流 Re=20 (Van Albada 限制函數，QUSG 網格) (a)流線圖；

(b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收歛情 形；(d)沿著圓柱表面的渦度及壓力分佈圖··· 234 圖7-127：圓柱低速流 Re=40 (Van Albada 限制函數，QUSG 網格) (a)流線圖；

(b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收歛情 形；(d)沿著圓柱表面的渦度及壓力分佈圖··· 235

圖7-128：圓柱低速流 Re=20 (Van Albada 限制函數，HUSG 網格) (a)流線圖； (b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收歛情 形；(d)沿著圓柱表面的渦度及壓力分佈圖··· 236 圖7-129：圓柱低速流 Re=40 (Van Albada 限制函數，HUSG 網格) (a)流線圖；

(b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收歛情 形；(d)沿著圓柱表面的渦度及壓力分佈圖··· 237 圖7-130：圓柱低速流 Re=20 (SUPERBEE 限制函數，HUSG 網格) (a)流線

圖；(b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收 歛情形；(d)沿著圓柱表面的渦度及壓力分佈圖 ··· 238 圖7-131：圓柱低速流 Re=40 (SUPERBEE 限制函數，HUSG 網格) (a)流線

圖；(b)沿著在圓柱下游的中心線(y=0)上之 u-速度分佈；(c)殘值收 歛情形；(d)沿著圓柱表面的渦度及壓力分佈圖 ··· 239 圖7-132：雙喉部噴嘴幾何外形··· 239 圖7-133：雙喉部噴嘴計算域網格(315×40 CV)··· 240 圖7-134：雙喉部噴嘴馬赫數等值圖(a)Reo= 100；(b) Reo= 400；(c) Reo= 1600 ··· 240 圖7-135：雙喉部噴嘴 Reo=100 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖··· 241 圖7-136：雙喉部噴嘴 Reo=400 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖··· 242 圖7-137：雙喉部噴嘴 Reo=1600 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖··· 243 圖7-138：雙喉部噴嘴，原始變數 SUPERBEE 限制函數法(a)壁面壓力分佈 圖；(b)壁面摩擦係數分佈圖；(c)對稱中心馬赫數分佈圖；(d)對稱 中心壓力分佈圖··· 244

圖7-139：雙喉部噴嘴，原始變數 Van Albada 限制函數法(a)壁面壓力分佈

圖；(b)壁面摩擦係數分佈圖；(c)對稱中心馬赫數分佈圖；(d)對稱 中心壓力分佈圖··· 245

(b) Van Albada··· 246 圖7-141：NACA 0012, M∞=0.8, α=10°, Re∞=500，原始變數 SUPERBEE 限

制函數(dm=0.2)(a)馬赫數等值圖；(b)流線圖；(c)收歛殘值；(d)翼 面壓力分佈比較圖··· 247 圖7-142：NACA 0012, M∞=0.50, α=0°, Re∞=5000，原始變數 Van Albada 限 制函數(a)馬赫數等值圖；(b)尾部流線圖；(c)收歛殘值；(d)翼面壓 力分佈比較圖··· 247 圖7-143：C-D 噴嘴內流場，特徵變數 MD1, Van Albada 限制函數法，噴嘴 中心軸馬赫數分佈圖··· 248 圖7-144：C-D 噴嘴內流場，特徵變數 MD1, SUPERBEE 限制函數法，噴嘴 中心軸馬赫數分佈圖··· 248 圖7-145：C-D 噴嘴內流場，特徵變數 MD2, Van Albada 限制函數法，噴嘴 中心軸馬赫數分佈圖··· 249 圖7-146：C-D 噴嘴內流場，特徵變數 MD2, SUPERBEE 限制函數法，噴嘴 中心軸馬赫數分佈圖··· 249 圖7-147：C-D 噴嘴內流場，無梯度修正之 SUPERBEE 與 Van Albada 限制

函數法結果比較(a)噴嘴中心馬赫數分佈；(b)震波局部放大圖·· 250 圖7-148：C-D 噴嘴內流場，梯度修正 1 次之 SUPERBEE 與 Van Albada 限

制函數法結果比較(a)噴嘴中心馬赫數分佈；(b)震波局部放大圖251 圖7-149：C-D 噴嘴內流場，梯度修正 2 次之 SUPERBEE 與 Van Albada 限

制函數法結果比較(a)噴嘴中心馬赫數分佈；(b)震波局部放大圖252 圖7-150：C-D 噴嘴內流場，特徵變數 MD1 之計算殘值比較(a)Pb=0.870Po； (b)Pb=0.769Po；(c)Pb=0.645Po··· 255 圖7-151：C-D 噴嘴內流場，特徵變數 MD2，計算殘值比較(a)Pb=0.870Po； (b)Pb=0.769Po；(c)Pb=0.645Po··· 258 圖7-152：Min=0.5 次音速渠道流場壁面馬赫數分佈比較圖 ··· 259 圖7-153：Min=0.675 穿音速渠道流場壁面馬赫數分佈比較圖 ··· 259 圖7-154：Min=1.4 超音速渠道流場壁面馬赫數分佈比較圖 ··· 260 圖7-155：Min=1.65 超音速渠道流場壁面馬赫數分佈比較圖 ··· 260

圖7-156：NACA 0012 次音速流場 M∞=0.63，α=2° (a)馬赫數等值圖；(b) 翼表面壓力分佈(特徵變數法，SUPERBEE 限制函數) ··· 261 圖7-157：NACA 0012 穿音速流場 M∞=0.8，α=1.25° (a)馬赫數等值圖；(b) 翼表面壓力分佈(特徵變數法，SUPERBEE 限制函數) ··· 261 圖7-158：NACA 0012 穿音速流場 M∞=0.85，α=1.0° (a)馬赫數等值圖；(b) 翼表面壓力分佈(特徵變數法，SUPERBEE 限制函數) ··· 261 圖7-159：NACA 0012 超音速流場 M∞=1.2，α=0.0° (a)馬赫數等值圖；(b) 翼表面壓力分佈(特徵變數法，SUPERBEE 限制函數) ··· 262 圖7-160：NACA 0012 穿音速流場 M∞=0.80，α=1.25° (a)馬赫數等值圖；(b) 翼表面壓力分佈(原始變數限制函數法)··· 262 圖7-161：NACA 0012 穿音速流場 M∞=0.80，α=1.25° (a)馬赫數等值圖；(b) 翼表面壓力分佈(特徵變數限制函數法)··· 262 圖7-162：NACA 0012 穿音速流場 M∞=0.80，α=1.25°，原始變數及特徵變 數通量限制法之翼表面壓力分佈比較圖··· 263 圖7-163：圓柱外流場四邊形(128×8 CV)與三邊形(18556 CV)混合網格·· 263 圖7-164：M∞=6 超音速圓柱外流場，原始變數 UDS 法(a)馬赫數等值圖； (b)殘值收歛情形；(c)圓柱表面之壓力分佈 ··· 264 圖7-165：M∞=6 超音速圓柱外流場，特徵變數 MD1, dm=0.9 (a)馬赫數等值 圖；(b) 殘值收歛情形；(c)圓柱表面之壓力分佈··· 264 圖7-166：M∞=6 超音速圓柱外流場，特徵變數 MD2, dm=0.8 (a)馬赫數等值 圖；(b) 殘值收歛情形；(c)圓柱表面之壓力分佈。··· 265 圖7-167：斜震波示意圖··· 265 圖7-168：斜震波θ-β-M 關係圖(摘自 Wikipedia 網站)··· 266 圖7-169：三角柱外流場計算網格··· 266 圖7-170：M∞=2.6 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬赫 數等值圖；(c)流線圖··· 267 圖7-171：M∞=2.54 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬 赫數等值圖；(c)流線圖··· 267 圖7-172：M∞=2.53 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬

赫數等值圖；(c)流線圖··· 267 圖7-173：M∞=2.52 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬 赫數等值圖；(c)流線圖··· 268 圖7-174：M∞=2.50 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬 赫數等值圖；(c)流線圖··· 268 圖7-175：M∞=2.40 三角柱超音速流場(a)馬赫數等值圖局部放大圖；(b)馬 赫數等值圖；(c)流線圖··· 268 圖7-176：三角柱附近流場之流線及速度向量分佈圖(M∞=2.53) ··· 269 圖7-177：穴流 Re=1000，特徵變數 MD-1, SUPERBEE (a)流線圖；(b)速度

向量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形(QUSG)270 圖7-178：穴流 Re=1000，特徵變數 MD-1, Van Albada 限制函數(a)流線圖；

(b)速度向量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形 (QUSG) ··· 271 圖7-179：穴流 Re=1000，特徵變數 MD-2, SUPERBEE (a)流線圖；(b)速度

向量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形(QUSG)272 圖7-180：穴流 Re=1000，特徵變數 MD-2, Van Albada (a)流線圖；(b)速度向

量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形(QUSG)273 圖7-181：穴流 Re=1000，特徵變數 MD1, SUPERBEE (a)流線圖；(b)速度向

量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形(TUSG)274 圖7-182：穴流 Re=1000，特徵變數 MD1, Van Albada (a)流線圖；(b)速度向

量分佈圖；(c)x,y 中心截面速度分佈；(d)殘值收歛情形(TUSG)275 圖7-183：圓柱低速流 Re=20，特徵變數 MD1, Van Albada 限制函數法(a)流

線圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈 ··· 276 圖7-184：圓柱低速流 Re=40，特徵變數 MD1, Van Albada 限制函數法(a)流 線圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈 ··· 276 圖7-185：圓柱低速流 Re=20，特徵變數 MD2, Van Albada 限制函數法(a)流 線圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈

··· 277 圖7-186：圓柱低速流 Re=40，特徵變數 MD2, Van Albada 限制函數法(a)流

線圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈 ··· 277 圖7-187：圓柱低速流 Re=20，特徵變數 MD1, SUPERBEE 限制函數(a)流線 圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈278 圖7-188：圓柱低速流 Re=40，特徵變數 MD1, SUPERBEE 限制函數(a)流線

圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈278 圖7-189：圓柱低速流 Re=20，特徵變數 MD2, SUPERBEE 限制函數(a)流線

圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈279 圖7-190：圓柱低速流 Re=40，特徵變數 MD2, SUPERBEE 限制函數(a)流線

圖；(b)尾流 u 速度分佈；(c)收歛殘值；(d)壁面渦度與壓力分佈279 圖7-191：雙喉部噴嘴 Reo=100 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖(特徵變數法) ··· 280 圖7-192：雙喉部噴嘴 Reo=400 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖(特徵變數法) ··· 281 圖7-193：雙喉部噴嘴 Reo=1600 (a)馬赫數等值圖；(b)壓力等值圖；(c)溫度 等值圖(特徵變數法) ··· 282 圖7-194：雙喉部噴嘴，特徵變數 MD1, SUPERBEE 限制函數法(a)壁面壓力 分佈圖；(b)壁面摩擦係數分佈圖；(c)對稱中心馬赫數分佈圖；(d) 對稱中心壓力分佈圖··· 283 圖7-195：雙喉部噴嘴，特徵變數 MD2, SUPERBEE 限制函數法(a)壁面壓力 分佈圖；(b)壁面摩擦係數分佈圖；(c)對稱中心馬赫數分佈圖；(d) 對稱中心壓力分佈圖··· 284 圖7-196：雙喉部噴嘴，特徵變數 MD1, Van Albada 限制函數法 (a)壁面壓

力分佈圖；(b)壁面摩擦係數分佈圖；(c)對稱中心馬赫數分佈圖； (d) 對稱中心壓力分佈圖··· 285

圖7-197：雙喉部噴嘴，特徵變數 MD2, Van Albada 限制函數法(a)壁面壓力

對稱中心壓力分佈圖··· 286 圖7-198：雙喉部噴嘴，特徵變數 MD1 數值收歛情形(a) SUPERBEE 限制函 數；(b) Van Albada 限制函數 ··· 287 圖7-199：雙喉部噴嘴，特徵變數 MD2 數值收歛情形(a) SUPERBEE 限制函 數；(b) Van Albada 限制函數 ··· 288 圖7-200：NACA 0012, M∞=0.8, α=10°, Re∞=500，特徵變數 MD1, SUPERBEE 通量限制函數(dm=0.7)(a)馬赫數等值圖；(b)尾部流線 圖；(c)收歛殘值；(d)翼面壓力分佈比較圖 ··· 289 圖7-201：NACA 0012, M∞=0.8, α=10°, Re∞=500，特徵變數 MD2, SUPERBEE 通量限制函數(dm=1)(a)馬赫數等值圖；(b)尾部流線 圖；(c)收歛殘值；(d)翼面壓力分佈比較圖 ··· 289 圖7-202：NACA 0012, M∞=0.50, α=0°, Re∞=5000，特徵變數 MD1, Van

Albada 通量限制函數(a)馬赫數等值圖；(b)尾部流線圖；(c)收歛殘 值；(d)翼面壓力分佈比較圖··· 290 圖7-203：NACA 0012, M∞=0.50, α=0°, Re∞=5000，特徵變數 MD2, SUPERBEE 通量限制函數(a)馬赫數等值圖；(b)尾部流線圖；(c) 收歛殘值；(d)翼面壓力分佈比較圖··· 290 圖7-204：一維 C-D 噴嘴之中心馬赫數分佈圖(遲滯壓力法)··· 291 圖7-205：一維 C-D 噴嘴之中心馬赫數分佈圖(遲滯壓力法)··· 291 圖7-206：一維 C-D 噴嘴之中心馬赫數分佈圖(遲滯壓力法，加密網格) · 292 圖7-207：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯壓力法) · 292 圖7-208：二維 C-D 噴嘴中心馬赫數分佈圖 (Pb=0.769Po，遲滯壓力法)293 圖7-209：二維 C-D 噴嘴中心馬赫數分佈圖 (Pb=0.769Po，遲滯壓力法)293 圖7-210：二維 C-D 噴嘴中心馬赫數分佈圖 (Pb=0.769Po，遲滯壓力法)293 圖7-211：渠道流場計算網格(256×64 CV，t/c=10%) ··· 294 圖7-212：Min=0.5 次音速流場壁面馬赫數分佈圖(遲滯壓力法)··· 294 圖7-213：Min=0.5 次音速流場馬赫數等值圖(遲滯壓力法)··· 294 圖7-214：Min=0.5 次音速流場 Van Albada 及 SUPERBEE 限制函數法殘值

圖7-215：Min=0.675 穿音速流場馬赫數等值圖(遲滯壓力法)··· 295 圖7-216：Min=0.675 穿音速流場隨 Mref 及 k 值變化之壁面馬赫數分佈圖(遲 滯壓力法)··· 296 圖7-217：Min=0.675 穿音速流場隨 Mref 及 k 值變化之壁面馬赫數分佈圖 (Van Albada，遲滯壓力法) ··· 296 圖7-218：Min=0.675 穿音速流場隨 Mref 及 k 值變化之壁面馬赫數分佈圖 (SUPERBEE，遲滯壓力法)··· 297 圖7-219：NACA 0012 次音速流場 M∞=0.63，α=2°，dm=2.0 (SUPERBEE) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯壓力法) ··· 297 圖7-220：NACA 0012 穿音速流場 M∞=0.80，α=1.25°，dm=0.8(SUPERBEE)

(a)馬赫數等值圖；(b)翼表面壓力分佈(遲滯壓力法)··· 298 圖7-221：NACA 0012 穿音速流場 M∞=0.85，α=1.0°，dm=1(Van Albada) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯壓力法) ··· 298 圖7-222：NACA 0012 超音速流場 M∞=1.20，α=0.0°，dm=1(Van Albada) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯壓力法) ··· 298 圖7-223：一維 C-D 噴嘴之中心馬赫數分佈圖(遲滯密度法)··· 299 圖7-224：一維 C-D 噴嘴之中心馬赫數分佈圖(遲滯密度法，加密網格) · 299 圖7-225：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯密度法) · 300 圖7-226：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯密度法) · 300 圖 7-227：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯密度法) · 300 圖 7-228：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯密度法) · 301 圖 7-229：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.870Po，遲滯密度法) · 301 圖 7-230：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.769Po，遲滯密度法) · 301 圖7-231：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.769Po，遲滯密度法) · 302 圖7-232：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.769Po，遲滯密度法) · 302 圖7-233：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.769Po，遲滯密度法) · 302 圖7-234：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.769Po，遲滯密度法) · 303 圖 7-235：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.645Po，遲滯密度法) · 303 圖 7-236：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.645Po，遲滯密度法) · 303

圖7-237：二維 C-D 噴嘴中心馬赫數分佈圖(Pb=0.645Po，遲滯密度法) · 304 圖7-238：二維 C-D 噴嘴中心馬赫數分佈圖(遲滯密度法)··· 304 圖7-239：Min=0.5 次音速流場壁面馬赫數比較圖(遲滯密度法)··· 305 圖7- 240：Min=0.5 次音速流場馬赫數等值圖(遲滯密度法)··· 305 圖7-241：Min=0.5 次音速流場 Van Albada 及 SUPERBEE 限制函數法殘值

收歛比較··· 305 圖7-242：Min=0.675 穿音速流場隨參數 dm 值變化之壁面馬赫數分佈圖(遲 滯密度法)··· 306 圖7-243：Min=0.675 穿音速流場隨參數 k 值變化之壁面馬赫數分佈圖(遲滯 密度法) ··· 306 圖7-244：Min=0.675 穿音速流場隨參數 Mref 變化之壁面馬赫數分佈圖(遲 滯密度法)··· 307 圖7-245：Min=0.675 穿音速流場馬赫數等值圖(遲滯密度法)··· 307 圖7-246：Min=0.675 穿音速流場壁面馬赫數分佈圖(遲滯密度法)··· 307 圖7-247：NACA 0012 次音速流場 M∞=0.63，α=2°，dm=2 (SUPERBEE) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯密度法) ··· 308 圖7-248：NACA 0012 穿音速流場 M∞=0.80，α=1.25°，dm=0.8 (SUPERBEE)

(a)馬赫數等值圖；(b)翼表面壓力分佈(遲滯密度法)··· 308 圖7-249：NACA 0012 穿音速流場 M∞=0.85，α=1.0°，dm=1(Van Albada) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯密度法) ··· 308 圖7-250：NACA 0012 超音速流場 M∞=1.20，α=0.0°，dm=1(Van Albada) (a) 馬赫數等值圖；(b)翼表面壓力分佈(遲滯密度法) ··· 309

符號說明

a 音速 AP，AC 線性方程矩陣之係數 C B 方程式(3.48)中的係數 C 相鄰網格控容體心之位置 c 音速，溝槽圓突物之弦長 cp 等壓比熱 cv 等容比熱 d 定義在δPC方向上之向量 D 差分方程中的壓力係數 dm 密度限制函數的上限值 e,

ε

單位體積之內能 eˆ 單位向量 f 熱通量向量 C F 對流通量 D F 擴散通量 h 單位質量之熱焓 k 熱傳導係數 K 音速平方的倒數(1/γRT) M 馬赫數 m 質量流率 m 每個控容體之邊數 n 邊界單位法向量 p 絕對壓力 P 網格控容體心之位置 Pb 背壓

Pr 浦朗特常數 Po 全壓(停滯點壓力) P′ _{壓力修正量} Qφ 線性方程矩陣右邊之係數(源項) r 兩連續梯度比 R 通用氣體常數 S 微分方程式的源項 S_{φ 傳輸方程式的源項} S 面積向量(向外) f S 控容體f 面之面積向量(向外) t 溝槽圓突物之高度(厚度) T 絕對溫度 T Jacobian 矩陣 To 全溫(停滯點溫度) u

，

v 在卡氏座標系統x，y 方向之速度分量 V 速度場 V_∞ 遠場速度 w 權重因子 x，y 卡氏座標系統 PC δ 從點P 到 C 的距離向量 BD δ _{從點B 到 D 的距離向量}

ρ

密度 Ω 體積 μ 動黏滯係數 κ 熱傳導係數 T ‧ 應力張量 τw，τn 剪應力，垂直應力 ij

γ 等壓與等容比熱之比值(

γ

= cp/cv) f γ 通量限制子(flux limiter) Γ 擴散係數 αφ 混合因子 β 入口流向與入口邊界向量之夾角，震波之偏移角 θ 壁面之偏移角 Θ 壁面或對稱面與卡氏座標x 軸之夾角 φ 傳輸方程式之變數，如速度或溫度 , , P D U φ φ φ 上游、下游及上上游格點之

φ

值 ϕ 特徵變數 , , P D U ϕ ϕ ϕ 上游、下游及上上游格點之ϕ值 ηφ 鬆弛因子 ψ 通量限制函數 ∇ 梯度應算符號 Δ 變量 ΔV 網格體積 Δt 時間間隔 上標 o 前一次疊代之值 *，** 預測值，修正後之值 x，y，z 在卡氏座標系統x，y，z 軸上之投影量 ′ 修正值 ″，⊥ 平行，垂直邊界 下標 1，2 座標點之標示 B 所考量局部面的上游格點

C 相鄰網格格點 D 所考量局部面的下游格點 f 控容體f 面之值 in 入口 i，j 卡氏座標系統在x，y 軸上之單位向量 P 主要網格格點 U 所考量局部面的上上游格點 x，y 在卡氏座標系統x，y 軸上之投影量 φ 任意原始變數，如速度、溫度、密度或壓力 o 停滯狀態 ∞ 遠場狀態 簡寫字

BICG Bi-Conjugate Gradient

CBC Convective Boundedness Criterion CDS Central Difference Scheme

CLAM Curved-Line Advection Method

COPLA Cmbination of Piecewise Linear Approximation

CPU Computer Procedure Unit

CUBISTA Convergent and Universally Bounded Interpolation Scheme CUS Cubic Upwind Scheme

CV Control Volume

ENO Essentially NonOscillatory

EULER Exponential Upwinding or Linear Extrapolation Refinement

FCT Flux-Corrected Transport

FVM Finite Volume Method

HOT High Order Term

HUSG Hybrid UnStructured Grid

LODA Local Oscillation-Damping Algorithm LUDS Linear Upwind Difference Scheme

LUS Linear Upwind Scheme

MINMOD Minimum Modulus

MUSCL Monotonic Upwind Scheme for Conservation Law NVD Normalized Variable Diagram

NVSF Normalized Variable and Space Formulation PRIME PRessure Implicit Momentum Explicit PISO Pressure Implicit with Split Operator QUSG Quadrilateral UnStructured Grid

SHARP Simple High-Accuracy Resolution Program

SIMPLE Semi-Implicit Method for Pressure-Linked Equations SIMPLEC SIMPLE Consistent

SIMPLEM SIMPLE-Modified SIMPLER SIMPLE-Revision SIMPLEST SIMPLE-ShorTened

SOUCUP Second-Order Upwind-Central differencing-first-order Upwind

SMART Sharp and Monotonic Algorithm for Realistic Transport STOIC Second and Third Order Interpolation for Convection TUSG Triangular UnStructured Grid

TVD Total Variation Diminishing UDS Upwind Difference Scheme

UMIST Upstream Monotonic Interpolation for Scalar Transport QUICK Quadratic Upstream Interpolation for Convective Kinematics WACEB Weighted-Average Coefficient Ensuring Boundedness

第

1 章 緒 論

1.1. 簡介 由於近年來電腦硬體科技的快速發展，快速且高記憶體容量之資料處 理能力，使得計算流體力學的發展有長足的進步，由簡單的幾何形狀進步 到複雜外形之流場模擬計算，進而包含邊界層紊流之模擬計算，並針對不 同之雷諾數與馬赫數之流場發展適用之數值計算方法，致力於提昇流場計 算準確度及效率的研究均有顯著的成效。 在實際應用上除了馬赫數小於 0.3 之低速流場採用不可壓縮流解子

(incompressible flow solver)外，餘均必須使用可壓縮流解子(compressible flow solver)，然而在可壓縮流場中往往都包含了次音速(subsonic)、穿音速 (transonic)及超音速(supersonic)等單獨或交互存在之複雜流場，由於次音速 之統御方程為橢圓型式(elliptical type)與超音速之雙曲線型式(hyperbolic type)流場特性不同，故可壓縮流解子必須具備自動偵測流場特性進而調整 統御方程型式之能力，另對於震波之模擬除要獲得精確的位置及強度外， 亦要避免在梯度變化劇烈、陡峭之區域或震波前後之數值振盪產生，方能 獲得準確的流場解析。 大部份發展用於求解可壓縮流場的數值計算方法大致可分為兩類：密 度基底法(density based)及壓力基底法(pressure based)，不同基底法的使用均 有其特定的馬赫速範圍，前者常用於解高馬赫速流場；後者則常用於解低 馬赫速流場，然而若要發展全速流(all speeds flow)的求解方法，使其有能力 完全處理不同的雷諾數或馬赫速流場範圍問題之模擬計算，則必須瞭解不 同基底算則設計架構上之困難點，尤其是壓力在可壓縮流場中所扮演的角

色。密度基底法傳統上都使用非穩態形式之Navier-Stokes 或 Euler 方程式，

這些方法皆使用密度作為主要變數，以連續方程式作為求解密度之傳輸方 程式，其壓力則藉由狀態方程式來求得，通常應用在解高速可壓縮流場

[1-2]，然此方法用於解不可壓縮流或低馬赫數流場時，因為密度的變量非 常小或幾近於常數，則其與速度之偶合關聯性會變得相當薄弱，因而連續 方程式無法成為求解密度之傳輸方程式且會抑制速度場，因此必須設計一 種機制能透過壓力場將連續方程與動量方程偶合起來，解決的方法為虛構 狀態方程式或在連續方程式中加入人工可壓縮因子[3-8]；壓力修正法是將 由動量方程式導出的壓力變量與速度變量的關係式代入連續方程式中，推 導得到壓力修正方程式進而解出壓力場，其密度則藉由狀態方程式來求 得，對於高速可壓縮流場而言，其速度的變量遠小於速度值本身，而密度 的變化則非常大，此時壓力的變量主要來自於密度變量，故必須將密度變 量轉變為壓力變量代入壓力修正方程式中，方能求解可壓縮流場。 由以上的討論可以知道，任何數值方法要有能力同時處理不可壓縮流 及可壓縮流流場不同特性之統御方程型式，壓力則扮演著相當重要的角 色，其影響速度與密度來滿足連續方程式，總之，經由使用所謂的虛擬或 人工可壓縮因子之技巧，已有多種密度基底法可模擬全速流流場之發展被 提出[3-4,8]，這些方法遭遇到的困難為如何避免艱難的矩陣求解效率而降低 了收歛率，為了克服這些問題及確保在整個全速流場範圍內均能收歛，則 致力於預先處理艱難矩陣的求解方法研究[4,8-12]。在另一研究領域許多學 者則進行壓力基底法延伸至高馬赫數可壓縮流之研究，由於壓力基底法之 壓力修正方程式本質上為橢圓型式之特性，故不適用於超音速或穿音速流 場之雙曲線型式或混合系統之特性，如果沒有修正而企圖直接應用將造成 計算過程不穩定或甚至發散現象，已有不少學者採用壓力基底法進行可壓 縮流計算，配合使用交錯網格(staggered grid)法[13-14]、重置變數(collocated variables)法[15-20]或多重網格法等已有不同程度的貢獻與成就。壓力法不 論在低速不可壓縮流或高速可壓縮流皆保有速度與壓力很強的偶合關聯 性，經過最近數十年的發展已被廣泛用於工業界設計與分析[21]，對使用者 而 言 ， 在 面 對 複 雜 流 場 分 析 之 工 具 需 求 為 簡 便(friendness)與功能強大 (robustness)，壓力基底法已受到關注並被期待致力於計算方法的持續發展，

俾利成為更有用的工程設計分析工具。

壓力法在流場計算之數值方法建構方式，依求解變數型式可分為三 類，一為以原始變數(primitive variables)求解穩態(steady state)或非穩態 (unsteady state)之統御方程式，為 Issa 及 Lockwood 等人 [22]所提出，其延 用求解不可壓縮流之方式，以連續方程式推導壓力修正方程式並作部份修 訂，將密度的變量轉換為壓力變量，在求面上質量通量(mass flux)時係將密 度及速度分別處理，此為壓力法解可壓縮流場之標準方式。一般對於面上 密度之計算均採用一階上風差分法來近似，雖較穩定但其對震波之解析較 模糊，為了得到較佳之震波解析，亦有人使用中央差分近似法[15-17]，但 易造成計算過程不穩定甚至發散；而 Karki 等人[14]及 Demirdzic 等人[16] 使用中央差分混合一階上風差分近似法，其混合因子則介於 0 與 1 之間， 主要的目的無非是要兼顧增加解的精確度及疊代過程的穩定性；Moukalled 與 Darwish[23]則使用與對流變數相同之 SMART 高階限制函數差分法來處 理面上密度。第二種為以守恆變數(conserved variables)求解非穩態之統御方 程式，其求解方式與不可壓縮流類似，唯在求解動量方程式時，使用遲滯 密度或遲滯壓力作為人工消散項代入係數中，藉以穩定此計算方法，McGuir

與Page[24]及 Lien 等人[15]分別引用 Wornom 等人[25]提出之”遲滯壓力”及” 遲滯密度”概念，結合 SIMPLE 壓力修正法，成功地應用於可壓縮流計算。 第三種為特徵變數(characteristic variables)，基本上仍是以原始變數求解穩態 或非穩態之統御方程式，只是運用 ENO 或 TVD 的技巧使用特徵變數並結 合Riemann 解子(solver)來計算面上的通量[26-27,48]，後續衍生性的研究在 文獻報告中並不多。 在處理包含對流(convection)與擴散(diffusion)性質的流場模擬計算問 題時，對於擴散項而言，採用二階準確之中央差分法幾乎是無異議的選擇， 然而對於對流項而言，則有各種不同的處理策略，而其最重要的考量因素 就是要能夠具有界限性(bounded)、精確性及穩定性的模擬方法。若從物理 上而言，當對流項使用界限性差分法離散時，則其近似值絕不會超越局部

解的極大或極小值，即可抑制非物理性或過與不及的數值振盪現象。以一 般流體問題純量性質的傳輸方程式為例，如多相流比率、紊流動能、質量 比率等，其界限性的重要性就非常清楚，如果在К-ε 紊流模式計算時產生負 的紊流消散值，則會導致負的紊流黏滯係數，通常會造成災害性的影響。 自從發展一階上風差分法用來離散對流通量之後，由於不能滿足應用 上的需求，故研究學者開始進行高階準確法之發展，由於「精確度」恰與 「穩定及界限性」相互矛盾及抗衡，可謂魚與熊掌不可兼得，必須在兩者 之間作適當地折衷，然而為了要成功地解決精確度的問題，則被迫必須面 臨數值穩定及物理界限性之困難議題。為了克服這些源自於準確度、穩定 性及界限性等三方面需求相互矛盾的問題，為了維持相當的準確度且抑制 振盪現象，不同的高階界限差分法組合程序被提出，最熱門的兩支家族為 全變量消減(TVD)法[28-34]及正規化變數(NV)法 [35-39]，這些方法均使用 已知之局部剖面解及不同的限制函數來調整對流項離散之差分法，其所使 用之限制函數型式係由一維解的局部分佈剖面分析所推導得到，通常參照 到的網格除了網格面相鄰之網格外，尚需一更上游之網格。許多研究者持 續致力於發展可準確模擬對流佔優流場之高階準確解析法，此類之方法大 概 可 分 為 通 量 混 合 法(flux blending method) 及 組 合 式 通 量 限 制 函 數 法 (composite flux limiter method)兩種類型，這兩種方法均企圖能在不影響準確 度之情況下，獲得抑制數值振盪之流場預測解。 通量限制函數之發展從1985 年由 Van Leer [40]首先提出，係以守恆變 數作為通量限制函數計算之關係變數，即通量限制函數為正規化守恆變數 之函數。後來亦有學者[18,37-38,41-46]以原始變數或特徵變數作為通量限制 函數計算之關係變數，發展高階單調之對流項差分方法，期獲得高精確且 無數值振盪之模擬計算法，大部份均以密度基底法求解。使用特徵變數作 為通量限制函數計算之關係變數在壓力基底法則不多見，例如 Kobayashi

與 Pereira [26] 以特徵變數內插方式代入 Roe 的 Riemann 近似解子用以計

方程之偶合；而在Issa 與 Javareshkian [27]的研究中，對流通量則使用與 Yee 等人[47]相同的方式並使用具 TVD 之 MINMOD 通量限制函數。Shyy 與 Thakur [48]則明確指出 TVD 差分法通常不完全適用於壓力基底法，其原因 有二：首先因為壓力基底法為動量方程與連續方程(壓力修正方程)分別求解 之順序疊代(sequential-iteration)求解法，其動量方程式中之壓力梯度係以顯 項來處理並將其置於源項中，不像密度基底法是將變數之梯度視為通量 (flux vector)的一部份，偶合連續方程、動量方程與能量方程同時求解統御 方程；第二個原因為缺乏作為TVD 限制函數法通量限制子(flux limiter)計算 基礎的流場局部特徵值之定義。雖然如此，仍有不少學者嚐試以壓力基底 法建構各種對流差分方法與機制，期能獲的類似密度基底法TVD 之精確解。 使用壓力基底法在建構限制函數時，一般均採用守恆或原始變數來計 算梯度比，此與適用於可壓縮流密度基底法所使用以特徵變數為基礎之限 制函數則有所不同。由於穿過震波時只有 Riemann 變數之變量最小，相較

之下守恆或原始變數之變量則非常大，經由Mulder 與 Van Leer [49]及 Lin

與Chieng [39]密度基底法的數值分析證明，至少在一維流場時使用 Riemann 變數確實可以獲得最高之精確度，Shyy 與 Thakur [48]嚐試使用特徵變數基 礎之限制函數，以標準之 SIMPLE 壓力修正法，應用結構性非正交網格、 重置變數之有限容積法求解可壓縮流，已可獲得不錯的結果，唯尚無發現 用於無結構性網格之研究文獻。 流場計算需要將整個流場離散為細小的網格，而網格產生的種類大致 可分為兩類：結構性網格(structured grid)與無結構性網格(unstructured grid) (參考圖 1-1、1-2、1-3)，各有其優劣點，結構性網格具有優先的方向及自 動指標等方便引用的特性；可是遇到比較複雜的幾何外形時，則會遭遇許 多限制及困難；早期大部分之流場解子均為使用結構性網格之有限體積法 (FVM)，最近已有許多研究發展可應用於無結構性網格之有限體積法，近年 來，無結構性網格產生技巧已被成功地應用在處理複雜幾何外形的流場計 算，由於其能用任意形狀的多邊形來產生複雜幾何外形之網格，而且可以 很容易地作局部加密，故其較結構性網格具有更大的彈性，但其在計算上

較複雜且需要較大的記憶空間及計算速度，然此問題則因電腦的快速發展 而獲得解決。 本文主要的目標為發展一套全速流流場解子，以有限體積法(FVM)、 應用重置變數(collocated variables)之任意邊形無結構性網格(unstructured grid)，引用可預測二維複雜幾何外形不可壓縮流的 SIMPLE 壓力修正方法 [50]，加以延伸至可壓縮流，以原始變數及守恆變數兩種型式求解統御方程 式並企圖嚐試以各種 TVD 或正規化變數(NV)之對流項高階準確差分法進 行數值方法研究。期能獲得具有高階準確、無數值振盪、穩定及收歛性良 好之全速流(all- speed flows)流場計算之數值方法與程式。

1.2. 文獻回顧 1974 年 Khosla 與 Rubin[51]提出一種空間二階準確中央差分近似法之 修改型，以一階隱項差分法配合顯項之二階差分法來離散空間一階導數微 分項，當收歛時則自然回復至二階準確中央差分法，此種方式可確保係數 矩陣為對角佔優(diagonal dominate)，為無條件穩定法，後來的學者均引用 此處理方式應用於有限容積法一階導數對流通量項之離散，並將這種處理 方式稱之為遲緩修正法(deferred correction method)。

1979 年 Hafez 等人[52]以求解守恆形式的全勢能(potential)方程式為基 礎 ， 藉 由 稍 微 修 改 在 超 音 速 區 域 的 密 度 來 產 生 人 工 黏 滯 項(artificial viscousity)，使用上風差分法處理修正密度，再以中央差分法將網格中心之 修正密度內插至網格面上，可解析流體流經圓柱及NACA 0012 翼型之二維 穿音速流場。 1983 年 Rhie 等人[53]提出控容面上速度計算之內插方式，將速度變量 以壓力梯度變量之型式導入壓力修正方程式中。 1984 年 Sweby[54]提出在一階上風差分法加上一個反擴散限制通量之 策略，推導TVD 及二階準確 TVD 特性之單調界限準則及具有二階準確 TVD

界限性質之通用限制函數(limited function)，獲得二階高階準確無數值振盪

之顯式純量微分法，同時對不同的限制函數，包括 Van Leer、Roe、

Chakravarthy 與 Osher 等人所提出之限制函數，進行理論推導與模擬計算之 比較，其使用通量限制函數ψ為 Y 軸；局部守恆梯度比(gradient ratio) r 為 X 軸之二維圖形表示方式，可清楚顯示通量限制函數(flux limiter function)是否

具有TVD 及二階準確 TVD 特性，即為後來常被人引用之 Sweby 的「TVD

關係圖」或「Φ-限制函數圖」。經由不同學者所提出之限制函數模擬計算

結果比較分析，Roe 的限制函數在線性情況下有優異的表現，Van Leer 的限

制函數表現與 Roe 相當並具有更好的可靠度，由此說明限制函數表現的優 略與流場系統的守恆律有關。 1984 年 Wornom[55]使用兩點式中央差分法來離散通量導數項，並應用 遲滯密度(retarded density)的概念在超音速點加入消散項，以一維漸縮-漸擴 噴嘴(C-D nozzle)測試此方法可執行次音速及超音速流之計算且有良好的效 率，其計算結果雖無數值振盪，但對於震波強度的捕捉則顯不足。 1986 年 Rhie[13]以壓力法使用多重網格求解 Navier-Stokes，可用於解 低速不可壓縮流及高速可壓縮流之全速流(all speeds flow)。

1986 年 Wornom 與 Hafez[25]使用遲滯密度(retarded density)或遲滯壓力 (retarded pressure)作為人工消散項，配合原解次音速流場之兩點法(two-point method)求解有震波之穿音速流場，其使用之方式藉由遲滯密度代入連續方 程式中以改變通量係數；而當使用遲滯壓力作為人工消散項時，於動量方 程式中以遲滯壓力取代真實壓力，只有在利用狀態方程計算密度時才使用 真實壓力。經由比較兩種方式計算結果，發現在使用遲滯壓力可得到不錯 的結果，而在使用遲滯密度時則有輕微抺平現象。若同時使此兩種人工消 散項，則在震波處會產生較嚴重的抺平現象。 1988 年 Gaskell 與 Lau[56]以正規化變數配合單調之限制條件推導出對 流界限準則(CBC)，將其顯示在 NVD 圖上，可清楚檢視各種差分法是否具 CBC 區域內提出以 QUICK 為架構配合曲率補償之三

階準確 SMART 對流項差分法，其具有界限性且在梯度變化劇烈處有不錯 的表現， 1988 年 Leonard[37]以用於穩態對流方程有不錯特色，然而在劇烈階梯 形剖面對流模擬時，則會產生非物理性的超估及數值振盪，尤其是在非線 性流場此現象會導致嚴重的問題之三階準確之 QUICK 上風差分法為基 礎，提出將網格面上對流變數使用正規化變數方式化為其相鄰網格及更上 游網格正規化變數之函數，並發展出單調的解析準則而不會犧牲解的精確 度。藉由具界限之SHARP 非線性正規化變數關係式應用於包含有薄剪力層 或混合層、震波及其他如前所述現象之流場模擬，SHARP 為顯性守恆控容 通量之差分公式，可用於一至三維之橢圓、拋物線、雙曲線或混合之流場 區域，經由斜階梯對流函數之測試並與一階上風法、二/三階非界限上風法 比較，有不錯的結果。 1988 年 Majumdar[57]探討鬆弛參數在非交錯網格法的流場計算中之動 量內插方法，在計算面上速度時可用相鄰網格中心點之速度值(無壓差效應 時)以線性內插之方式，求出之值再與該面相鄰網格中心點之壓差效應結 合，來近似面上速度，結果都能得到穩定的唯一值。 1988 年 Peric 等人[58]以有限容積之數值方法來作交錯網格與重置網格 在二維不可壓縮流流場計算之比較，計算結果在收歛率、鬆弛參數之依附 性、計算上的能力與準確性等方面，一般而言兩者均相當。而在某些情況 重置網格法收歛較快，如果要考慮作非正交或多區塊網格時，重置網格較 具優越性。對於大部份的不可壓縮流場計算，建議最佳的速度鬆弛因子為 0.8，壓力鬆弛因子為 0.3。

1988 年 Zhu 與 Leschziner [59]提出 QUICK 及一階上風混合之局部振盪 緩和法(LODA)求解不可壓縮穩態流場，混合因子以類似 Gosman 與 Lai 的 方式決定，由流場依動態調整。

序及Roe 的通量函數，發展格心(cell-center)及格點(mesh-vertex)兩種高階上 風有限體積差分法求解Euler 流場，任何高階單調準確法皆以一階上風為基 礎，再加上消散項，當在極值產生時，為了保持單調性避免振盪產生，就 必須建立一套機制由高階上風回歸到一階上風法。其使用二階上風法並在 二次項中加入一個限制因子來重建網格面上之近似值，再計算通量，限制 因子之值介於 0 與 1 之間，此限制因子的計算方式為先行使用格點共用之 網格中心面積平均值求得該格點上之近似值，利用流動方向及梯度比求得 該格點之限制因子，再由局部網格所共構格點之限制因子中找出極小值， 而此值即作為網格面上二階上風重建法之限制因子，以確保網格面上之近 似值具界限性，其值介於局部網格與相鄰網格中心之極值範圍內。梯度之 計算策略採用網格格點所包含之所有相鄰網格中心連線之積分路徑。 1989 年 Karki 等人[14]使用結構性之交錯網格，以壓力法求解任意外形 之黏性全速流流場，為了要處理不可壓縮流與可壓縮流，在壓力與密度相 依關係中，以壓力作為優先之主要變數，壓力場由可壓縮形式之SIMPLER 算則中解出，壓力解出後再求出密度。 1990 年 Peraire 等人[61]，指出結構性網格雖然具有高效率的計算方 法，但對於複雜幾何外形，則必須作座標曲線轉換及採用多區塊網格產生 法，須花費相當多的時間精力去處理；而無結構性網格則可以很快地直接 在物理空間產生複雜幾何外形之計算網格，並具有容易作局部加密及網格 調適等功能，非常有彈性，由於需要額外儲存網格之點、線、面及其相鄰 網格之關係，故計算程序較複雜且需要較大的記憶體及CPU 時間。非結構 性網格之離散方程式均使用積分方式之有限體積法或有限元素法。

1990 年 Mcguirk 及 Page [24]在解準一維非穩態的歐拉(Euler)方程式

時，使用不可壓縮流之 SIMPLE 壓力修正算則，將遲滯壓力同時用於動量

方程式及壓力修正方程式，以獲得雙曲線型式之統御方程式特性，遲滯壓 力轉換為真實壓力之值係由局部馬赫數決定，計算結果對於震波的捕捉及 計算時間上都有不錯的表現。另亦執行紊流軸對稱衝擊噴流之模擬計算，

結果可以精確地捕捉正震波。 1991 年 Leonard [38]提出一種新的界限對流差分法 ULTIMATE 來模擬 暫態對流傳輸方程，其以顯性守恆的控制體積差分公式為基礎，加入一種 時間平均通用正規化限制函數，其目標是要獲得一個健全的方法，能夠得 到精確單調的階梯型函數解而不會產生變形。首先以標準的測試問題，即 以方波、正弦函數波及半橢圓形波等形函數，使用暫態純對流傳輸方程各 別對一階上風差分法、二階中央差分法、二階上風差分法、Fromm 差分法 及三階以上之差分法等11 種不同線性差分方法進行測試，模擬結果除了一 階上風差分法因具有強的擴散性而可獲得單調無振盪的解外，其餘方法均 會產生數值振盪現象。另亦以相同問題對 MINMOD、Chakravarthy 與

Osher、MUSCL、Van Leer 的 CLAM、Roe 的 SUPERBEE、Super-C、Hyper-C

等7 種非線性”震波捕捉”或”TVD”之差分方法進行測試，除了 Super-C 外，

餘均滿足 Sweby 所提出的 TVD 準則，模擬結果就整體性能而言，MUSCL

有不錯的表現，Super-C、Hyper-C 對方形波之表現最佳，正弦波則以 Super-C

及 SUPERBEE 之表現最佳，半橢圓形波則以 MUSCL 及 CLAM 之表現最

佳。Leonard 提出 ULTIMATE 之策略主要是將面上空間正規化變數差分法 改以時間平均正規化變數差分法，經由測試結果顯示，對於高階之數值振 盪問題有不錯的改善，可應用於任何高階差分法。

1991 年 Lin 與 Chieng [39]使用密度基底法，以 Leonard 之正規化變數 (NV)概念並以特徵變數為基底，發展一新的通量限制函數，使對流項面上 之近似差分法可具有三階之準確度，其發展之通量限制函數為上風基底 MUSCL 型式之改良型，延伸自 SMART 與 SHARP 差分法之控容體通量公 式概念。以守恆變數(conserved variables)、原始變數(primitive variables)及特 徵變數(characteristic variables)等三種不同基底之通量限制函數，驗證一維震 波管、二維斜震波階梯渠道(3 及 10 馬赫，震波角 59 度)及 NACA 0012 翼 型(0.8 馬赫，攻角 0 度)等情況，結果採用特徵變數基底之通量限制函數較 其他兩種方式有較小的數值振盪現象，且其精確度與收歛性均優於二/三階