碩士論文中華大學

(1)

中華大學碩士論文

題目：An LU-SGS upwind scheme for the 3D incompressible flow

系所別：機械與航太工程研究所碩士班學號姓名：M09108024 歐陽弘道指導教授：蔡永培博士

中華民國九十四年十一月

(2)

(3)

(4)

(5)

(6)

(7)

中文摘要

本文主要的研究目的乃建立起一三維數值模擬程式，藉以分析不可壓縮流之流場。數值方法的選取在時間離散方面，採用 LU-SGS (Lower-Upper Symmetric-Gauss-Seidel)隱式解法以增強數值之穩定性及加速程式之收斂，在空間離散方面，無黏性通量採用 UP-WIND 曲型算則，黏性通量則是採中央差分法。

文中選取了兩個典型不可壓縮流之測試範例-空穴流、背向階梯流，並與文獻所載之結果相驗證及比對其數值模擬解。

(8)

Abstract

A numerical simulation program to analyze 3D incompressible flow has been developed. Several numerical solution schemes were used to solve the governing Navier-Stokes equations. The LU-SGS ( Lower-Upper Symmetric-Gauss-Seidel ) implicit scheme was adopted in time to enhance numerical stability and to achieve fast convergence. UP-WIND hyperbolic conservation was used in space to obtain high resolution for non-viscous flux, while central-difference method was applied in space for viscous flux.

Two classical incompressible flows, cavity flow and backward-step flow, were tested to validate the numerical solution program. Good agreement between the numerical simulation and the published experimental data in the literature was obtained.

(9)

誌謝

本文得以順利完成，首要感謝我的指導老師蔡永培博士於我研究生涯中不辭辛勞的悉心指導，得以令我在學術研究及待人處事中獲益良多，此外承蒙口試委員黃和順博士與本校老師蔡博章士熱心指導，

並惠賜許多寶貴意見及建議，使本文更臻理想，在此深表感謝。

在我求學期間，承蒙學長黃耀德、鐘荃福、韋立興、董維龍等人給我的指導，同學徐明照、林長毅等人給予我的寶貴意見，讓我在研究過程中更具信心，伴我渡過難忘的求學生涯；此外，在研究過程中，

學弟周義騰、周中祺等人幫助我收集資料，並給予我電腦軟體上的幫助，令我在研究過程中無後顧之憂，特此感謝。

最後，特別感謝我最敬愛的父母歐陽憶農先生蔡淑惠女士及我的妹妹歐陽婉麗與女友雅芳，在求學過程中全程陪伴著我，不時給予我支持與鼓勵，讓我擁有無比的信心得以順利完成學業；願將本文獻給關心我的每一個人，共同分享我的喜悅與成果。

(10)

中文摘要 ……… Ⅰ 英文摘要 ……… Ⅱ 誌謝 ……… Ⅲ 目錄 ……… Ⅳ 圖目錄 ……… Ⅵ 符號說明 ……… XI

第一章緒論 ……… 1

1.1 研究動機 ……… 1

1.2 文獻回顧 ……… 2

1.3 研究方法與目的 ……… 5

1.4 章節介紹 ……… 6

第二章基本理論與統御方程式 ……… 8

2.1 前言 ……… 8

2.2 基本物理特性 ……… 9

2.3 統御方程式 ……… 10

2.4 座標轉換 ……… 12

(11)

第三章數值方法 ……… 19

3.1 空間離散 ……… 19

3.2 迎風算則 ……… 20

3.3 時間離散 ……… 21

3.4 LU-SGS ……… 24

3.5 網格建立 ……… 27

3.5.1 幾何公式格點產生法 ……… 29

3.6 邊界條件 ……… 32

3.6.1 空穴流 ……… 33

3.6.2 背向階梯流 ……… 34

第四章結果與討論 ……… 36

4.1 三維穩態流之應用 ……… 36

4.2 空穴流 ……… 36

4.3 背向階梯流 ……… 38

第五章結論與未來展望……… 39

5.1 結論 ……… 39

5.2 未來展望 ……… 40

參考文獻 ……… 42

(12)

圖目錄（附錄）﹕

三維空穴流部份：

【圖 5.1】三維空穴流之流場示意圖及其邊界條件……… 44

【圖 5.2】三維空穴流計算外型及網格分布圖……… 44

【圖 5.3】三維空穴流等壓力圖(Re=100)……… 45

【圖 5.4】三維空穴流等 u 速度圖(Re=100)……… 45

【圖 5.5】三維空穴流等 v 速度圖(Re=100)……… 46

【圖 5.6】三維空穴流等 w 速度圖(Re=100)……… 46

【圖 5.7】三維空穴流流線-等 U 速度圖(Y 方向切面) (RE=100)…… 48

【圖 5.8】三維空穴流流線-等 W 速度圖(Z 方向切面) (RE=100)…… 49

【圖 5.9】雷諾數為 100 時三個中剖面之向量流線圖……… 50

【圖 5.10】雷諾數為 100 時之等壓力圖……… 51

【圖 5.11】雷諾數為 100 時之等渦度圖……… 52

【圖 5.13】雷諾數為 400 時之等壓力圖……… 54

【圖 5.14】雷諾數為 400 時之等渦度圖……… 55

【圖 5.15】三維空穴流等壓力圖(Re=1000)……… 56

【圖 5.16】三維空穴流等 u 速度圖(Re=1000)……… 56

