中 華 大 學 碩 士 論 文
超音速混合層流場物理特性之數值模擬 A Numerical Investigation for the Physical
Properties of the Supersonic Mixing Layer
系 所 別：機械工程學系碩士班 學號姓名： M09708027 陳咸佑 指導教授： 蔡 永 培 博士
中 華 民 國 100 年 7 月
中文摘要
本論文之研究為建立一數值模擬程式，藉以分析超音速可壓縮流之剪切混合 層流場，研究此流場之捲起、成對及融混合併的現象。
剪切混合層是在分離板上藉由兩個速度不同，密度不同的流體再分離版末端 處產生捲起、成對及融混合併的現象。這種特殊的流場外型，事實上在許多的工 程應用上都會遭遇到，因此對這種流場模式要徹底研究，並運用在其它的工程科 學中。
在本研究中，吾人以可壓縮流領域中發展出的加權型基本不震盪(Weighted Essentially Non-Oscillatory) 算則直接求解 Navier-Stokes 方程式。而程式撰寫技 巧上採用具增強數值穩定性的Lower-upper Symmetric successive over relaxation ( LU-SSOR method ) 以求更快的收斂。文中包含了超音速剪切混合層內之渦漩運 動及流場內紊流特性，其結果顯示都有不錯的精確度，例如：混合層厚度、雷諾 剪應力。
關鍵詞:剪切混合層、加權型基本不震盪算則、混合層厚度、雷諾剪應力。
Abstract
The basic idea of this research is to establish a numerical simulation, which analyzes the compressibility of supersonic mixing layer, its rolling up, paring and mixing phenomenon.
The physical properties of plane, free shear layer generated by coflowing of two flows with different velocity ratio and density ratio are numerically investigated. It is important in mixing processes which are encountered in many engineering application.
For comparison purposes, we can utilize the “ Weighted Essentially Non-Oscillatory ＂ developed from compressible field to directly explain the configuration of Navier-Stokes equation. The programming scheme of which employs
“Lower-upper Symmetric successive over relaxation” to increase numerical stability.
This study also reveals the vortex growth and turbulence within supersonic mixing layer and henceforth a faster means of contraction, all of which exhibits accurate findings, ex: mixing layer thickness, Reynolds shear stress.
Keyword: Mixing layer, Weighted Essentially Non-Oscillatory, Layer thickness, Reynolds shear stress.
誌謝
天下沒有不散的延席，畢業不是一個結束，而是一個新的旅途的開始。
三年，在邁向人生的下一個階段之前，細數過去點滴，所有的辛苦皆化成歡笑的 淚水，成為永恆的美好回憶。
研究所的生涯轉瞬間消逝，感謝這段期間指導教授蔡永培博士給予我在學業 及人生觀念上的教導。並感謝蔡博章博士及鐘振忠博士在論文上的指教。
謝謝周中祺學長，在課業及生活中的協助畢生難忘。還有畢業的小老虎及小 朋友學長在忙碌中還要請他們給予指導。謝謝研究所的同學、朋友，在沉悶之及，
一起陪我度過無數的夜晚。
內心的感動實非筆墨得以形容，對於我敬愛的父母有說不出的歉疚與感激。
沒有他們的耐心支持與鼓勵就沒有今天的我。
最後，再次感謝以上所有幫助過我的人。由於你們的參與為我的人生增添 許多光彩，讓我活得更出色。謹將這分喜悅與你們一同分享。
目錄
第一章 導論...1
1-1 研究動機...1
1-2 文獻回顧...2
1-3 研究方法與目的...5
第二章 基本理論與統御方程式...7
2-1
前言...7
2-2
基本物理特性...8
2-3
統御方程式...9
2-4
座標轉換...10
第三章 數值方法...15
3-1
空間離散...15
3-2
WENO算則 ...16
3-2-1 WENO2...17
3-2-2 WENO3...19
3-3
時間離散...22
3-4
LU-SSOR ...24
3-5
網格建立...27
3-6
初始條件...28
3-7
邊界條件...30
3-7-1 上底及下底邊界條件...30
3-7-2 入口及出口邊界條件...30
第四章 結果與討論...34
4-1
網格測試...34
4-2
數值模擬結果...34
第五章 結論...41
參考文獻...42
表目錄
表4-1 混合剪切流之計算條件 ...59
圖目錄
圖4-1 次音速流體剪切混合層之瞬間等密度線圖 ...44
圖4-2 次音速流體剪切混合層之瞬間等渦漩線圖 ...44
圖4-3 次音速流體剪切混合層之瞬間等密度線圖 ...44
圖4-4 次音速流體剪切混合層之瞬間等渦漩線圖 ...44
圖4-5 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-6 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-7 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-8 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-9 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-10 超音速流體剪切混合層之瞬間等渦漩線圖 ...45
圖4-11 超音速流體剪切混合層之瞬間等渦漩線圖 ...46
圖4-12 超音速流體剪切混合層之瞬間等渦漩線圖 ...46
圖4-13 超音速流體剪切混合層之瞬間等渦漩線圖 ...46
圖4-14 超音速流體剪切混合層之瞬間等渦漩線圖 ...46
圖4-15 超音速流體剪切混合層之瞬間等渦漩線圖 ...47
圖4-16 超音速流體剪切混合層之瞬間等渦漩線圖 ...47
圖4-17 超音速流體剪切混合層之瞬間等渦漩線圖 ...47
圖4-18 超音速流體剪切混合層之瞬間等渦漩線圖 ...47
圖4-19 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4-20 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4-21 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4-22 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4-23 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4- 24 超音速流體剪切混合層之瞬間等密度線圖 ...48
圖4-25 次音速流體剪切混合層之瞬間等密度線圖 ...49
圖4-26 超音速流體剪切混合層之瞬間等密度線圖 ...49
圖4-27 超音速流體剪切混合層之瞬間等密度線圖 ...49
圖4-28 次音速流體剪切混合層之瞬間等渦漩線圖 ...49
圖4-29 超音速流體剪切混合層之瞬間等渦漩線圖 ...49
圖4-30 超音速流體剪切混合層之瞬間等渦漩線圖 ...49
圖4-31 次音速流體剪切混合層之瞬間等渦漩線圖 ...50
圖4-32 超音速流體剪切混合層之瞬間等渦漩線圖 ...50
圖4-33 超音速流體剪切混合層之瞬間等渦漩線圖 ...50
圖4-34 次音速流體剪切混合層之瞬間等密度線圖 ...50
圖4-35 超音速流體剪切混合層之瞬間等密度線圖 ...50
圖4-36 超音速流體剪切混合層之瞬間等密度線圖 ...50
圖4-37 無因次的剪切混合層成長率與對流馬赫數之關係圖 ...51
圖4-38 超音速流體對流馬赫數M =0.51 之平均流線方向速度輪廓之自我相似說明圖C ...52
圖4-39 超音速流體對流馬赫數M =0.45 之平均流線方向速度輪廓之自我相似說明圖C ...52
圖4-40 超音速流體對流馬赫數M =0.86 之平均流線方向速度輪廓之自我相似說明圖C ...53
圖4-41 超音速流體對流馬赫數M =0.86 之平均流線方向速度輪廓之自我相似說明圖C ...53
圖4-42 超音速情況下瞬時之密度與壓力圖 ...54
圖4-43 超音速流體(密度比值λ =1.55,速度比值λ =0.43,時間Time=0.005s, M =1.4)ρ u 1 ...55
圖4-44 超音速流體(密度比值λ =0.5,速度比值λ =0.53,時間Time=0.002s, M =1.6)ρ u 1 ...56
圖4-45 超音速流體(密度比值λ =0.64,速度比值λ =35,時間Time=0.0019s, M =1.81)ρ u 1 ...57
圖4-46 超音速流體(密度比值λ =1.55,速度比值λ =0.57,時間Time=0.004s, M =1.91)
ρ u
1 ...58
符號說明
X 物理域直角座標系統之X 方向 Y 物理域直角座標系統之Y 方向 ξ 計算域之座標軸
η 計算域之座標軸
Δ 微變量
ΔX X 座標空間間隔 ΔY Y 座標上之空間間隔 Δξ ξ座標上之空間間隔 Δη η座標上之空間間隔 t 時間座標
Δt 時間間隔
P 壓力
u X 方向之速度 v Y 方向之速度 ρ 密度
μ 動力黏度 υ 運動黏度
ij
黏滯剪應力張量Re 雷諾數
Q 物理域之無因次守恆變數向量 E X 方向之無因次非黏性通量向量 F Y 方向之無因次非黏性通量向量
E v
X 方向之黏性通量
F v
Y 方向之黏性通量
Qˆ
計算域之守恆變數向量Eˆ
ξ 方向之通量向量Fˆ
η 方向之通量向量Eˆ v
ξ 方向之黏性通量
Fˆ v
η 方向之黏性通量V 控制體積
ˆA
Eˆ
之Jacobian 矩陣 ˆBFˆ
之Jacobian 矩陣 λ 特徵值R 右特徵向量之矩陣
R 1
左特徵向量之矩陣rs 對應於第s 個特徵值之右特徵向量 ls 對應於第s 個特徵值之左特徵向量
A
m 擾動之強制振幅
G y
( ) 高斯分佈函數 c 聲速U
1 上板速度（較快之速度）
U
2 下板速度（較慢之速度）
U
兩自由流速差
U
c 兩自由流之平均速度f
m 頻率
熱容量比值
0 初始動量厚度 m 隨機相位
密度比值
u 速度比值 剪切混合層厚度
u
( , ) x y
無因次平均流線方向速度
相似座標U
統計時間內對時間作平均之流線方向速度V
統計時間內對時間作平均之跨流線方向速度'2
u
雷諾正向應力'2
v
雷諾正向應力u v
' ' 雷諾剪應力第一章 導論
1-1 研究動機
“渦漩”在自然界中是非常引人注意的特殊物理結構。許多科學家紛紛利用各 種不同的方法想探索其真正的性質，但是人類相較於自然之下是顯得相當微小 的，對於它的了解實在有限，特別是它所衍伸出的紊流特性。而剪切混合層是最 能反映出渦漩結構及特色的流場模式。
當兩個不同流速的流體均勻匯合時，其接觸面會發展出一個剪應變極大的區 域，流體質點在此作大量的動量交換，使得兩個均勻流在下游不遠處便發生混 合。這片區域沿著下游逐漸擴大，我們稱其為剪切混合層。
藉 由 這 兩 個 均 勻 速 度 的 平 行 流 之 間 的 紊 流 混 合 所 產 生 之 空 間 發 展，平面自由剪切層具有簡單的幾何性質。剪切混合層在許多工程應用上扮 演著重要角色。例如﹕混合耳、熱交換器及一般飛行器之燃燒室、後燃器等裝置 都涉及了剪切混合層的運用。剪切混合層控制混合的比率，而多種流的混合通常 對整個設計裝置之長度有決定性之影響，剪切混合層亦被認為係推進系統中，過 量噪音產生的重要因素。在理論上，剪切混合層的行為被認為是最基本的。在足 夠高之雷諾數（即黏滯力值與慣性力值作比較時，比值甚小）以及相當長之下游 距離的情況下，對於平面紊流混合層而言，適當的運動方程式與邊界條件，能於 數值模擬解。
目前我們對於超音速混合層流動力機構、控制混合與燃燒過程的了解，實在 非常有限。近一、二十年來的實驗研究多以平均流場的量測為主，少數研究定量 分析超音速混合層流的動力機構，以及其與次音速流的不同點。而在可壓縮剪切 層中紊流結構並不像不可壓縮流那樣的有規律。所以可壓縮剪切層的研究不但有 基礎研究的價值，也有很大的實際應用價值。
1-2 文獻回顧
紊流混合層的研究在這幾年來是相當受到歡迎的，除了他們在實際的空氣動 力學之重要性外，另外它們在基本的自由剪切層的研究上也是相當重要。在實際 應用上，混合層控制著燃燒室的混合效率，它們也反映出在推進系統上所產生推 進時的噪音問題，特別是在噴流(Jet)的狀況下[1,2]，因此剪切混合層在許多工程 應用上扮演重要角色。
剪切層流是藉由一分離板所產生出來的，主要的原因是因為由
Kelvin-Helmholtz 不穩定性所產生之以二度空間為主之巨大渦漩結構(Large eddy)所控制支配[3]，在流場中小的微擾(Perturbation)發展成非線性的波動；同 時並破碎分離且捲起(Roll Up)，將原始的剪切層渦漩轉化成一獨立”團狀”渦漩 [4]。在相當高之雷諾數之情況下，相似的渦漩融合亦已被人觀察得知了[3]。通 常 可 將 此 現 象 說 明 為 ﹕ 將 流 體 吸 入 剪 切 混 合 層 中 ； 較 大 的 應 變 率 (Strain Rate)通常發生在成對的過程，且兩流體介面會在此過程中逐漸增加，在較 大的應變率發生時所產生的小尺度微擾也就是造成渦 漩 成 對 交 互 作 用
(Vortex Pairing)的 原 因 。 而 渦 漩 成 對 混 合 也 可 視 為 在 次 調 合
(Subharmonic)不穩定下之產物。相關研究發現，當渦漩捲起(Roll Up)及渦漩交互 成對到流場下游時，流體速度增加，此效應傳回到流場上游，因而產生自我供給
（Self-Sustaining）不穩定性，此一不穩定性在剪切混合層中會自我刺激
（Self-Excitation），使得整個流場充滿不穩定性機制[5]。這種不穩定性機制是 Kelvin-Helmholtz 不穩定性，它的性質基本上是非黏性不穩定（Invisid Unstable）， 且與速度有著密切的關係。
連續的渦漩混合構成了流場下游混合層的擴散率及控制渦漩結構成對交互 作用下混合層的成長。在實驗方面證實了當馬赫數增加時，混合層的成長率將會 下降，此外，混合層的成長率也與二流體的密度比存在著關係[3]。
從一九七Ｏ年代開始的一系列實驗分析結果顯示，在超音速噴流中其尾流的
大小比一般低速流的區域小，針對理論的背景探討剪力層擴散率與噴流馬赫數的 關係，剪力層擴散率與流場馬赫數成反比。在不可壓縮的流場，利用不同的氣體 的組合，雖然混合剪力層的成長受到密度比的影響，但是此效應並沒有在可壓縮 流場觀察的明顯，因此可壓縮流剪力層成長率減少的現象，不能很單純的由密度 變化來解釋，與流體可壓縮的程度也有關。
許多由實驗所得之證據顯示，主要的流線方向渦漩（Streamwise vortices）存 在於Kelvin-Helmholtz 捲浪（Roller）之間的帶狀區域。帶狀區域係由較高應變 值與較低渦漩值所形成的帶狀結構組織，藉此帶狀區域可連接鄰近的翼展方向渦 漩。
在混合層方面的實驗，Lasheras 力用微小擾動的方法作用於混合層實驗中的 翼 展(Spanwise) 方 向 ， 同 時 並 在 分 離 板 處 改 變 此 擾 動 波 常 來 觀 察 流 線 渦 漩 (Streamwise Vortex)與翼展渦漩(Spanwise Vortex)之間的交互作用，以及在成對相 反方向流線渦漩下蕈形(Mushroom)的成長情況[6]。
Bogdanoff 引用參考馬赫數M+，這是根據兩流體速度值得二分之ㄧ極具幾何 特性之平均聲速所定義出來的。他將實驗結果與Ｍ＋做成一關係比較後發現Ｍ＋小 於0.8，若Ｍ＋大於 0.8，則他推測流場中會發生斜震波；且與流體流動方向呈現 的
夾角關係(7)。Papamoschou 和Roshko 依循著Bogdanoff 的概念並建立了一個由兩自由流 體速度差及兩流體平均聲速這兩個參數所構成的新參數稱之為:對流馬赫數Mc， 對流馬赫數可以描述因留場中可壓縮性(Compressibility)所引起的效應。同時他們 將剪切層成長率與對劉馬赫數畫成一關係圖比較後發現；當對流馬赫數大於 0.006 與小於 1 馬赫時，渦漩成長率隨著對流馬赫數增加而下降；當馬赫數大於 1 時，責成長率呈現”漸近線”式得接近低水平趨勢[8]。
Samimy 在超音速混合層實驗中發現；可壓縮效應除了與渦漩成長率有關之 外，也與雷諾應力有這密切關係。他同時也說明；當對流馬赫數增加時，流線與
跨流線流體的紊流強度會遞減，同時也會隨著流線位置增加而遞減，特別是在完 全發展流區域之後。此外，紊流強度方面與不可壓縮混合層的強度相比較下也降 低許多[9]。
Grinstein 將特徵邊界法[10]利用在邊界條件中，同時也証明了在渦漩結構 中，速度較快的流體分子比起速度較慢的流體分子在混合層中擁有較多的成份，
這個現象在其他的數值模擬中也有相同發現[9]，此外也證明了當渦漩捲起及成 對合併過程中，在流場出口處會因為流體速度加快而進一步產生壓力擾度傳遞到 上游流場且在入口處產生不穩定。
Tsai 在其研究中發現，當使用＂強制＂(Forced)頻率條件再流場中時，剪切層會 停止成長，且在這非成長(Non-Growing)區域中，渦漩成對現象無法發生。在使 用次調合頻率條件於留場中時，渦漩合併的情形會發生，混合層發生速度快，但 接著動量後度停止成長[11]。
SOH 在研究中說明了史特赫數(Strouhal Number)St 與流場的關係。當中說明 了在流場初始過程中St 值約為 0.021，在末期時約為 0.0086，這說明了低史特赫 數下，低頻率的波點波長較長，且停留在流場下游的時間也會比高頻波長。在流 場中，渦漩所在位置會存在一較小靜壓值(Static Pressure)，而整個流場內存在著 一列的局部靜壓極大與極小值，基本上這是分析渦漩結構內對樓速度之依據 [12]。
Ansari 利用 DNS（Direct Numerical Simulation）模擬較高與較低速度的初始 入口層流條件，並計算剪切混合層內的性質。其結果發現使用層流條件時，數值 模擬解的流場特性會較早在流場內發生，且使層流導入之入口速度曲線與使用紊 流導入之曲線其結果相符合[13]。
U
1-3 研究方法與目的
在一般流體力學的研究中，實驗方法，解析法和數值模擬是三種最主要而常 見的方法。一般而言，解析法提供了迅速而完整的解﹔然而，這類的分析方法需 要有較多的假設及邊界條件，使得解析法受到相當的限制。實驗方法最具有提供 真實資料的優點，可惜由於花費昂貴，且實驗過程費時，量測不易，以致無法廣 泛使用。近年來資訊科技突飛猛進，電腦的效能十分強大，因此利用電腦來模擬 流體力學等自然現象之計算流體力學成為了在一般實驗方法下無法觀察及了解 流體行為的重要輔助工具。
剪切混合層為兩流體具有相同的流動方向，流體之間通常存在著速度及密度 上的差異。當流體由分離板面分離之後，藉由兩流體之間存在的不穩定性使流 場之中產生渦漩捲起及合併等多種的物理現象。因此造成流場下游處充滿了 紊流的特性，這些紊流的特性成為一個重要的研究領域。
如圖（1-1）所示之剪切混合層流體的幾何外形，圖中所示，剪切混合層最 初係被一分離板隔開，且在分離板尾端混合著不同速度或密度之平面流體。在上 半部的流體具有較快速之流體 1而密度為ρ1，而下半部流體具有較慢速度U2而 密度為ρ2。假設上板部分的密度為固定值，改變下板部分的密度，則可得到不同 的密度比值。而分離板假設為水平的放置，位置在計算區域左方邊界的中間。分 離板可考慮成不具厚度，且分離板並不在計算模擬中。本模擬以各數據如表附錄 一所示。
本文利用LU-SSOR（Lower-upper Symmetric successive over relaxation）與 加權型基本不振盪算則（Weighted Essentially Non-Oscillatory）[14, 15]，求解可 壓縮流體在剪切混合層之奈維爾–史托克方程式（Navier-Stokes Equation）。而在 進行數值模擬中，計算網格的選取影響其結果甚大，應慎選格點，以便能正確估 計出各種相關的物理現象及其發生的位置。
對於分離板後的初始流線方向速度分佈，則採用雙曲線正切函數分佈來進行 數值模擬。利用這樣的速度分佈曲線，以物理角度而言是具有意義，同時更能模 擬真實的流場狀況[16]。
使用之統御方程式為奈維爾–史托克方程式，並配合不同方法的特徵邊界條 件，模擬二維次音速及超音速等剪切混合層流場，藉 此 了 解 渦 漩 動 力 學 及 渦 漩 的 變 遷 過 程，特 別 是 在 流 場 中 因 為 Kelvin-Helmholtz 不穩定性下，使得渦 漩在分離板尾端捲起、渦漩成對的交互作用及渦漩在流場下游處產生的擴大效 應、混合層成長率、巨大渦漩結構、雷諾應力及混合層流場中的相似，以上均 為本研究中所討論到的主題。
第二章 基本理論與統御方程式
2-1 前言
在流體力學中，流場分析是一項非常重要的工作，為探討流體的運動狀況，
一般有二種研究的途徑﹕一為拉格朗奇描述法（Lagrangian description）﹔另一為尤 拉描述法（Eulerian description）。拉格朗奇描述法是以粒子的運動為概念，針對固定的 流體粒子，以初始位置―時間（
x
0，t）座標系來描述粒子的運動軌跡。這種描 述法的優點在於可追蹤每一流體粒子的軌跡及運動性質，而且是自始至終追隨某 一固定的粒子，故質量守恆的關係必然成立。但是在流體力學的領域中，吾人感 興趣的是流體運動的狀態與現象，而非某一流體粒子的自我表現，因此除了特殊 的問題外，一般均採用尤拉描述法而不採用拉格朗奇描述法。而尤拉描述法則是 以場為概念，以空間而只探討在某一位置流體粒子的速度隨時間而變化的狀況，
即著重在某一特定位置處的流體粒子的表現，此流體運動的探討方式，與一般試 驗室中風洞及水洞的測試段（test section）量測原理相同，故較為方便，為大多 數的研究人員所採用。
尤拉描述法雖然較為方便，但是流體力學的三大定律―質量守恆、動量守 恆、能量守恆―卻是針對某一固定流體粒子而言，亦即屬於拉格朗奇描述法的 觀念範疇。因此，為了求解流場問題而推導描述流體運動的統御方程式時，必須 結合尤拉描述法與拉格朗奇描述法的觀念，才能得到一組代表流體運動的統御方 程式。
2-2 基本物理特性
在物理學的領域中，物質一般被區分為固體、液體及氣體三種狀態，此三者 之熱力特性截然不同。但在力學的範疇內，卻可簡化為固體與流體二大類別。兩 者之差異在於靜止情況下固體可以抵抗剪應力（shear force），而流體卻不行。
一般而言，流體具有四類性質，即為運動性質（kinematic properties）
―包括速度、角速度、渦度、加速度及應變率（strain rate）等﹔傳輸性質（transport properties）―包括黏滯度、熱傳導及質量擴散等﹔熱力性質（thermodynamic properties）―包括壓力、密度、溫度、比熱、普朗特常數、熵、焓等。基於上述 之性質而建構出奈維爾―史托克程式，可完全描述流體流動的各種物理現象，如 震波、渦漩、邊界層、分離流及紊流等。吾人即由此出發，利用數值模擬分析剪 切混合層之結構，以了解其流場與速度或密度之間的相互關係。
2-3 統御方程式
在理想氣體狀況之假設條件下，考慮暫態，可壓縮，具黏滯性，不考慮熱傳 導及重力因素之二維奈維爾―史托克方程式，在直角座標系下，經過無因次化可 寫成守恆律形式之偏微分方程組如下﹕
0
y F F x
E E t
Q
v v（2-1）
其中Q 為守恆變數向量定義為
e
v Q u
（2-2a）
E，F 分別為 x，y 方向之無黏性通量向量定義為
u p e
uv p u
u
E
2
（2-2b）
v p e
p v
uv v
F
2
（2-2c）
而黏性通量向量Ev，Fv定義為
x yx xx v
E M
0Re （2-2d）
y yy xy v
F M
0Re （2-2e）
以上各流場變數是以適當的自由流值做為無因次化之依據，而空間座標x，y及時 間t則分別以物體之特徵長度L及L/C∞為無因次化單位。通常無因次化密度ρ，速 度分量u，v，總能e，聲速c及黏滯係數μ分別定義為該物理量與對應之自由流值
，c
， c ，
2c
，
等 比 ，而壓力是以 c
2為無因次化單位，假設為理 想氣體則無因次化壓力可表示為﹕之 值
2 2
2
1
e
1u v
p
（2-3）其中γ=1.4 為空氣之比熱比值。在（2-2b,c,d,e）中剪應力及熱通量項可表示為﹕
y x
yx xy
y y
x yy
x y
x xx
v u
v v
u
u v
u
2
2（2-4a）
Pr 1
Pr 1
x xx xy
y y
T u v x
T u v y
x yy
（2-4b）
λ 則根據 Stokes 之假設等於
32 ，黏滯係數μ 則由 Sutherlan’s law 而得
5 .
1
1C T T
C
（2-5）其中T為無因次化溫度，而 C為Sutherlan’s常數其值為 110.4/T∞。M∞為自由流馬 赫數，雷諾數定義為
M c L
Re
，普朗特數（Prandtl number）定義t p
k
c
Pr 定壓比熱而
k 為熱傳導係數。在層流
t （laminar flow）中 P，
其中
c
p為 r 值一般取 0.72。-4 座標轉換
2
與有效性，吾人經常使用座標轉換，將主導方程式於卡式 直角座標系寫映至一般曲線座標系，此寫映的目的在於將所有的物體表面映成座 標表面，而將預期會且有較大變化梯度，如邊界層及尾流等處之格點加密，而將 原本不規則的物理域轉換成均一的計算域，如此有助於離散化公式的應用及簡化 邊界條件之處理，而轉換後的方程式又可寫成守恆律形式。
座標轉換令 為增加數值精確度
x y t , ,
x y t , ,
（2-6）則方程式（2-1）可改寫成
ˆ ˆ
ˆ ˆ
0ˆ
v vF F E
Q E
（2-7）
其中
（2-8a）
（2-8b）
（2-8c）
U，V，定義為
t
e
v J u
Q
ˆ 1
p U p e
p vU
p uU
U J
E
t y x
ˆ 1
p V p e
p vV
p uV
V J
F
t y x
ˆ 1逆變速度分量
v u V
v u
y x t
y x t
U
（2-9）
而黏性通量為
y y x x
yy y yx x
xy y xx x
v
J
E M
0ˆ Re （2-10a）
y y x x
yy y yx x
xy y xx x
v
J
F M
0ˆ Re （2-10b）
其中座標轉換之Jacobian 及各偏微分幾何項之定義如下﹕
Jx Jy Jx Jy
x y x
（2-11）
以上（2-7）式至（2-11）式即為二維非穩態可壓縮奈維爾―史托克方程式在一般
y x y
x J
1 曲線座標系上之形式。
在此先說明非黏性通量
Eˆ
、Fˆ
的Jacobian 矩陣定義為﹕Q
Eˆ A ˆ ˆ
，Q
Fˆ B ˆ ˆ
其通式可以下式表示0 0
(2 ) ( 1) ( 1)
ˆ ( 1) (2 ) ( 1)
( 2 ) ( 1) ( 1) ( 1) ( 1)
x y
x x x y
i
y y x y
t x t y t
k k
u k p k u k v k u k
A v k p k u k v k v k
e p k e p u k e p v
x
y
，
i 1 , 2
其中的
A
A ˆ
1
，A ˆ
2 B
y x
i k
k
v k u k
y i y x
i x
y x
2
1
,
2 , 1 ,
,
而Jacobian 矩陣的四個特徵值分別如下:
i1, i2, ,i3 i4
, , c
,
c
在此定義c 為聲速
A ˆ
，Bˆ
存在如下之相似轉換ˆ i R i i R i 1
A
則矩陣A ˆ
，Bˆ
的特徵值所構成的對角矩陣，可表示如下0 0 0 0
0 0 0
0 0 0
0 0 0
i
C
C
R
i是特徵值對應的右特徵向量構成的矩陣2 2
2 2
1 0
/ ( / ) ( /
/ ( / ) ( /
( / ( /( 1) ( /( 1)
1( )
2 / ) ( / / )) ( / / ))
y x x
x y y
i
y
x x x x
x
u k u k c u k c
v k v k c v k c
R
k u p c p c
u v
c k u k v c k u k v k v
) )
則為其反矩陣，為左特徵向量構成之矩陣
2 2 2
1
1 / ( 1) / ( 1) / ( 1) /
( / / ) / ( / ) / ( / ) / 0
/ ( ) ( / ( 1) ) ( / ( 1) ) / ( 1) / ( ) ( / ( 1) ) ( / ( 1) ) / ( 1)
y x y y
i
x y
x y
p c c u c v
2k u k v k k
R p c k u k v
p c k u k v
c
其中的
2 2
/ 2
/ /
x y
x y
c k k k u k
v
第三章 數值方法
3-1 空間離散
在引用有限體積法來處理插分方程式時，將數值通量定義在單元表面上，而 守恆量則定義在單元形心上。因此方程式（2-7）的半離散(semidiscrete)插分型式 可寫成﹕
t
Qˆ
1,
1,, 2 2
1
v v
i j i j
i j
E E S E E S
V
, 1
, 1, 2 2
1
v v
i j i j
i j
F F S F F S
V
（3-1）式中
V
ij：第(i.j)個控制單元的體積 S ：控制面積j i
E , 2 1
~
，2 , 1
~
j i
F
:
、
方向非黏性項之數值通量j v i
E
, 2)
1( ~
，
2 , 1
~ )
( F
v i j :
、
方向黏性項之數值通量關於非黏性通量部分，吾人以加權型基本不震盪算則（WENO）來求得數值 通量，並配合二階、四階中央插分所得黏性項之數值通量。
3-2 WENO 算則
關於WENO 算則，以式（3-1）的半離散差分式著手。首先把物理 通量向量
E ˆ ˆ Q
總體或局部地分離成正、負兩部分
E ˆ Q ˆ E ˆ
Q ˆ E ˆ
Q ˆ
（3-2a）當中的
E ˆ
Q ˆ / Q ˆ 0 E ˆ
Q ˆ / Q ˆ 0
，分離方式採用Lax-Friedrich 的通量 分離法求解流場。總體分離形式
E Q E ˆ Q ˆ Q ˆ
2 ˆ 1 ˆ
max
（3-2b）
局部分離形式
E Q E ˆ Q ˆ Q ˆ
2 ˆ 1
ˆ
（3-2c）
1 2 3 4
max
diag
max,
max,
max,
max
其中 是將所有計算單元內不同特徵值
的最大值所組合而成的對角矩陣，而
diag
1 , 2 , 3 , 4
則代表局部之特徵值的絕對值組合而成的對角矩陣。
將物理通量分離成正、負之後，以
2 1
~
i
E
及 2 1
~
i
E
分別表示 的正、負數值通量。將兩正、負數值組合即為非黏性項的數值通量。
ˆ ) ˆ Q ( E
即
2 1 2
1 2
1
~
~
~
i i
i
E E
E
（3-3）必須注意的是，當求取數值通量時，在
2
1
i
面的數值通量必須在特徵場下求得，再投射回物理空間。
定義 ， 分別為第
r
sl
s s 個特徵值所構成的右特徵向量及左特徵向量。3-2-1 WENO2
第s 個特徵場的正數值通量表示成
1 0, 0 1, 1
2 ,
( 1, 2,3, 4)
S S
i s
E q q s
q
0 、q
1為不同組合的stencils，定義如下:
l
s E
i l
s E
iq ˆ 3 ˆ
2 1
1 0
11
ˆ ˆ
2 1
i s i
s
E l E
l q
而
0 , s
， 1 , s
為第 s 個特徵場的加權值，表示如下:1 , 0 ,
, 1 , 0
,
,
k
s s
s k s
k
1,
1,
22 , 0 ,
0
3
, 2 3
1
s
IS
s
s IS
s
1,
1
22 1 ,
0
ˆ ˆ , ˆ
ˆ
s
l
s E
i l
s E
iIS
s l
s E
i l
s E
iIS
式中
於算例中選取 10
6，主要是避免分母為零的情況發生，文獻亦提出取值與 範圍內，不會影響數值解的真確性。同樣地第s 特徵場之負數值通量 列如下定義﹕
10
510
7
0, 0
1, 12,
1
q q
E
S Ss
i
q
0 、q
1為不同組合的stencils，定義如下:
10
ˆ ˆ
2 1
i s i
s
E l E
l q
1
21
3 ˆ ˆ
2 1
i s i
s
E l E
l q
而
0,s，
1,s 為第 s 個特徵場的加權值，表示如下:1 , 0 ,
, 1 , 0
,
,
k
s s
s k s
k
1,
1,
22 , 0 ,
0
3
, 1 3
2
s
IS
s
s IS
s
1,
2 1
22 1
,
0
ˆ ˆ , ˆ ˆ
s
l
s E
i l
s E
iIS
s l
s E
i l
s E
iIS
最後將正、負數值映射回物理空間，組合為非黏性之數值通量，如下表示﹕
4
1 1
1 ,
2 2
i i s s
s
E
E
r
4
1 1
1 ,
2 2
i i s s
s
E
E
r
2 1 2
1 2
1
~
~
~
i i
i
E E
E
3-2-2 WENO3
同樣地第s 個特徵場的正數值通量表示成﹕
S
S
S Ss
i
q q q
E
0, 0 1, 1 2, 2,2,
1
q
0 、q
1、q
2為不同組合的stencils，定義如下:
l
s E
i l
S E
i l
s E
iq 2 ˆ 7 11 ˆ
6 1
1 2
0
1
11
ˆ 5 2 ˆ
6 1
i s i S i
s
E l E l E
l q
1
22
2 ˆ 5 ˆ
6 1
i s i S i
s
E l E l E
l q
而
0 , s
、
1,s 、
2,s為第s 個特徵場的加權值，表示如下:, ,
0, 1, 2,
, 0,1
k s k s
s s S
k
0,
2,
0
10
1
s
IS
s
， 1,
1,
210
1
s
IS
s
， 2,
2,
210
1
s
IS
s
2 1
2
2 1
2 ,
0
ˆ ) ˆ 3
ˆ 4 4 (
1
ˆ 2 ˆ
12 ˆ 13
i S i
S i
S
i S i s i
s s
E l E l E
l
E l E l E
l
IS
2 1 1
2 1 1
, 1
ˆ ) ( ˆ
4 1
ˆ 2 ˆ
12 ˆ 13
i S i S
i S i s i
s s
E l E l
E l E l E
l IS
2 2 1
2 2 1
, 2
ˆ ) 4 ˆ
3 ˆ 4 ( 1
ˆ 2 ˆ
12 ˆ 13
i S i S i S
i S i s i s s
E l E l E l
E l E l E l IS
第s 特徵場之負數值通量如下定義﹕
S
S
S Ss
i
q q q
E
0, 0 1, 1 2, 2,2,
1
q
0 、q
1、q
2為不同組合的stencils，定義如下:
1
10
ˆ 5 2 ˆ
6 1
i s i S i
s
E l E l E
l q
1
21
2 ˆ 5 ˆ
6 1
i s i S i
s
E l E l E
l q
1
2
32
11 ˆ 7 2 ˆ
6 1
i s i
S i
s
E l E l E
l q
而
0,s、
1,s 、
2,s為第s 個特徵場的加權值，表示如下:1 , 0 ,
, 2 , 1 , 0
,
,
k
S s s
s k s
k
0,
2,
0
10
1
s
IS
s
， 1,
1,
210
1
s
IS
s
， 2,
2,
210
1
s
IS
s
2 1 1
2 1 1
, 0
ˆ ) ˆ 3
ˆ 4 4 (
1
ˆ 2 ˆ
12 ˆ 13
i S i S i
S
i S i s i
s s
E l E l E
l
E l E l E
l IS
2 2
2 2 1
, 1
ˆ ) ( ˆ
4 1
ˆ 2 ˆ
12 ˆ 13
i S i S
i S i s i s s
E l E l
E l E l E l IS
2 3 2
1
2 3 2
1 ,
2
ˆ ) 4 ˆ
3 ˆ 4 ( 1
ˆ 2 ˆ
12 ˆ 13
i S i S i
S
i S i s i
s s
E l E l E
l
E l E l E
l IS
同樣地最後將正、負數值映射回物理空間，組合為非黏性之數值通量，如下﹕
4
1 1
1 ,
2 2
i i s s
s
E
E
4
1 1
1 ,
2 2
i i s
s
E
E
s
r
r
1 1
2 2
i i
E E E
1i 2
WENO2 之精度為 3 而 WENO3 之精度為 5。
3-3 時間離散
在每一網格的通量和計算出後，須利用時間積分求得下一時間物理變數值，
直到得到穩態解或欲觀察的時間為止。通常時間積分分為顯式法及隱式法，本 研究中採取隱式解法。
首先對時間n+1 時間步階的通量作線性化處理
1 2 1 2 1 2 1 2
ˆ ˆ ˆ
ˆ ˆ
ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ
ˆ ˆ ˆ
Q Q
B F F
Q Q
A E E
Q Q
B F F
Q Q
A E E
n v n v n
v
n v n v n
v
n n n
n n n
（3-4）
Aˆ
，Bˆ
:非黏性通量Eˆ
、Fˆ
之Jacobian 矩陣Aˆ v
，Bˆ v
:黏性通量Eˆ v
，Fˆ v
之Jacobian 矩陣 Q
:時間步階之守恆變數的微變量依據非黏性通量Jacobian 的特徵值正負號將其分解如下
A ˆ
i A ˆ
i A ˆ
i R
i
iR
i1 R
i
iR
i1 （3-5）其中