【圖 5.17】三維空穴流等 v 速度圖(Re=1000)……… 57

(13)

【圖 5.18】三維空穴流等 w 速度圖(Re=1000)……… 57

【圖 5.19】三維空穴流流線-等 U 速度圖(Y 方向切面)(RE=1000)… 59 【圖 5.20】三維空穴流流線-等 W 速度圖(Z 方向切面)(RE=1000)… 60 【圖 5.21】雷諾數為 1000 時三個中剖面之向量流線圖……… 61

【圖 5.22】雷諾數為 1000 時之等壓力圖……… 62

【圖 5.23】雷諾數為 1000 時之等渦度圖……… 63

【圖 5.25】雷諾數為 3200 時之等壓力圖……… 65

【圖 5.26】雷諾數為 3200 時之等渦度圖……… 66

【圖 5.28】雷諾數為 5000 時之等壓力圖……… 68

【圖 5.29】雷諾數為 5000 時之等渦度圖……… 69

【圖 5.31】雷諾數為 7500 時之等壓力圖……… 71

【圖 5.32】雷諾數為 7500 時之等渦度圖……… 72

【圖 5.33】三維空穴流等壓力圖(Re=10000)……… 73

【圖 5.34】三維空穴流等 u 速度圖(Re=10000)……… 73

【圖 5.35】三維空穴流等 v 速度圖(Re=10000)……… 74

【圖 5.36】三維空穴流等 w 速度圖(Re=10000)……… 74

(14)

【圖 5.37】三維空穴流流線-等 U 速度圖(Y 方向切面) (RE=10000)… 76

【圖 5.38】三維空穴流流線-等 W 速度圖(Z 方向切面) (RE=10000)… 77

【圖 5.40】雷諾數為 10000 時之等壓力圖……… 79

【圖 5.41】雷諾數為 10000 時之等渦度圖……… 80

【圖 5.42】通過中心軸的速度分布比較(Ex.- UP-WIND)……… 81

【圖 5.44】通過中心軸的速度分布比較(WENO3 - UP-WIND)……… 82

【圖 5.48】不同算則通過中心軸的速度分布比較圖……… 84

(15)

【圖 5.58】人工壓縮因子對收斂歷線影響圖……… 89

【圖 5.59】人工壓縮因子對收斂歷線影響放大圖……… 89

【圖 5.60】時間離散方法對收斂歷線影響圖(Lusgs)……… 90

【圖 5.61】時間離散方法對收斂歷線影響圖(Lussor)……… 90

三維背向階梯流部份：

【圖 5.62】三維背向階梯流之流場示意圖及其邊界條件……… 91

【圖 5.63】三維背向階梯流計算外型及網格分布圖……… 91

【圖 5.64】三維背向階梯流等壓力圖(Re=100)……… 92

【圖 5.65】三維背向階梯流等壓力圖(Re=100)……… 92

【圖 5.66】三維背向階梯流等 u 速度圖(Re=100)……… 93

【圖 5.67】三維背向階梯流等 u 速度圖(Re=100)……… 93

【圖 5.68】三維背向階梯流等 v 速度圖(Re=100)……… 94

【圖 5.69】三維背向階梯流等 v 速度圖(Re=100)……… 94

【圖 5.70】三維背向階梯流等 w 速度圖(Re=100)……… 95

【圖 5.71】三維背向階梯流等 w 速度圖(Re=100)……… 95

【圖 5.72】三維背向階梯流流線圖(Y 方向切面) (Re=100)………… 96

(16)

【圖 5.73】三維背向階梯流流線圖(Z 方向切面) (Re=100)………… 97

【圖 5.74】雷諾數為 100 時中剖面之等速度圖……… 98

【圖 5.75】雷諾數為 100 時中剖面之向量流線圖……… 98

【圖 5.82】雷諾數－流離接觸點關係圖………102

(17)

符號說明

X 物理域直角座標系統之 X 方向 Y 物理域直角座標系統之 Y 方向 Z 物理域直角座標系統之 Z 方向 ξ 計算域之座標軸

η 計算域之座標軸 ζ 計算域之座標軸 Δ 微變量

ΔX X 座標上之空間間隔 ΔY Y 座標上之空間間隔 ΔZ Z 座標上之空間間隔 Δξ ξ座標上之空間間隔 Δη η座標上之空間間隔 Δζ ζ座標上之空間間隔

Δt 時間間隔 u X 方向之速度 v Y 方向之速度 w Z 方向之速度

(18)

Q 物理域之無因次守恆變數向量 E X 方向之無因次非黏性通量向量 F Y 方向之無因次非黏性通量向量 G Z 方向之無因次非黏性通量向量

E

v X 方向之黏性通量

F

v Y 方向之黏性通量

G

v Z 方向之黏性通量

Qˆ

計算域之守恆變數向量

Eˆ

ξ 方向之非黏性通量向量

Fˆ

η 方向之非黏性通量向量

Gˆ

ζ方向之非黏性通量向量

Eˆ

v ξ 方向之黏性通量

Fˆ

v η 方向之黏性通量

Gˆ

v ζ方向之黏性通量 V 控制體積

ˆA Eˆ 之 Jacobian 矩陣

ˆB Fˆ 之 Jacobian 矩陣

Cˆ Gˆ

之 Jacobian 矩陣 λ 特徵值

(19)

R

與矩陣 Aˆ 、 Bˆ 、 Cˆ 之特徵值相對應之右特徵向量所 構成的矩陣

R

1

與矩陣 Aˆ 、 Bˆ 、 Cˆ 之特徵值相對應之左特徵向量所 構成的矩陣

r^s 第 s 特徵值所構成之右特徵向量 l^s 第 s 特徵值所構成之左特徵向量 Re 雷諾數

υ 運動黏度

ij 表層黏滯剪應力張量 μ 動力黏度

l 表層流運動黏滯度

t 表紊流渦漩黏滯度

T

表層流與紊流的黏滯度結合在一起之黏滯度

（ _T _l _t ）人工壓縮因子

P

參考壓力

V

x 參考速度 L 參考長度 P 壓力

(20)

□* 表□的無因次化 h 座標轉換之 Jacobian

LU-SSOR 之參數值

Rˆ

殘餘量（residual），亦即通稱之右手項（right hand side）

D

向前差分運算子

D

向後差分運算子

Aˆ

i 表

Aˆ

_i矩陣處理過後之矩陣特徵值為非負

Aˆ

i 表

Aˆ

_i矩陣處理過後之矩陣特徵值為非正

Aˆ 表矩陣 Aˆ 之譜半徑（spectral radius of Jacobians）

Bˆ 表矩陣 Bˆ 之譜半徑（spectral radius of Jacobians）

Cˆ 表矩陣 Cˆ 之譜半徑（spectral radius of Jacobians）

A

ˆ 表矩陣 Aˆ 之特徵值

Bˆ

表矩陣 Bˆ 之特徵值

Cˆ

表矩陣 Cˆ 之特徵值 c 聲速

A

m 擾動之強制振幅

( )

G y

高斯分佈函數

f

m 頻率

(21)

m 隨機相位

S

ij 表應變率張量

R

ij 表雷諾應力

ij 表 Kronecker delta 函數

R

kk 表法向應力

(22)

第一章緒論

1.1 研究動機

早期的計算流體力學(Computational Fluid Dynamics)研究主要是借助於理論分析與實驗，然而傳統的理論分析方法由於有許多假設與簡化，所以其能解決的問題通常有限。隨著近代電腦硬體急速的發展，速度的不斷提升，使得每一步運算的時間大量的減少，計算流體力學所能解決問題的尺度與複雜度也逐漸加大，時至今日，計算流體力學已成為學界研究流體力學的主要利器之一，與理論流力和實驗流力構成現代研究流體力學之三大主流。

計算流體力學的內容主要是流體力學、數學、數值方法及電腦運算等的整合，而他的應用範圍也非常之廣泛如機械，航太、汽車、船舶、土木、化工、醫工、電工、軍事、大氣與海洋等均涵蓋在內，如飛機與汽車之外形設計，各類引擎燃燒室及冷凍空調系統設計，空氣及水污染物擴散預測，建築結構物受風及水的影響，氣候及溫度的改變，原子彈爆炸，心臟與血管的血流流動，高速火車進出隧道的噪音問題等，都可利用計算流體力學來研究與解決。

(23)

1.2 文獻回顧

奈維爾－史托克方程式為描述流體流動之最廣泛之方程式，各種不同型態流場的結構如震波、渦漩、邊界層、分離流及紊流等物理現象皆可由求解奈維爾－史托克方程式而獲得。在描述流體流動之偏微分方程式，大都是採用雷諾平均奈維爾－史托克（Reynolds Averaged Navier－Stokes；簡稱 RANS）方程式來模擬。RANS 方程式為目前廣泛使用的數學模式，其紊流結構是由所謂的紊流模式（turbulence model）來分析，但可惜的是，到目前為止，尚未能發展出可適用於任何流場且令人滿意的紊流模式，因而紊流模式的發展一直被視為計算流體力學最重要的領域之一，而常見的紊流模式之介紹與各種數學模式之推導細節可參考 Anderson et al.（1984）[1]。

不可壓縮與可壓縮奈維爾－史托克方程式之間最大的不同點，在於不可壓縮流的連續方程式中缺乏對時間的微分項，因此，「滿足質量守恆方程式」，則成為求解不可壓縮奈維爾－史托克方程式最重要的關鍵。若以物理觀點看來，壓力波的速度在真實不可壓縮流埸中為無窮大，此一壓力波的橢圓行逕，乃是不可壓縮流的主要物理特性之一，而在奈維爾－史托克方程式中，壓力場也是解的一部分，但不可壓縮流之連續方程式中並無壓力項，故壓力無法由統御方程式中求

(24)

取，因此，利用數值模擬方法求解流場時，最難處理的即為壓力項之未知，也正是如此，為了解決動量方程式中之壓力項的困擾，許多的求解方式則因運而生。

早期是用流函數－渦度法，此法對於二維流場問題的處理是利用交微分的方式，以消去兩動量方程式中之壓力項，然後產生一組渦度傳輸方程式，再引入二維之流函數以配合求解。另外，渦度－速度法則是從渦度傳輸方程式出發，將渦度傳輸方程式取漩度（curl）運算並引用連續方程式而得到一組速度的 Poisson 方程式。但是上述這種方法對於求解三維流埸時都有其困難處，並不容易解決。

而在以壓力與速度當作變數的原始變數法中，質量守恆及其與壓力的關係必須適當的掌控才可達到計算的效率，因此，發展出了許多利用原始變數來求解不可壓縮奈維爾－史托克方程式的技巧與研究，玆介紹如下：

MAC 方法（Marker－And－Cell method）：

此法是由 Harlow ＆ Welch（1965）[2]所發展出，其出發點是將壓力當做影射參數以滿足連續方程式，首先是藉由對動量方程式作散度（divergence）運算以得到壓力的 Poisson 方程式，這種必須直接先求解壓力的 Poisson 方程式的方法，其實用性僅限於簡單的二維流場計算，而對於三維流場的計算，則需採用疊代方式求解才是最佳

(25)

的選擇，因 MAC 方法必須在每一計算時步內就得先求得正確的壓力值以滿足無散度速度場的要求，此種要求的限制，會嚴重的降低計算效率。畢竟對穩定流而言，當流場達到穩定解暨滿足收斂需求時，正確壓力場的需求才具有意義。

SIMPLE 方法（Semi－Implicit Method for Pressure－Linked Equations）：

因為顯式法存在著時間步伐的先天性嚴格要求，所以許多的研究都朝向不受時步限制的隱式（implicit）解法發展，而 Caretto et al.

（1972）、Patankar ＆ Spalding（1972）[3]及 Patankar（1980）

[4]所提出之 SIMPLE 算則，即是基於疊代原理所提出的隱式解法。

SIMPLE 方法的基本概念在於求解強藕合之動量與連續方程式時，將連續方程式藉由增量型之動量方程式的引入，化為壓力修正方程式來配合疊代求解。其最大特點在於對增量型動量差分方程式之假設，而忽略了主格點速度修正量的影響，故其引用之動量差分方程式與原始之動量微分方程式，便具有相當程度的不一致性，導致求解壓力場時收斂狀況不佳，使收斂速度變慢。於是有人為了提升 SIMPLE 收斂速率而提出了 SIMPLER、SIMPLEC 及 PISO 等算則。

人工壓縮法（Artificial－compressbility scheme）：

在可壓縮流的領域中，發展了許多穩健的高解析算則，因此，是

(26)

否能應用這些可壓縮流的高解析算則，來求解不可壓縮流的流場問題，已為大家所注目發展之焦點。為了達此目的，Chorin（1967）[5]

提出了虛擬壓縮性法或稱為人工壓縮法，它是將壓力對時間的微分項加在連續方程式中，使得連續方程式與動量方程式能藕合在一起，成為一組雙曲線的方程組，而本文就是採用此法。

1.3 研究方法與目的

在解決流體力學的問題上，實驗方法，解析法和數值模擬是三種最主要而常見的方法。一般而言，解析法提供了迅速而完整的解﹔然而，這類的分析方法需要有較多的假設及邊界條件，使得解析法受到相當的限制。實驗方法最具有提供真實資料的優點，可惜由於花費昂貴，且實驗過程費時，量測不易，以致無法廣泛使用。近年來資訊科技突飛猛進，電腦的效能十分強大，因此利用電腦來模擬流體力學等自然現象之計算流體力學成為了在一般實驗方法下無法觀察及了解流體行為的重要輔助工具。

本文在時間積分上利用 Lusgs，空間差分上運用一階精度的迎風算則 (UPWIND) 來求解不可壓縮流之奈維爾 – 史托克方程式

(27)

（Navier-Stokes equation）。而在進行數值模擬時，計算網格的選取影響其結果甚大，應慎選格點，以便能正確估計出相關的物理現象發生的位置。

本文使用之統御方程式為雷諾平均奈維爾–史托克方程式，並配合不同方法的特徵邊界條件，模擬三維空穴流及背向階梯流流場，藉此了解不可壓縮流的特徵。

使用之計算時間與設備分述如下:

軟體設備：(一) 程式編輯：”Microsoft Fortran Power Station”

6.5 版

(二) 後處理：Amtec 公司所發展之”Tecplot” 9.0 版硬體設備：(一) Pentlum IV 2.4G 電腦一部，RAM 為 1024KB

(二) Pentlum IV 3.0G 電腦一部，RAM 為 1024KB 計算時間：以使用格點數 32*32*32，需費時約 24 小時。

1.4 章節介紹

本文第一章說明此研究的動機與目的，並將本研究中撰寫程式

(28)

時，所選用的算則方法做一文獻上的回顧。第二章中將說明研究流體的相關物理性質，說明物理性質。然後再將其性質作一正確的分析假設。第三章吾人將以第二章中所作之物理假設，建立出適用於三維不可壓縮流的統御方程式，即為奈維爾-史托克方程式。此外將所運用的數值方法，作逐一說明介紹。第四章中將空穴流、背向階梯流，作一比較分析。最後第五章說明對於此研究的未來發展方向。

(29)

第二章基本理論與統御方程式

2.1 前言

在流體力學中，流場分析是一項非常重要的工作，為探討流體的運動狀況，一般有二種研究的途徑 ﹕ 一為拉格朗奇描述法

（Lagrangian description）﹔另一為尤拉描述法（Eulerian description）。拉格朗奇描述法是以粒子的運動為概念，針對固定的流體粒子，以初始位置─時間座標系來描述粒子的運動軌跡。這種描述法的優點在於可追蹤每一流體粒子的軌跡及運動性質，而且是自始至終追隨某一固定的粒子，故質量守恆的關係必然成立。但是在流體力學的領域中，

吾人感興趣的是流體運動的狀態與現象，而非某一流體粒子的自我表現，因此除了特殊的問題外，一般均採用尤拉描述法而不採用拉格朗奇描述法。而尤拉描述法則是以場為概念，以空間而只探討在某一位置流體粒子的速度隨時間而變化的狀況，即著重在某一特定位置處的流體粒子的表現，此流體運動的探討方式，與一般試驗室中風洞及水洞的測試段（test section）量測原理相同，故較為方便。

(30)

恆、動量守恆、能量守恆─卻是針對某一固定流體粒子而言，亦即屬於拉格朗奇描述法的觀念範疇。因此，為了求解流場問題而推導描述流體運動的統御方程式時，必須結合尤拉描述法與拉格朗奇描述法的觀念，才能得到一組代表流體運動的統御方程式。

2.2 基本物理特性

在物理學的領域中，物質一般被區分為固體、液體及氣體三種狀態，此三者之熱力特性截然不同。但在力學的範疇內，卻可簡化為固體與流體二大類別。兩者之差異在於靜止情況下固體可以抵抗剪應力

（shear force），而流體卻不行。

一般而言，流體具有下列性質，即為運動性質（kinematic properties）─包括速度、角速度、渦度、加速度及應變率（strain rate）等﹔傳輸性質（transport properties）─包括黏滯度、熱傳導及質量擴散等﹔熱力性質（thermodynamic properties）─包括壓力、密度、溫度、比熱、普朗特常數、熵、焓等。基於上述之性質而建構出奈維爾─史托克程式，可完全描述流體流動的各種物理現象，

如震波、渦漩、邊界層、分離流及紊流等。吾人即由此出發，利用數

(31)

值模擬分析剪切混合層之結構，以了解其流場與速度或密度之間的相互關係。

2.3 統御方程式

三維含虛擬壓縮性之雷諾平均奈維爾－史扥克方程式之守恆律積分型式可表示如下

¹ ¹ ^F ^d ^S ⁰ V V

d V Q

t

^V ^S

其中，V 為任意控致單元之體積，S 為任意控制面之面積， S

d 的方向

為垂直向外，而

Q

則為守恆變數，

F E E

_v

i F F

_v

j ( G G

_v

) k

為流通向量，在卡氏座標系中，三維含虛擬壓縮性之雷諾平均奈維爾－

史扥克方程式之守恆律微分型式可表示如下

z Gv y

Fv x

Ev z

G y

F x

E t

Q

（2-1）

其中 Q 為守恆變數向量定義為

(32)

w v u p

Q

（2-2a）

E，F，G 分別為 x，y，z 方向之無黏性通量向量定義為

uw uv

p u

u E

2

（2-2b）

vw p v

uv v

F

₂ （2-2c）

p w

vw uw w G

2

（2-2d）

而黏性通量向量 E^v，F^v，G^v定義為

(33)

x z

x y

x

w u

v u Ev 2 u

0

（2-2e）

y z

y y x

w v

v u v

Fv 2

0

（2-2f）

z z y

z x

w v w

u Gv w

2 0

（2-2g）

2.4 座標轉換

為增加數值精確度與有效性，吾人經常使用座標轉換，將主導方程式於卡式直角座標系寫映至一般曲線座標系，此寫映的目的在於將所有的物體表面映成座標表面，而將預期會且有較大變化梯度，如邊界層及尾流等處之格點加密，而將原本不規則的物理域轉換成均一的計算域，如此有助於離散化公式的應用及簡化邊界條件之處理，而轉換後的方程式又可寫成守恆律形式。

座標轉換令

(34)

, ,

, , z z

y y

x x

（2-3）

則含虛擬壓縮性之不可壓縮奈維爾-史托克方程式，可改寫成在一般座標係

, ,

之守恆律形式如下：

ˆ ˆ ˆ ˆ ˆ ˆ ˆ 0

v v

v

F F G G

E E t

Q

（2-4）

其中守恆變數向量

Qˆ

定義為

w v u p h

Qˆ

（2-5a）

非黏性通量向量

Eˆ

，

Fˆ

，

Gˆ

定義為

p wU

p vU

p uU

U h

E

z y

ˆ

x （2-5b）

p wV

p vV

p uV

V h

F

z y

ˆ

x （2-5c）

(35)

p wW

p vW

p uW

W h

G

z y

ˆ

x （2-5d）

而黏性通量為

Eˆ

_v，

Fˆ

_v，

Gˆ

_v則定義為

z z y

y z x x z

z z y y

y x x y

z z x v

y x x x v

w w

v w

u

v w v

v u

u w u

v u

h E

2 2

2

0 ˆ

（2-6a）

z z y

y z x x z

z z y y

y x x y

z z x y

y x x x v

w w

v w

u

v w v

v u

u w u

v u

h F

2 2

2

0 ˆ

（2-6b）

z z y

y z x x z

z z y y

y x x y

z z x y

y x x x v

w w

v w

u

v w v

v u

u w u

v u

h G

2 2

2

0 ˆ

（2-6c）

U，V，W 為逆變速度分量，定義為

(36)

w v u W

w v u V

w v u U

z y x

（2-7）

而 h 為座標轉換之 Jacobian，其實際物理意義為控制單位之體積，

可以寫成

z z z

y y y

x x x

h

(2-8)

z y x z

y x z

y x

非黏性通量 Eˆ ， Fˆ ，

Gˆ

的 Jacobian 矩陣

Aˆ

，

Bˆ

，

Cˆ

定義為

Q A E ˆ

ˆ ˆ

，

Q B F

ˆ

ˆ ˆ ，

Q C G

ˆ

ˆ ˆ

(2-9)

其共用通式可以下列表示

w k w

k w k k

v k v k v

k k

u k u

k u k k

k k

k A

z y

x z

z y

x y

z y

x x

z y

x

i

0 ˆ

, (i=1,2,3) (2-10a)

(37)

其中

A ˆ A ˆ

1 ，

A ˆ B ˆ

2 ，

A ˆ C ˆ

3 (2-10b)

而

3 2

1

, ,

3 , 2 , 1 , ,

, k k i

k

w k v k u k

i z y z

i x y

i x

z y x

(2-10c)

矩陣

Aˆ

_i的四個相異時數特徵值分別為

) ,

, , ( ) , , ,

(

¹_i _i² ³_i _i⁴

c c

(2-11)

此處 c 為虛擬聲速定義為

c

² (2-12)

矩陣 Aˆ ， Bˆ ，

Cˆ

存在如下之相似轉換

A ^ˆ

_i

R

_i _i

R

_i¹ (2-13)

(38)

其中 _i為矩陣 Aˆ， Bˆ ，

Cˆ

的特徵值所構成的對角矩陣之通式，可表示如下

c A

_i

c

0 0 0

0 0

0 0 0

0

(2-14)

而

R

_i是與 Aˆ， Bˆ ，

Cˆ

之特徵值相對應的右特徵向量所構成的矩陣，可以寫成

z z

y y

x x

i

k u k w z z

k u k v y y

k u k u x x

c c

R

1 2

0 0

(2-15)

1

R

i 則為其反矩陣，可以寫成

z u

x

z y

x i

k k

k

k k

k

d y d x d

x d z d

z d y a

z a y a x

d x d y d

z d x d

y d z a

z a y a x

R c

1 1

2 2

2

1

2 2 2 3 2 1 2 3 2 2 2 1 2 3 2 2 2 1

2 1 1 1 1

1 3 1 3

1 2 1 1

1 3 1 2 1

2 1

(2-16)

(39)

其中的

w k

d v k

d u k

d

w k u k a v k w k a u k v k a

c c

z z y y

x x

z z y y

x x

z t

x

x z z

y y

x

i i

i

i i

i

3 2

1

3 2

1

2 1

3

2 2

1 3

2

1 1

, ,

, , ,

, ,

(2-17)

(40)

第三章數值方法

3.1 空間離散

在引用有限體積法來處理插分方程式時，將數值通量定義在單元表面上，而守恆量則定義在單元形心上。因此方程式（2-7）的半離散(semidiscrete)插分型式可寫成﹕

k j v i

k j i

S E E S

E V E

t Q

, 2, , 1

2, 1 ,

,

~

~ ˆ 1

_v _i _j _k _v _i _j _k

k j i

S F F S

F

V F

2^,

, 1 2,

, 1 ,

,

~

~ 1

2 , 1 2 ,

, 1 , ,

,

~

~ 1

k j v i k

j v i k

j i

S G G S

G

V G

（3-1）

式中

V

_i_,_j_,_k為第(i,j,k)個控制單元的體積，S 則為控制面之面積，

k j i

E

, , 2 1

~

，

k j i

F

,

2 , 1

~

，

2 , 1 ,

~

k j i

G

為、、ζ方向非黏性項之數值通量，

K j v i

E

, , 2

)

1

( ~

，

K j v i

F

.

2 , 1

~ )

(

，

2 , 1

)

,

( ~

k j v i

G

則為、、ζ方向黏性項之數值通量。關於非黏性通量部分，吾人以迎風算則（UP-WIND）來求得數值通量，並配合中央插分所得黏性項之數值通量。

(41)

3.2 迎風算則(up-wind)

對於雙曲型方程式，不同特徵值會使得波的傳播方向不同，原始的中央差分法即忽略了這種傳播方向的影響，從物理上考慮為了計算點 i 和 i+1 之間的單元形心 i+

2

1

上的數值通量，應算及各特徵向量分

量的方向性，對於向右傳播的分量，應該使用左邊的數值來計算，因為影響來自左邊，即應該使用向左差分法，對於向左傳播的分量，應該使用右邊的數值來計算，因為影響來自右邊，即應該使用向右差分法，這就是迎風算則的基本概念。

假設一維常係數線性雙曲線方程式如下

x 0 a u t

u

（3-2）

其中 a 為常數，此方程式之離散解須符合波傳遞方向，故一階精度之迎風算則可表示如下

)]

( ) (

[

₁ ₁

1 n

i n i n

i n

i

u a u u a u u

u

(3-3)

(42)

x

t

)

2 (

1 a a

a

(3-4)

)

2 (

1 a a

a

上式有另一變化型式表示如下

) 2

2 ( )

2 (

¹ ¹ ¹ ¹

1 n

i n i n

i n

i n i n

i n

i

u a u u a u u u

u

(3-5)

將式(3-3)寫成如此型式的目的為確保數值解的穩定性及符合波的傳遞性質，其優點為不需如同中央差分法增加人工黏滯項以保持波的傳遞性質。

3.3 時間離散

在每一網格的通量和計算出後，須利用時間積分求得下一時間物理變數值，直到得到穩態解或欲觀察的時間為止。通常時間積分分為顯式法及隱式法，本研究中採取 LUSGS﹝Yoon&jameson（1987)﹞

隱式時間積分法，首先對時間 n+1 時間步階的通量作線性化處理

(43)

1 2 1 2 1 2 1 2 1 2 1 2

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

ˆ ˆ ˆ

Q Q

C G G

Q Q

B F F

Q Q

A E E

Q Q

C G G

Q Q

B F F

Q Q

A E E

n v n v n v

n n n

（3-6）

其中 Aˆ，Bˆ，

Cˆ

分別為非黏性通量 Eˆ ，Fˆ ，

Gˆ

之 Jacobian 矩陣，

A ˆ

v，

Bˆ

_v，

Cˆ

_v為黏性通量

Eˆ

_v，

Fˆ

_v，

Gˆ

_v之 Jacobian 矩陣， Q =

Q

ⁿ ¹

Q

ⁿ

為時間步階之守恆變數的微變量。

依據非黏性通量 Jacobian 的特徵值正負號將其分解如下

A ^ˆ

_i

A ^ˆ

_i

A ^ˆ

_i

R

_i _i

R

_i ¹

R

_i _i

R

_i ¹ （3-7）

其中 _i為非負特徵值構成，而 _i為非正特徵值構成。

式（3-1）之空間離散式以隱式尤拉法(Euler implicit)時間離散可寫成

(44)

1 21 , ,

~ 1 ~

21 , ,

~

1 2, , 1

~ 1 ~

2, , 1

~

1 , 2, 1

~ 1 ~

, 2, 1

~

~ ,

ˆ, ,1 ˆ, , ,

n k j Si Gv n G

k j Si Gv G

n k j Si Fv n F

k j Si Fv F

n k j Si Ev n E

k j Si Ev t E

n k j Qi n

k j Qi k j Vi

.. （3-8）

其中上標 n 為時間指標。將式（3-6）線性化之通量與式（3-7）分解過的 Jacobain 代入式（3-8）之隱式部份，並忽略二次及高次項，則可得一尚未分解且對角線主控(diagonal dominance)的隱式分解法如下:

k j i k j

i

I Q

t V

, , ,

,

ˆ

_i _j_k

k j v i

k j k i j

v

S

i

Q A A S Q

A

_, _, ₁_, _,

2 , 1

, , 2,

1

ˆ ˆ ˆ ˆ

ˆ ˆ

+

^A ^A

_v

^S

_i _j_k

^Q

_i _j_k

^A ^A

_v

^S

_i ,_j,_k

^Q

_i,_j,_k 2

, 1 , , 1

2,

1

ˆ ˆ ˆ ˆ

ˆ ˆ

+

B B

v

S

i j k

Q

i j k

B B

v

S

i j _,k

Q

i_,j ₁_,k 2

, 1 ,

, , 2

, 1

ˆ ˆ ˆ ˆ

ˆ ˆ

+

n k j k i j v i k

j k i j

v

S

i

Q B B S Q

B

_, _, _,

2 , 1 ,

1 , ,

2

, 1

ˆ ˆ ˆ ˆ

ˆ ˆ

+ , , 1

2 , 1 , ,

2 , , 1

,

ˆ ˆ ˆ ˆ

ˆ

ˆ C

v

S

i jk

Q

i j k

C C

v

S

i jk

Q

i jk

C

+

n k j k i

j v i k

j k i

j

v

S

i

Q C C S Q

C

_, _,

2 , 1 1 ,

, 2 , , 1

,

ˆ ˆ ˆ ˆ

ˆ ˆ

=

E E

_v

S

_i _j_k

E E

_v

S

_i _,_j_,_k

2 , 1

2,

1

~ ~

~

(45)

F F

_v

S

_i _j _k

F F

_v

S

_i _j _,_k

2 , 1 2,

, 1

~

2 , 1 2 ,

, 1 ,

~

k j v i k

j

v

S

i

G G S

G

（3-9）

RHS

其中 I 為單位矩陣。在隱式法中我們加入黏性通量的 Jacobian 矩陣，

將有助於增加收斂速率，特別在高雷諾數的流場中，於固體邊界加密格點時尤然。

3.4 LU-SGS

Jiang＆Shu 所提出之 LU-SSOR 隱式法中，將非黏性通量 Jacobian 矩陣分解為非正和非負兩部分。文獻指出以下列所示之分解，有不錯的結果。

I C C

I B B

I A

A

_A_ˆ _B_ˆ

ˆ

_C_ˆ

2 ˆ 1 ˆ ,

2 ˆ 1

(3-10)

Jacobian 之頻譜半徑（spectral radius of Jacobian）

Aˆ為

C k

B k

A

k

_B _C

A

max ˆ ˆ ,

max ˆ ,

max

_ˆ _ˆ

ˆ (3-11)

(46)

其中

Aˆ 、 Bˆ 及 Cˆ

分別代表 Jacobian 矩陣 Aˆ、Bˆ 與

Cˆ

之特徵值。

k 為略大於 1 的數值，以確保 Jacobian 矩陣分解成對角線主控型式。

採用兩點單向（one-sided）差分則可將式（3-9）中現對網格界面計算的通量改為相部網格中心來來計算。且令網格變化不大，故對角線上相鄰網格界面有下列所示之關係﹕

2 , 1 2 , , 1 2 ,

, 1 2 ,

, 1 ,

2, , 1 2, , 1 2,

, 1 2,

, 1

, 2, , 1 2, , 1

2, , 1

2, 1

5 . 0

k j i k j K i

k j i k j i

k j i k j J i

k j i k j i

k j i k j I i

k j i k j i

S S

S

S S

S

S S

S

(3-12)

由定義可知

C B

A

B B C C

A

A ˆ ˆ

_ˆ

, ˆ ˆ

_ˆ

, ˆ ˆ

_ˆ (3-13)

而黏性通量之 Jacobian 則以其頻譜半徑來近似，即

V I C S

V I B S

V I

A

_v _A

S

^I _v _B ^J _v _C ^K

v v

v ˆ ˆ

ˆ

, ˆ , ˆ

ˆ

(3-14)

(47)

對於時間微分項，是由隱式尤拉法得到，為一階時間準度。

以上述之關係式代入式（3-9），可將隱式部份近似分解為 LU-SSOR 隱式算則如下﹕

n n n

RHS Q

U D D L

D

¹

ˆ

(3-15)

其中

2 , 1 1 , , ˆ , 2.

, 1 . 1 ˆ , .

2, . 1 , ˆ 1

, ˆ ,

ˆ ˆ

ˆ ,

,

2 , 1 1 , , ˆ , 2,

, 1 , 1 ˆ , ,

2, , 1 , ˆ 1

ˆ ˆ

ˆ

2 2

2

ˆ ˆ

ˆ

k j k i j C i k

j k i j B i k

j k i j A i

k j K i C J C

B I B

A A k

j i

k j k i j C i k

j k i j B i k

j k i j A i

S C

S B

S A

U

S S

S t I

D V

S C

S B

S A

L

v v

v

v v v

v v

V

(3-16)

此時若 t 趨近無限大，則使得此算則之計算時間步階與 CFL (courant friedrichs lewy)條件無關，成為無條件穩定的牛頓迭代隱式算則，即 LU-SGS 隱式解法。而三個運算子則分別如下﹕

碩 士 論 文 中 華 大 學

中 華 大 學 碩 士 論 文

題目：An LU-SGS upwind scheme for the 3D incompressible flow

系 所 別：機械與航太工程研究所碩士班 學號姓名：M09108024 歐 陽 弘 道 指導教授：蔡 永 培 博 士

中華民國 九十四 年 十一 月

中文摘要

Abstract

誌謝

目錄

中文摘要 ……… Ⅰ 英文摘要 ……… Ⅱ 誌 謝 ……… Ⅲ 目 錄 ……… Ⅳ 圖 目 錄 ……… Ⅵ 符號說明 ……… XI

第一章 緒 論 ……… 1

1.1 研究動機 ……… 1

1.2 文獻回顧 ……… 2

1.3 研究方法與目的 ……… 5

1.4 章節介紹 ……… 6

第二章 基本理論與統御方程式 ……… 8

2.1 前言 ……… 8

2.2 基本物理特性 ……… 9

2.3 統御方程式 ……… 10

2.4 座標轉換 ……… 12

第三章 數值方法 ……… 19

3.1 空間離散 ……… 19

3.2 迎風算則 ……… 20

3.3 時間離散 ……… 21

3.4 LU-SGS ……… 24

3.5 網格建立 ……… 27

3.5.1 幾何公式格點產生法 ……… 29

3.6 邊界條件 ……… 32

3.6.1 空穴流 ……… 33

3.6.2 背向階梯流 ……… 34

第四章 結果與討論 ……… 36

4.1 三維穩態流之應用 ……… 36

4.2 空穴流 ……… 36

4.3 背向階梯流 ……… 38

第五章 結論與未來展望……… 39

5.1 結論 ……… 39

5.2 未來展望 ……… 40

參考文獻 ……… 42

圖目錄（附錄）﹕

三維空穴流部份：

三維背向階梯流部份：

符號說明

E

F

G

Qˆ

Eˆ

Fˆ

Gˆ

Eˆ

Fˆ

Gˆ

ˆA Eˆ 之 Jacobian 矩陣

ˆB Fˆ 之 Jacobian 矩陣

Cˆ Gˆ

R

P

V

Rˆ

D

D

Aˆ

Aˆ

Aˆ

Aˆ

A

Bˆ

Cˆ

A

( )

G y

f

S

R

R

第一章 緒論

1.1 研究動機

1.2 文獻回顧

1.3 研究方法與目的

1.4 章節介紹

碩士論文中華大學

中華大學碩士論文

系所別：機械與航太工程研究所碩士班學號姓名：M09108024 歐陽弘道指導教授：蔡永培博士

中華民國九十四年十一月

中文摘要 ……… Ⅰ 英文摘要 ……… Ⅱ 誌謝 ……… Ⅲ 目錄 ……… Ⅳ 圖目錄 ……… Ⅵ 符號說明 ……… XI

第一章緒論 ……… 1

第二章基本理論與統御方程式 ……… 8

第三章數值方法 ……… 19

第四章結果與討論 ……… 36

第五章結論與未來展望……… 39

第一章緒論

第二章基本理論與統御方程式

¹ ¹ ^F ^d ^S ⁰ V V