國 立 交 通 大 學
機械工程學系
碩士論文
利用 LES 法研究矩形管道流場之變化
An investigation of duck flow by Large Eddy
Simulation method
研 究 生：林煒峰
指導教授：傅武雄 教授
利用 LES 法研究矩形管道流場之變化
An investigation of duck flow by Large Eddy Simulation method 研究生：林煒峰 Student：Wei-Fong Lin 指導教授：傅武雄 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 June 2007 Hsinchu,Taiwan,Republic of china 中華民國九十六年六月
利用 LES 法研究矩形管道流場之變化
研究生：林煒峰 指導教授：傅武雄 教授
摘要
本文目的主要利用有限差分數值方法探討可壓縮流體的紊流三維渠道流場。採 用 LES 模 式 計 算 雷 諾 應 力 項 ， 非 黏 滯 性 項 則 用 Roe scheme 計 算 。 並 利 用 Preconditioning 法增加計算的收斂速度，使得不論在高速或低速流區域均可求得 正確結果。探討所得之雷諾應力分布，平均速度分布之變化，並和以往結果比較。
An investigation of duck flow by Large Eddy Simulation method Student：Wei-Fong Lin Advisor：Wu-Shung Fu
Department of Mechanical Engineering National Chiao Tung University
Abstract
The aim of this study is to investigate a turbulent channel flow of compressible fluid numerically. The numerical method is finite difference method. LES(Large Eddy Simulation) method is adopted to display turbulent stress term and an approximation Roe scheme is used to solve an inviscid term. Besides, the preconditioning method is added to improve convergence speed of iteration procedure and cause accurate results to be obtained in all speed ranges. The results of Reynolds stresses , average velocity and vorticities of the flow field are examined and validated with existing paper.
誌謝 我由衷的感謝指導老師傅武雄教授在這兩年來給予課業和論文上的指導，以及在 生活各方面上的關心與教誨，在此謹致最高的謝忱與敬意。同時也感謝機械系諸 師長在課業方面的指導。以及交大統計所的洪慧念老師，台大機械系黃美嬌老 師，中央大學地球科學系的蔡武廷老師在研究上的指導。學長李崇綱在數值模擬 方面的協助，實驗室同窗陳鴻儀、黃玠超、陳沅佑、學弟、學長們在精神上的鼓 勵，讓我得以突破難關。 更要感謝父母含辛茹苦的養育之恩，與家人和親友的支持與關懷，今日方能順利 完成學業。最後謹將此喜悅與所有關心我的人共同分享。
目錄
中文摘要………i 英文摘要………ii 致謝………iii 目錄………iv 表目錄………v 圖目錄………vi 符號表………viii 一、緒論………1 二、物理模式………4 2-1 Governing equation………4 2-2 邊界條件、初始條件………5 三、數值模式………12 3-1 MUSCL 法……… 12 3-2 Roe scheme………15 3-3 LES model………22 3-4 Preconditioning……… 25 3-5 LUSGS scheme………29 四、結果與討論………34 4-1 初始條件確認……… 34 4-2 Channel flow………45 五、結論………61 參考文獻………62 附錄………64表目錄
表 3-1 精度係數值………14 表 4-1 流場物理條件………47
圖目錄
圖2-1 channel flow幾何形狀示意圖………11 圖 3-1 黎曼問題特徵值結構圖………21 圖 4-1 初始條件確認，平均速度U ………36 圖 4-2 初始條件確認，平均速度V ………37 圖 4-3 初始條件確認，平均速度W ………38 圖 4-4 初始條件確認， Rxx uτ ………39 圖 4-5 初始條件確認， Ryy uτ ………40 圖 4-6 初始條件確認， Rzz uτ ………41 圖 4-7 初始條件確認，Rxy2 uτ ……… 42 圖 4-8 初始條件確認，R ……… 43 xz 圖 4-9 初始條件確認，R ……… 44 yz 圖 4-10 channel flow 無因次平均速度與理論值對照………48 圖 4-11 channel flow， Rxx uτ ………49 圖 4-12 channel flow， Ryy uτ ………50 圖 4-13 channel flow， Rzz uτ ………51 圖 4-14 channel flow， R2xy uτ − ………52 圖 4-15 channel flow， xy xx yy R R R − ………53圖 4-16 channel flow， Wx uτ ………54 圖 4-17 channel flow， Wy uτ ………55 圖 4-18 channel flow， Wz uτ ………56 圖 4-19 Rxx uτ ，初始條件亂數……… 57 圖 4-20Rxy2 uτ ，初始條件亂數………58 圖 4-21 channel flow 不同ε比較， Rxx uτ ………59 圖 4-22 channel flow 不同ε比較， R2xy uτ − ………60
Nomenclature ρ density[ 3 kg m⋅ − ] P pressure[ 2 N m⋅ − ] ij τ stress[ 2 N m⋅ − ] p
C constant-pressure specific heat[ 1 1
J kg⋅ − ⋅k− ]
v
C constant-volume specific heat[ 1 1
J kg⋅ − ⋅k− ] e internal energy[ 1 J kg⋅ − ] T Kelvin temperature[ K ] k thermal diffusivity[ 1 1 W m⋅ − ⋅k− ]
u velocity component in x -direction[ 1
m s⋅ − ]
v velocity component in y -direction[ 1
m s⋅ − ]
w velocity component in z -direction[ 1
m s⋅ − ] β pressure gradient[ 3 N m⋅ − ] Ac area[m ] 2 t Δ time difference[ s ] τ
Δ large eddy turn over time[ s ]
ij R Reynolds stress[ 2 N m⋅ − ] a sound speed[ 1 m s⋅ − ] R gas constant[ 1 1 J kg⋅ − ⋅k− ] 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[ 1
m s⋅ − ]
第一章 緒論
在自然界中最常發生的流場行為是紊流，在工程應用上最常遇到的流體問題也 多為紊流。大氣邊界層中的流場、河川中的水流、燃燒室裡的氣流、流動在飛機機 翼上邊界層中的氣流等等，紊流無所不在。因此為了改善或是解決紊流所以發生的 種種問題，必須努力研究它。然而研究紊流雖然可依靠實驗的觀測，但是有時因為 環境以及外在的條件因素下的限制，無法有效的依賴實驗的數據。此時便必須仰賴 數值計算的模擬。紊 流 數 值 計 算 大 致 上 分 為 三 類 ， DNS(Direct Numerical Simulation) 、 RANS(Reynolds Average Navier-Storkes Equation) 以 及 LES(Large Eddy Simulation)三類。DNS 主要是直接求解 Navier-Stokes Equation，黏滯係數則使 用動黏滯係數而不使用任何的模式化，直接計算流體在流場中的動態變化。因此其 網格的尺寸須小於最小渦流的尺寸，以求得較佳的計算解。故在計算高雷諾數時需 要大量的網格數涵蓋其流場，以使精確度提高，因此對電腦資源和計算時間上的需 求過大，使得運用性降低。RANS 則是以特定的紊流模式去產生流場中的等效黏滯 係數，利用此等效黏滯係數計算出流場的平均解，亦即產生的解為時均化的結果。 雖然 RANS 較 DNS 在時間上以及電腦資源上的需求為低，但是其仍存在運用上的缺 點。第一、就其時均化的結果來說，工程上的應用有時需要非時均化的結果來表現 紊流的特性。舉例來說，在噪音工程中，需要壓力的變動來計算聲音的強度，因此 需要的是非時均化的結果。RANS 便無法滿足其需求。第二、RANS 中的紊流模式化係 數需要從實驗中取得，與其流場的幾何尺寸有關。因此需要依靠實驗的量測數據方 可得到所需的係數，但有時計算的模型不易實行實驗，甚至無法實驗的，因此應用 上仍有其缺陷。LES 將流場分為大尺度(Large scale)以及次格點尺度(Subgrid scale)。對於大尺度的流場直接用 Navier-Stokes Equation 求解，而次格點尺度則 需要模式化。LES 所計算出來的結果為暫態的結果，而且所計算的雷諾數可比 DNS 還高。網格數也比 DNS 較少，因此 LES 有著計算時間較 DNS 迅速。且 LES 可計算暫 態結果，因此現實面上的應用比 RANS 較為廣。所以 LES 在現今運用上十分的廣泛。 然而使用 LES 計算紊流時，進口和初始條件對於整個 LES 模式計算會有很大的影響。
Najjar et al.[1]在其論文中利用LES模式計算模擬不可壓縮的channel flow流場。 其邊界條件在流動方向為循環性邊界，壁面為不可滑移條件(no-slip condition)。 其中LES模式中的SGS(subgrid scale) model，採用的是dynamic model。Dynamic model可隨著流場的變化調整其SGS model中的係數，使得程式可計算流場中energy back scatter的情形。計算結果與J.Kim and P.Moin[2]以DNS模擬的流場比對，不 論是在平均速度或是Renolds stress都非常的相近。顯示出利用LES法亦可模擬出準 確的數據。在[1]中計算模擬的是不可壓縮流場， Alkishriwi[3]在其論文中利用LES 並加上Preconditioning來計算模擬可壓縮channel flow。。其計算結果與Moser RD[4] 所使用DNS計算結果比對，在平均速度以及Renolds stress十分的相近。顯示出LES 與Preconditioning結合計算流場的準確性。在[1]與[3]中模擬的是具有壁面邊界的 流場。因為進出口的幾何形狀相同，所以[1]與[3]皆是使用循環性邊界，壁面也是 不可滑移條件。Mare[5]在其論文中利用LES模擬不可壓縮的plane jet流場。因為其 進出口幾何形狀並不相同，所以無法使用循環性邊界。在[5]中利用在時間域具有正 交性質的隨機亂數，將隨機亂數以線性組合的方式模擬出進口條件。然而其線性組 合的係數則必須求解homogeneous turbulence的two point correlation而得到。或 是利用已知的correlation求得。其出口條件在速度以及壓力使用黎曼邊界條件。模 擬出的進口條件經過平均後與所預定的進口條件相比十分的接近，可知利用此法可 以準確的產生所需的進口條件。 最後之結果在Reynolds stress與平均速度和 P.Spalart[6]中的DNS結果比對。平均速度在進口或是出口都很準確。但是Reynolds stress則在進口處附近有些許的誤差，不過整體來看仍是與[6]中的數據吻合。 綜合以上的論文，可知利用LES模式或是以其他方法計算紊流時，使用的邊界條 件有以下兩種： 1.為週期性邊界。即為將所算出的出口速度改變為進口的條件。此種方法十分的簡 便，但卻有其缺點之處。首先週期性邊界容易使流場產生衰減，使紊流衰退成層流。 因此比需外加source term來驅動流體。其次容易受幾何形狀的限制，出口和進口的 幾何形狀須一致，才能夠使用此種方法，如計算方形管道的流場、channel flow等。 而在Lund[7]中發展出另一種方法來改善此缺點。將出口的速度利用座標的轉換，因
入其入口條件。此方法為 semi-periodic而不是完全的週期性條件，且此法無法使 用在計算噴流等流場上。 2. 利 用 其 他 的數 值 方 法， 模 擬 出 人 工 的 紊 流。 如 在Lee[8] 中 則 先求 得 所 需 的 homogeneous turbulence的速度曲線在波數空間上的係數，在將此利用Fourier inverse transform轉換回直角座標空間，以求得所需的速度曲線。然而此種方法需 要進行Fourier transform的過程，增加了計算以及程式寫作的難度。在Klein et al.[9]中假設紊流的擾動為一連串隨機數的線性組合，而其中的組合係數為代求定 的解。利用homogeneous turbulence的two point correlation來求解其中的係數。 係數一經求出後即可套用於所有的時階，因此方便於程式的寫作。除此法以外，在 Noda[10]中利用正弦、餘弦以及homogeneous turbulence 的energy spectrum結合 以算出所需的速度擾動。但是[9]、[10]其所得到的紊流進口條件為homogeneous的 紊流。假如需要特定的紊流形式，則使用[9]、[10]中的方法將無法求得。如欲得到 特定的紊流形式，則需利用[5]中所使用的方法。為直接在一預定的平均速度曲線上 加上隨機的擾動，而將此隨機的擾動符合所需的Reynolds stress tensor或是其他 的物理條件以作為所需的暫態紊流進口條件。在[8]中便以此種方法作為計算 compressible isotropic turbulence的進口條件。而在Lee and Moin[11]、[12] 中則利用此種方法去計算背階(Backward-facing step)流場。此種方法不僅不受進 口的幾何形狀的限制，而且產生快速。最重要的是所產生的紊流可符合所需的條件。 李[13]已成功的將LES與Precondition法結合，但其美中不足的地方為進口條 件。在[13]中進口條件為一平均曲線加上隨機亂數擾動，此亂數擾動並不符合任何 的物理意義。如Reynolds stress或是homogeneous two point correlation等。因 此本研究將進行對模擬的紊流流場作驗證。利用隨機亂數和Reynolds stress tensor 來產生流場的初始速度，用以模擬channel flow流場。並利用Reynolds stress以及 平均速度來測試初始條件的正確性。最後計算其模擬結果之Reynolds stress、平均 速度、vorticity以驗證結果之正確性。
第二章、物理模式
本研究主要模擬的是 channel flow，其流場的幾何形狀如圖 2-1。其幾何外觀 如圖 2-1 所示。x 方向為 stream wise 的方向、y 為垂直壁面方向最後 z 為 span wise 的方向。 u 、 v 、 w 分別為其所對應的速度。在本章的其他小節中則介紹包括 governing equation、邊界條件、初始條件以及流場的一些基本假設。在第一節介 紹計算時所用到的方程式，以及流場基本假設。第二節為邊界條件、初始條件的取 法。 2-1、Governing equation： 本程式主要是求解 Navier-Stokes 方程式，主要內容如下： 3 1 2 F 0 F F U t x y z ∂ ∂ ∂ ∂ + + + = ∂ ∂ ∂ ∂ (2-1) 其中， 1 2 3 ( , , , , )T U = ρ ρ ρu u ρu ρe 。 1 1 1 2 2 2 3 3 3 ( ) i i i i i i i i i i i i ij j i u u u P F u u P u u P T e P u u x ρ ρ δ τ ρ δ τ ρ δ τ ρ τ λ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − ⎟ ⎜ ⎟ =⎜ + − ⎟ ⎜ + − ⎟ ⎜ ⎟ ∂ ⎜ + − − ⎟ ⎜ ∂ ⎟ ⎝ ⎠ 1, 2,3 i ∀ = 。 (2-2) ρ為密度， P 為壓力。u 、1 u 、2 u 分別為 x 、 y 、 z 方向的速度。3 λ ρ= C kp ， k 為 thermal diffusivity。 2 2 2 1 2 3 1 ( ) 2 v e=C T+ u +u +u ，C 為 等 容 比 熱 。v τij為 stress tensor。 在計算時所採取的基本假設為以下： 1. 假設工作流體為空氣，理想氣體。物理性質皆取在 298k 時的值。 2. 假設流體為牛頓流體，黏滯性為等方向性。
4. 所有壁面皆為不可滑移條件〈no-slip condition〉。 2-2、邊界條件、初始條件： A、入口條件： 在紊流計算中，入口條件的取法是很重要的關鍵。在本程式中使用兩種取法， 一為循環性邊界〈periodic condition〉，另一個則使用[5]在其論文中所提到的方 法。其以一特定的平均速度曲線加上線性組合之擾動隨機擾動，此線性組合係數須 符合特定的物理性質，如雷諾數等。將為此作詳細的解說。 循環性邊界： 循環性邊界為一種較為簡單的邊界條件設法。其優點為快速不耗費計算時間， 缺點為當計算的流場若進出口幾何形狀不相同時便無法使用。所以較常被用來計算 平板流或是 square channel flow 等完全發展的流場。其設置的方法如下：
(0, , ) ( , , ) P j k =P nx j k (0, , ) ( , , ) u j k =u nx j k (2-3) (0, , ) ( , , ) v j k =v nx j k (0, , ) ( , , ) w j k =w nx j k 其中 0 代表進口的邊界虛格點〈ghost cell〉， nx 則為出口的邊界格點。 P 為壓力， u 、 v 、 w 分別為 x 、y 、 z 、方向的速度。 然而在計算完全發展紊流時，需要以壓力梯度來當作 driving fore。因此需要一線 性壓力梯度，使流場具有動力來源而不會衰減。加入的方法如下： 將壓力分解為兩部份，一為梯度項，一為循環項： ( , , ) p( , , ) P i j k =βx+P i j k (2-4) 其中β為線性壓力梯度，P 為循環性壓力。在計算時將 xp β 梯度項分離出來，將它 視為 source term。因此(2-3)式中壓力項改變為： (0, , ) ( , , ) p p P j k =P nx j k (2-5) 但是在程式計算中，流量會因為壁面摩擦或是數值本身計算上的消散而使得流量不 守恆。因此必須在每個時階下調整β值，利用β值來固定流量。其方程式如下：
1 1 0 1 [( ) 2( ) ( ) ] n n m m n m n t Ac Ac Ac β + =β − − + − Δ (2-6) 其中 Ac 為流體流動之截面， m 為流量。 tΔ 為時階大小。 0 () 代表為初始狀態所算出 的值，也就是所希望流場固定的流量。()n 為第 n 時階所計算出的值，同理 1 ()n+ 為n+1 所計算出的值。 隨機變動邊界： 當進出口幾何形狀不同時，或是需計算非完全發展的流場時。就必須使用會隨 時間變化的進口條件。 首先對應網格 i 需產生三組高斯隨機變數a 、1 a 、2 a 。假設需要 N 組的進口條件，3 則a 、1 a 、2 a 皆需有 N 個高斯隨機變數。 3 2 2 n a =a ∀ = ∼n 1 N (2-7) 1 1 n a =a ∀ = ∼n 1 N (2-8) 3 3 n a =a ∀ = ∼n 1 N (2-9) 令所產生的高斯隨機亂數之平均值為 0 ，亦即： 1 1 1 0 N n n a a = =
∑
= (2-10) 2 2 1 0 N n n a a = =∑
= (2-11) 3 3 1 0 N n n a a = =∑
= (2-12) 接下來再使用 G.S 正交化使三組亂數互相單位化正交。因此此三組亂數便具有以下 特性。 單位化： 2 2 1 1 1 ( ) 1 N n n a a = =∑
= (2-13) 2 2 2 2 1 ( ) 1 N n n a a = =∑
= (2-14)2 2 3 3 1 ( ) 1 N n n a a = =
∑
= (2-15) 正交性： 1 2 1 2 1 0 N n n n a a a a = =∑
= (2-16) 1 3 1 3 1 0 N n n n a a a a = =∑
= (2-17) 2 3 2 3 1 0 N n n n a a a a = =∑
= (2-18) 將(2-13)~(2-18)整理如下： 0 m a = (2-19) m n a a = 0 1m n
m n
∀ ≠
∀ =
(2-20) 接下來定義 x 、 y 、 z 方向的速度為： 1 1 11 1 12 2 13 3 u = +u a u +a u +a u (2-21) 2 2 21 1 22 2 23 3 u =u +a u +a u +a u (2-22) 3 3 31 1 32 2 33 3 u =u +a u +a u +a u (2-23) 其中u 為 x 、i y 、 z 方向的平均速度，此為給定的曲線。u 為符合特定條件之隨機i 擾動。a 為隨機擾動線性組合之係數。 ij 因此必須解出a 以求得進口速度。在此利用 Reynolds stress 求解出係數。由ij (2-21)、(2-22)、(2-23)可知 x 、 y 、 z 方向的速度擾動〈velocity fluctuation〉 為： 1 11 1 12 2 13 3 u′ =a u +a u +a u (2-24) 2 21 1 22 2 23 3 u′ =a u +a u +a u (2-25) 3 31 1 32 2 33 3 u′ =a u +a u +a u (2-26) 由 Reynold stress 的定義可知： 2 1 xx R = u′ (2-27)1 2 xy R = u u′ ′ (2-28) 1 3 xz R = u u′ ′ (2-29) 2 2 yy R = u′ (2-30) 2 3 yz R = u u′ ′ (2-31) 2 3 zz R = u′ (2-32) 由(2-24)和(2-27)式可得： 2 11 1 12 2 13 3 ( ) xx R = a u +a u +a u (2-33) 若令上式之a12 =a13 = 則可得： 0 2 11 1 ( ) xx R = a u =a112 u12 (2-34) 由(2-19)式可知， 2 1 1 u = 。因此上式變為： 2 11 xx R =a (2-35) 因此可解出： 1/ 2 11 ( xx) a = R (2-36) 同理，利用(2-24)、(2-25)式代入(2-28)式可得： 11 1 12 2 13 3 21 1 22 2 23 3 ( )( ) xy R = a u +a u +a u a u +a u +a u (2-37) 因為a12 =a13= ，再令0 a23= 。則(2-37)式可變為： 0 2 11 21 1 11 22 1 2 xy R = a a u +a a u u 2 11 21 1 11 22 1 2 a a u a a u u = + (2-38) 由(2-19)及(2-20)式可知，(2-38)式中的 2 1 1 u = 且 u u1 2 = 。故(2-38)式可變為： 0 11 21 xy R a a ⇒ = 21 xy/ 11 a R a ⇒ = (2-39) 因此可由(2-39)式，利用R 及xy a 計算出11 a 。 21 利用同樣的方法，將(2-24)~(2-26)式代入(2-29)~(2-32)可得：
2 2 2 1/ 2 22 21 22 ( 21) yy yy R =a +a ⇒a = R −a (2-40) 31 11 31 11 xz xz R R a a a a = ⇒ = (2-41) 32 22 21 31 32 ( 21 31) / 22 yz yz R =a a +a a ⇒a = R −a a a (2-42) 2 2 2 2 2 2 1/ 2 31 32 33 33 ( 32 33) zz zz R =a +a +a ⇒a = R −a −a (2-43) 因此藉由(2-36)、(2-39)以及(2-40)~(2-43)式，將已知的 Reynold stress 代入就 可一步一步將係數解出。係數解出後代入(2-21)~(2-23)式即可得到 x 、 y 、 z 方向 的進口速度。 再將(2-36)、(2-39)及(2-40)~(2-43)整理如下： 11 12 13 21 22 23 31 32 33 ij a a a a a a a a a a ⎛ ⎞ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 1/ 2 2 1/ 2 11 21 2 2 1/ 2 11 21 31 22 31 32 0 0 / ( ) 0 / ( ) / ( ) xx xy yy xz yz zz R R a R a R a R a a a R a a ⎛ ⎞ ⎜ ⎟ =⎜ − ⎟ ⎜ − − − ⎟ ⎝ ⎠ (2-32) B、出口條件： 出口條件的取法可分為兩種情況。一為使用循環性邊界時，二為使用隨機變動 邊界時。 循環性邊界： 使用循環性邊界時，出口條件為： ( , , ) (1, , ) p p P nxx j k =P j k ( , , ) (1, , ) u nxx j k =u j k (2-33) ( , , ) (1, , ) v nxx j k =v j k ( , , ) (1, , ) w nxx j k =w j k 其中 nxx 代表出口的邊界虛格點〈ghost cell〉，1則為出口的邊界格點。P 為循環p 壓力項， u 、 v 、 w 分別為 x 、 y 、 z 、方向的速度。
隨機變動邊界： 超音速流：並不需給予出口條件。 次音速流：壓力為大氣壓力，速度與密度則由計算 Domain 外插得到。 C、壁面邊界： 邊界速度：不可滑移條件，u1=u2=u3 = 0 邊界壓力：Neumann 條件（利用 Poisson 方程式來求得邊界壓力）。 根據高斯散度定理，可將 Poisson 方程式寫成如下形式 Poisson 方程式： 2 2 A C p p dy dl y n ∂ = ∂ ∂ ∂
∫∫
∫
(2-34) 所以 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 + − Δ 因此在邊界的壓力可以設定成py+1= py D、初始條件： 初始條件可以採用產生隨機進口條件的方法。產生與網格數一樣多的初始條件。可 以節省計算的時間。圖2-1 channel flow幾何形狀示意圖。 X
Z Y
第三章、數值模式
本章主要說明本論文數值計算所使用之所有模式。首先第一節介紹的是 MUSCL 法， 主要是處理網格之間的物理量。第二節為黎曼解〈Riemann solver〉中的 ROE 法〈Roe scheme〉，利用 ROE 法計算出非黏滯的通量。第三節為 LES model，利用 LES Model 中的 Smargorinsky model 求出次格點尺度〈subgrid scale〉的黏滯性項。第四節 為 Preconditioning 法，當計算低速可壓縮流時，因速度和音速的數量及上差距過 大，所以為彌補此一缺點須使用 Preconditioning 法。最後第五節為 LUSGS Scheme， 在使用 Preconditioning 時已破壞了整個統御方程式。因此需使用 dual time stepping 疊代使其在 artificial domain 收斂時才能進入下一個真實時階，將在此 小節做詳細的解說。綜合上述，本論文在數值上的計算過程為，利用 MUSCL 法算出 網 格 間 的 物 理 量 ， 再 利 用 ROE 法 求 出 非 黏 滯 的 通 量 ， 在 計 算 通 量 時 加 入 Preconditioning 法。接下來使用 LES 模式求出黏滯項，其與 ROE 法求出的通量結 合得到真正的物理通量。最後使用 LUSGS 疊代以求出下一時階的物理量。
3-1、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-1) 1/ 2 1/ 2 1/ 2 R R i i i u+ = −u Δu+ (3-2) 1 / 2 (1 )( 1 ) ( 1) L i i i i i u+ β u+ u β u u− Δ = − − + − +θc(−ui−1+3ui−3ui+1+ui+2) +θd(−ui−2+3ui−1−3ui+ui+1) (3-3) 1/ 2 (1 )( 1 ) ( 2 1) R i i i i i u+ β u+ u β u+ u+ Δ = − − + − +θc(−ui−1+3ui−3ui+1+ui+2) +θd(− +ui 3ui+1−3ui+2+ui+3) (3-4) 其中(3-3)、(3-4)式中的β、θc、θd 值可由表(3-1)中查得。代入不同的值可以得
表 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 3 4 4 5 6
3-2、Roe scheme： 在此節首先介紹一維線性黎曼方程式，並導出其 exact solution。其次將其推 廣至非線性的黎曼方程式，同時也利用 Roe 法求出其近似解。 首先考慮一維線性黎曼方程式如下： 0 U U A t x ∂ + ∂ = ∂ ∂ ，其中 A 為一常數 Jacobian 矩陣。 (3-5) 初始條件為 (0) (0) (0) 1 ( , , m )T U = u u 。左上方括號代表 t=0。 求出 A 之特徵值矩陣以及特徵向量。 1 A= ΛK K− ，其中Λ 為特徵值矩陣： 1
0
0
m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 。 (1) ( ) , , m T K= ⎣⎡K K ⎤⎦ 為特徵向量，故 ( )i ( )i i AK =λK 。 接著定義特徵變數W 〈characteristic variables〉，其定義如下： ( 1) W =K − U或U =KW。因此Ut =KWt且Ux =KWx，將此結果代入(3-5)式中可得： 0 t x KW +AKW = ，可再繼續簡化成： 0 t x W + ΛW = (3-6)方程式(3-6)稱為 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-7) (3-7)可由特徵曲線法求得其解為： (0) ( , ) ( ) i i i w x t =w x−λt = i i α β 0 0 i i x t x t λ λ − < − > (3-8) 其中，αi與βi為初始值的特徵變數。由於U =KW ，可以得到 (0) ( ) ( , ) ( ) m 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-9) 除此之外，還可決定出 ( , )u x t 中的 jump uΔ ： ( ) 1 m i R L i i U U U αK = Δ = − =∑
(3-10) 其中αi =β αi− 。 i 在一維線性黎曼問題中，雖然有 exact solution，但在非線性問題裡需利用疊 代等方法，這些動作將耗費大量的時間，因此在實際應用上並不廣泛。為了解決此 問題，一般皆求解近似黎曼問題〈approximation Riemann problem〉解而不直接求 其 exact solution。在求解近似黎曼問題中最被廣泛應用的方法為 Roe[15]所提出， 亦即為 Roe scheme，其內容如下： 假設一維尤拉方程式： 0 U F t x ∂ +∂ = ∂ ∂ (3-11) 根據 chain rule，可將方程式(3-11)改寫如下： 0 U F U t U x ∂ +∂ ∂ = ∂ ∂ ∂ 再令 ( ) F A U U ∂ = ∂ ，於是方程式(3-11)可以表示成： ( ) 0 U U A U t x ∂ ∂ + = ∂ ∂ (3-12) 其中， ( )A U 就稱為 joacobian 矩陣。而 Roe scheme 將原本的 Jacobian 矩陣 ( )A U 用一常數 Jacobian 矩陣〈constant Jacobian matrix〉 (A U UL, R)代替，因此本來的黎曼問題可以改寫成近似黎曼問題： ( ) 0 U U A U t x ∂ + ∂ = ∂ ∂ (3-13) ( , 0) U x = L R U U 0 0 x x < > 於是前述方法可以得到(3-11)的近似解。由以上的原理可得知，在近似黎曼問 題上，Roe 利用常數 Jacobian 矩陣取代原本的 Jacobian 矩陣使方程式由非線性轉
變成線性，但是初始條件並沒有改變，因此可以得到方程式(3-11)的近似解。為了 要求得合理的常數 Jacobian 矩陣，須合乎 Roe 所提出的四項條件： 1. U 與 F 之間，存在著線性轉換的關係。 2. 當UR −UL→ ，則 ( ,U A U UL R)→A U( )，此處 F A U ∂ = ∂ 。 3. (A UL−UR)=FL−FR。 4. 矩陣 A 的特徵向量必須線性獨立。 這四項條件都是雙曲線方程式所需具備的，這同時也說明了 Roe 所提出的常數 Jacobian 矩陣必須有實數特徵值，其所對應的特徵向量必須線性獨立。除此之外， 條件 3.則是為了符合守恆定律（conservation law）與 Rankine-Hugoniot 條件。
線性黎曼問題的解析解，可以直接從(3-8)與(3-9)式得到， 1 2 ( / ) i U x t + 的解可以 利用下面的方程式計算： ( ) 1 0 2 ( / ) i i L i i U x t U K λ α + = +
∑
< (3-14) 或 ( ) 1 0 2 ( / ) i i R i i U x t U K λ α + = −∑
> (3-15) 其中 1 2 i+ 表示網格與網格之間的交界面(face)。 而黎曼問題的近似解，則須從解近似黎曼問題著手： ( ) 0 U F U t x ∂ +∂ = ∂ ∂ ，根據(3-13)式可得知 F =AU 為了符合守衡的條件，因此下式必須成立： ( R) ( L) ( R) ( L) F U −F U =F U −F U (3-16) 接著在固定體積的條件下，積分近似解 1 2 (0) i U + ，可得到通量(flux)的數值公式： 1 1 2 2 ( (0)) ( R) ( R) i i F F U F U F U + = + − − (3-17) 再從 F =AU的關係中可進一步求得： 1 1(0) ( R) R F AU F U AU + = + − − (3-18)再根據(3-14)式與(3-15)式可以推導出： ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i R i i i i F F U A K F U K λ α λ α+ + = −
∑
> = −∑
= (3-19) 或 ( ) ( ) 1 0 1 2 ( ) ( ) i m i i R i L i i i i F F U A K F U K λ α λ α− + = +∑
> = +∑
= (3-20) (3-19)與(3-20)所指的λi−與λi+分別是代表負的特徵值與正的特徵值，接著再利用 平均的方法將 1 2 i F + 更進一步表示成： ( ) 1 1 2 1 ( ) ( ) 2 m i R L i i i i F F U F U λ αK + = ⎡ ⎤ = ⎢ + − ⎥ ⎣∑
⎦ (3-21) 再由(3-9)式可再次改變 1 2 i F + 的形式如下： 1 2 1 ( ) ( ) 2 R L i F F U F U A U + = ⎡⎣ + − Δ ⎤⎦ (3-22) 其中Δ =U UR−UL、 1 A =A+−A− =K Λ K− ， 10
0
m λ λ ⎛ ⎞ ⎜ ⎟ Λ = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 。 當進行模擬計算時可在(3-22)式中加入一調節係數ε，使(3-22)變成以下形式： 1 2 1 ( ) ( ) 2 R L i F F U F U ε A U + = ⎡⎣ + − Δ ⎤⎦ (3-23) 其中 0≤ ≤ 。利用ε 1 ε可改變 artificial viscous term A UΔ 在計算時所佔的比 例。 接下來需找出 A 中所需的物理量，必須利用下列方法： 現考慮一維等溫尤拉方程式： ( ) 0 t x U +F U = (3-24) 其中 1 2 u U u u ρ ρ ⎡ ⎤ ⎡ ⎤ =⎢ ⎥ ⎢= ⎥ ⎣ ⎦ ⎣ ⎦ ； 1 2 2 2 f u F f u a p ρ ρ ⎡ ⎤ ⎡ ⎤ =⎢ ⎥ ⎢= ⎥ + ⎣ ⎦ ⎣ ⎦ ， a 為聲速 方程式(3-24)的 Jacobian 矩陣與其對應的特徵值與特徵向量如下所示：2 2 0 1 ( ) 2 F A U a u u U ⎡ ⎤ ∂ = = ⎢ ⎥ − ∂ ⎣ ⎦ (3-25) 特徵值：λ1= − ，u a λ2 = + u a 特徵向量： (1) 1 K u a ⎡ ⎤ = ⎢ − ⎥ ⎣ ⎦， (2) 1 K u a ⎡ ⎤ = ⎢ + ⎥ ⎣ ⎦ 接著選定 parameter vector 1 2 q U Q q u ρ ρ ρ ⎡ ⎤ ⎡ ⎤ =⎢ ⎥= ⎢ ⎥ ⎣ ⎦ ⎢⎣ ⎥⎦ (3-26) 再將 F 與U 利用 Q 表示： 2 1 1 1 2 1 2 u q U q Q u q q ⎡ ⎤ ⎡ ⎤ =⎢ ⎥= = ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ (3-27) 1 1 2 2 2 2 2 2 1 f q q F f q a q ⎡ ⎤ ⎡ ⎤ =⎢ ⎥ ⎢= ⎥ + ⎣ ⎦ ⎣ ⎦ (3-28) 為了表示出 UΔ 與 FΔ 需在定義 averaged vector Q ： 1 2 1 1 ( ) 2 2 L R L R L L R R q Q Q Q q u u ρ ρ ρ ρ ⎡ + ⎤ ⎡ ⎤ =⎢ ⎥= + = ⎢ ⎥ + ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ (3-29) 再找出B=B Q( )與C=C Q( )使得 U B Q Δ = Δ ； F C QΔ = Δ (3-30) 將(3-30)結合可得 1 ( ) F CB− U Δ = Δ (3-31) 再根據條件 3 求出近似 Jcaobian 矩陣 1 A=CB− (3-32) 為了滿足(3-30)，可以求得 1 2 1 2q 0 B q q ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦； 2 1 2 2 1 2 2 q q C a q q ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ (3-33) 再帶入(3-32)可得
2 2 0 1 2 A a u u ⎡ ⎤ = ⎢ − ⎥ ⎣ ⎦ (3-34)
u 為 Roe averaged velocity
L L R R L R u u u ρ ρ ρ ρ + = + (3-35) 因此可以用同樣方法得到以下物理量： L L R R L R v v v ρ ρ ρ ρ + = + (3-36) L L R R L R w w w ρ ρ ρ ρ + = + (3-37) L L R R L R H H H ρ ρ ρ ρ + = + (3-38) 1/ 2 [( 1)( 1/ 2 )] a= γ − H− V
(3-39) 其中 u 、v、w 分別代表 x 方向、 y 方向、 z 方向的速度。 H 、a 則分別為焓和音速。 (3-35)～(3-38)式中的U 以及L U 則是利用 MUSCL 法求出。 R
y x< 0 x= 0 x> x 0 圖 3-1 黎曼問題特徵值結構圖。 m λ 1 m λ − i λ 1 p λ + i λ p λ 1 λ 2 λ
3-3、LES model〈Smagorinsky model〉：
LES 將流體物理量區分為大尺度〈large-scale〉與次格點尺度〈subgrid-scale〉 兩部分。對於大尺度的物理量在 LES 中直接由 Navier-Stokes 方程式求解，而在次 格點尺度內的物理量則需要模式化。
目前大部分的次格點尺度流體剪應力模式〈subgrid-scale stress model〉是 以流體剪應力假設為基礎，其中最常被引用的模式為 Smagorinsky 次格點尺度流體 剪應力模式。 首先應用 LES 解決紊流問題的是 Deardorff[16]。他所模擬的為管道流體運動， 雷諾數相當的大。其最大貢獻在於闡述了利用 LES 方法來模擬三維流場的可行性。 由 Kolmogrov 的相似法則得知，小渦流的性質和雷諾數無關，因此小渦流部分 使用模式計算，而大渦流部分使用數值計算。其中渦流大小的劃分採用過濾方法 〈filter〉，用 filter function 將計算範圍確定：
( , ) ( , ) ( , ) U r t =U r t +U r t′ (3-39) 在(3-39)式中 ( , )U r t 為大尺度量值，而U r t′( , )為小尺度量值。 ( , )U r t 可經由過濾函 數而得到： ( , ) ( , ) ( ) U r t =
∫
U x t G r−x dx，其中 (G r− 為 filter function。 x) 其中較常見的 filter function 有下列幾種：(1) Gussian filter function：
( , ) 1 exp 6( 2 ) 6 i i i i i i x x G x x π ⎡ − ′ ⎤ ′ = × ⎢− ⎥ Δ ⎢⎣ Δ ⎥⎦ (3-40)
(2) Top hat filter function： 1 i Δ 2 2 i i i i i x −Δ <x′ < +x Δ ( ,G x x ′i i)= (3-41)
(3) sharp cutoff filter function： ( ) 2 ( , ) ( ) i i in i i i i i x x S G x x x x π π ⎡ − ′ ⎤ − ⎢ ⎥ Δ ⎢ ⎥ ⎣ ⎦ ′ = ′ − (3-42) 其中 Δ ：filter width i 利用(3-40)式對可壓縮三維之 Navier-Stokes 方程式過濾，且定義 f ρf ρ = 則可得： 3 1 2 0 F F F U t x y z ∂ ∂ ∂ ∂ + + + = ∂ ∂ ∂ ∂ (3-43) 其中， 1 2 3 ( , , , , )T U = ρ ρ ρu u ρu ρe ， 1 1 1 2 2 2 3 3 3 2 2 2 ( ) 2 i i i i i i i i i i i i ij i i u u u P A F u u P A u u P A T e P u A u x ρ ρ δ μ ρ δ μ ρ δ μ ρ μ λ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − ⎟ ⎜ ⎟ ⎜ ⎟ = + − ⎜ ⎟ + − ⎜ ⎟ ⎜ ∂ ⎟ ⎜ + − − ⎟ ⎜ ∂ ⎟ ⎝ ⎠ 1, 2,3 i ∀ = 。 (3-44) 其中λ ρ= C kp ， k 為 thermal diffusivity。 2 2 2 1 2 3 1/ 2 ( ) v e C T u u u ρ =ρ + ρ + + 。 1 2 ( ) 2 3 j i ij ij i j u u A u x x δ ⎡∂ ∂ ⎤ = ⎢ + − ∇ ⎥ ∂ ∂ ⎢ ⎥ ⎣ i ⎦， P=ρRT 。 現定義一物理量 subgrid-stress tensorϒ ，ϒ = −ij ρu ui j+ρu ui j。將其分離成兩部 分： 1 1 1 3 3 3 ij ij llδij iiδij τij iiδij ϒ = ϒ − ϒ + ϒ = + ϒ ， (3-45)
1 3 ij ij ll ij τ = ϒ − ϒ δ 。 由(3-45)式可將(3-44)式整理成以下的結果： 1 1 1 1 2 2 2 2 3 3 3 3 2 2 2 ( ) 2 i i i i i i i i i i i i i i i i ij j ij j i u u u A F u u A u u A T e u Q u A u x ρ ρ ϖδ τ μ ρ ϖδ τ μ ρ ϖδ τ μ ρ ϖ τ μ λ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − − ⎟ ⎜ ⎟ ⎜ ⎟ = + − − ⎜ ⎟ + − − ⎜ ⎟ ⎜ ∂ ⎟ ⎜ + − − − − ⎟ ⎜ ∂ ⎟ ⎝ ⎠ 1, 2, 3 i ∀ = 。 (3-46) 其中： ( ) ( ) i i i Q = −ρ e+P u + ρe+ϖ u (3-47) 2 2 2 1 2 3 1 1 ( ) 2 2 v ll e C T u u u ρ =ρ + ρ + + − ϒ 2 2 2 1 2 3 1 1 ( ) ( ) 2 2 v ll v C T u u u C ρ ρ ρ = − ϒ + + + 2 2 2 1 2 3 1 ( ) 2 v C u u u ρ ϑ ρ = + + + (3-48) 1 3 ll P ϖ = − ϒ (3-49) (3-48)式中 1 2 v ll T C ϑ ρ = − ϒ (3-50) 且由 P=ρRT和〈3-49〉式可得： 1 ( ) 2 v 3 ll R RT C ϖ ρ= + − ϒ 3 5 6 ll RT γ ρ − = + ϒ (3-51) 現在定義次格點尺度馬赫數〈subgrid Mach number〉 2
sgs M ： 2 2 ll sgs M c ρ ϒ = (3-52)
因為 2 c =γRT 所以(3-52)經移項後可以改寫為： 2 ll γMsgsP ϒ = (3-53) 在低馬赫數時(3-52)式可視為極小，因此在大部分的情況底下可視為忽略。所以由 (3-49)以及(3-50)式可知ϖ P、ϑ T。將此結果代入(3-46)式可得： 1 1 1 1 2 2 2 2 3 3 3 3 2 2 2 ( ) 2 i i i i i i i i i i i i i i i i ij j ij j i u u u P A F u u P A u u P A T e P u Q u A u x ρ ρ δ τ μ ρ δ τ μ ρ δ τ μ ρ τ μ λ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎜ + − − ⎟ ⎜ ⎟ ⎜ ⎟ = + − − ⎜ ⎟ + − − ⎜ ⎟ ⎜ ∂ ⎟ ⎜ + − − − − ⎟ ⎜ ∂ ⎟ ⎝ ⎠ 1, 2,3 i ∀ = 。 (3-54) 為 了 使 (3-54) 式 形 成 close form ， 必 須 對τij 以 及Q 進 行 模 式 化 。 使 用i subgrid-stress model 可得： ij tAij τ ρν (3-55) Pr t i p t i T Q C x ν ρ ∂ ∂ (3-56) 其中νt使用 Smagorinsky model 可將其模式化為： 1 2 2 1 2 ( ) ( ) 2 j i t s j i u u v C x x ⎡ ∂ ∂ ⎤ = Δ ⎢ + ⎥ ∂ ∂ ⎢ ⎥ ⎣ ⎦ 1 3 1 2 3 ( ) Δ = Δ Δ Δ (3-57) s C ：0.1 壁面紊流；0.25 等方向性紊流；0.15 介於上述兩者之間 Δ ：網格間距 1：X 2：Y 3：Z 3-4、Preconditioning 法： 為了提高程式的應用範圍，加入 preconditioning 法，讓程式不論在高速或低 速流體內，皆可獲得精確的計算結果。本程式採用 Weiss and Smith 的
preconditioning method [17]，方程式改變如下： 0 U F G H t x y z ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-58)
(3-58)為原始方程式，接著將保守形式（conserved variables）轉變成主要變數形 式（primitive variables），其形式如下： 0 p U F G H M t x y z ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-59) 其中 [ ]T p U = P U V W T ， M 為轉換矩陣： 0 0 0 0 0 0 0 0 0 1 p T p T p T p p T p T p u u U v v M U w w H u v w H C ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ∂ ⎢ ⎥ = = ∂ ⎢ ⎥ ⎢ ⎥ ⎢ − + ⎥ ⎣ ⎦ (3-60) 其中 p p ρ ρ =∂ ∂ ； T T ρ ρ =∂ ∂ 接著將(3-59)式的方程式乘上矩陣 K 2 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 ( ) 1 u v K w H V u v w ⎡ ⎤ ⎢ − ⎥ ⎢ ⎥ ⎢ − ⎥ = ⎢ − ⎥ ⎢ ⎥ ⎢− − − − − ⎥ ⎢ ⎥ ⎣ ⎦ (3-61) 再將 K 與 M 相乘 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 p T p KM C ρ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢− ⎥ ⎣ ⎦ (3-62) 將(3-62)式帶入(3-59)式，觀察連續方程式： ( ) 0 p p u v w t x y z ρ ρ ρ ρ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-63) 在理想氣體中可將(3-63)再表示成 2( ) 0 p u v w C t x y z γ ∂ +∂ρ +∂ρ +∂ρ = ∂ ∂ ∂ ∂ (3-64)
其中 C 為聲速 由(3-64)式可以看出，在等密度條件下，由於ρp為零，(3-63)式將變成 0 u v w x y z ρ ρ ρ ∂ +∂ +∂ = ∂ ∂ ∂ (3-65) 上式即為不可壓流的連續方程式。 綜上所述，可以得知只要改變(3-62)式中的ρp項，利用當地流場速度（local velocity）的倒數取代，即可轉換系統中的特徵值，藉此改變低速情況下流場的聲 速，使聲速與流場速度冪次級數(order)相同，系統不再受到 CFL 條件的限制，提高 程式的計算效率。 利用θ取代ρp項： 2 1 1 ( ) r p U TC θ = − (3-66) ε×Umax if u < × ε C Ur = u if ε× <C u < (3-67) C C if u > C 其中ε為一極小的值，約等於 5 10−，其主要是用來防止停滯點（stagnation point） 在計算時所造成的奇異點（singular point）現象。對於黏制性流體而言，U 必須r
大於流體的當地擴散速度（local diffusion velocity），因此U 還需加入下列限r
制： max( , ) r r U U x ν = Δ 將θ帶入(3-62)式後，可得到一新矩陣Γ nc 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T nc C θ ρ ρ ρ ρ ρ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ Γ = ⎢ ⎥ ⎢ ⎥ ⎢− ⎥ ⎣ ⎦ (3-68)
經過上述推導之後，方程式從(3-58)式轉變如下： ( ) 0 p nc U F G H K t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-69) 為了讓(3-69)式中的通量項再度轉換成保守形式，在乘上 1 K− 1 (K nc) Up F G H 0 t x y z −Γ ∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ (3-70) 根據(3-70)式，定義 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-71) 最後方程式簡化成如下形式： 0 p U F G H t x y z ∂ ∂ ∂ ∂ Γ + + + = ∂ ∂ ∂ ∂ (3-72) 由於方程式在時間項經過改變，因此必須重新推導 Roe 所提出的近似黎曼解。在 (3-22)式中，可以觀察到 1 2 i F + 項，是由 1 ( ( ) ( )) 2 F UR +F UL 的中央差分法加上為了解 決不連續面問題的 artificial viscosity term 1
2 A UΔ 所組成。加入
preconditioning 的方程式只需在 artificial viscosity term 做改變即可，其推 導如下：
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-73) 其中 1 1 AM KA DA KA − − Γ = × × 3-5、LUSGS scheme： 方程式(3-72)中的 Navier-Stokes 方程式在時間項方面遭到修改，因此利用修 改後的方程式來計算暫態結果並不恰當，因此本程式再加入 dual time stepping[18] 的模組，不僅讓程式在計算暫態結果方面較準確，更提高程式的效率，縮短計算時 間。
首先，先在原始 Navier-Stokes 方程式加入一虛擬時間項，稱為 artificial time term。方程式改變如下： 0 U U F G H t x y z τ ∂ +∂ +∂ +∂ +∂ = ∂ ∂ ∂ ∂ ∂ (3-74) 其中τ 即為 artificial time， t 為 physical time
接著在 artificial time term 加入 preconditioning 的方法： 0 p U U F G H t x y z τ ∂ ∂ ∂ ∂ ∂ Γ + + + + = ∂ ∂ ∂ ∂ ∂ (3-75) 最後對 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-76)
此處 k 為 artificial time 中的疊帶次數， n 為 physical time 的計算階數。 上述方程式，當 artificial time term 收斂至
1 0 k k p p U U τ + − Γ = Δ 時，方程式即會回復 到原始的 Navier-Stokes 方程式，並且包含著時間項，故程式可以計算暫態結果。 接著整理(3-76)式，先將其線性化 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-77) 其中 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 δ + − ⎡ ⎤ = ⎢ − ⎥ Δ ⎣ ⎦。 再將ΔUp項，放置於等號左邊，其餘則於等號右邊： 1 3 1 [ ( )] 2 k k k k x p y p z p p I M A B C U R t δ δ δ τ − − + Γ + Γ + + Δ = Δ Δ (3-78) 此處 1 3 4 ( ) ( ) 2 k n n k k k k x y z U U U R F G H t δ δ δ − − + = − − + + Δ ， I 為單位矩陣。 現在令 1 p p A = Γ− A Bp 1Bp − = Γ 1 p p C = Γ−C (3-79) 並將其分成兩部分： p p p A =A ++A − Bp =Bp++Bp− Cp =Cp++Cp− (3-80) 其中 1 ( ) 2 p p A A ± = A ±λ I 1( ) 2 p p B B ±= B ±λ I 1( ) 2 p p C C ± = C ± λ I (3-81) A λ 、λB、λC分別為A 、p B 、p C 中最大的特徵值。 p 由(3-81)式可導出：
1 , , 1, , , , 2 ( p p) ( p p i j k) ( p p i) j k i j k A U A + U A − U + + Γ Δ = Γ Δ + Γ Δ 1 1, , , , , , 2 ( p p) ( p p i) j k ( p p i j k) i j k A U − A + U − A− U Γ Δ = Γ Δ + Γ Δ 1 , , , 1, , , 2 ( p p) ( p p i j k) ( p p i j) k i j k B U B + U B − U + + Γ Δ = Γ Δ + Γ Δ (3-82) 1 , 1, , , , , 2 ( p p) ( p p i j) k ( p p i j k) i j k B U − A+ U − A− U Γ Δ = Γ Δ + Γ Δ 1 , , , , 1 , , 2 ( p p) ( p p i j k) ( p p i j k) i j k C U C + U C − U + + Γ Δ = Γ Δ + Γ Δ 1 , , 1 , , , , 2 ( p p) ( p p i j k) ( p p i j k) i j k C U C + U − C − U − Γ Δ = Γ Δ + Γ Δ 接著將(3-82)式代入(3-78)式可得 (L+ +D U)ΔUp =Rk (3-83) 其中， 1 1, , , 1, , , 1 1 1 1 ( p )i j k ( p )i j k ( p )i j k L A B C x y z − + + + − − − ⎡ ⎤ = −Γ ⎢ Γ + Γ + Γ ⎥ Δ Δ Δ ⎣ ⎦ (3-84) 1 1 , , , , , , 3 1 1 (( ) ( ) ) (( ) 2 p i j k p i j k p i j k D M A A B t x y − − ⎡ + − + = Γ + Γ ⎢ Γ − Γ + Γ Δ ⎣Δ Δ (3-85) , , , , , , 1 ( Bp )i j k) (( Cp )i j k ( Cp )i j k) z − + − ⎤ − Γ + Γ − Γ ⎥ Δ ⎦ 1 1, , , 1, , , 1 1 1 1 ( p )i j k ( p )i j k ( p )i j k U A B C x y z − − − − + + + ⎡ ⎤ = Γ ⎢ Γ + Γ + Γ ⎥ Δ Δ Δ ⎣ ⎦ (3-86) (3-83)式又可以近似為 1 ( ) ( ) k p L+D D− D U+ ΔU =R (3-87) (3-87)式可以用以下的步驟解出： 1. (L+D)ΔUp∗=Rk 1 ( k ) p p U ∗ D− R L U ∗ Δ = − Δ
2. (D U) Up D Up ∗ + Δ = Δ 1 p p p U U ∗ D U U− Δ = Δ − Δ 3. k 1 k p p p U + =U + ΔU 第三個步驟求出的 k 1 p
程式流程圖：
開始 設定邊界條件 MUSCL 設定參數 RoeDeal time stepping
Residual Preconditioning 小於 0.001 大於 0.001 到達預定時階 Output 是 否 Viscid terms
第四章 結果討論
為了驗證程式的正確性，本研究計算模擬channel flow的流場。本章首先驗證 利用第二章所產生的其初始條件正確性。其取法是利用Mare[5]中的方法。其次介紹 模擬的物理模式，包括流場的初始條件、物理條件。並與Kim[2]中所模擬的流場數 據做比較。 4-1初始條件確認 首先，先對產生的初始條件作確認。利用Mare[5]中方法所產生的初始條件。這 些初始條件皆依照所規定的Reynold stress R 和 u 、 v 、 w 的平均速度所產生。ij ij R 與u 、 v 、 w 的定義如下： ( , , ) ( , , , ) ( , , , ) ij i j R x y z =<u′ x y z t u′ x y z t > (4-1) 2 1 1 ( , , , ) t t u u x y z t dt t = Δ∫
(4-2) 2 1 1 ( , , , ) t t v v x y z t dt t = Δ∫
(4-3) 2 1 1 ( , , , ) t t w w x y z t dt t = Δ∫
(4-4) 其中Δ = − 。 t t2 t1 先對其Reynold stress 和平均速度作檢驗。如圖4-1~4-3分別為平均速度 u 、 v 、 w 。而其中黑線部分為所規定的數值，紅點部分為程式所計算出的數值。在平均速度 u 利用fully developed channel flow的理論值給定，而在平均速度 v 、 w 則 給定為零。各圖的一致性相當良好。圖4-4~4-9為Reynolds stress圖。圖4-4~4-6 為R 、xx R 、yy R 分別取根號後對磨擦速度 uzz τ作無因次化。各圖的一致性相當良好。 其中 uτ的定義如下： w uτ = τ ρ (4-5) w τ 為單位壁面所受的摩擦力，ρ為密度。
而圖4-7為R 並對摩擦速度 uxy τ之平方做無因次化。以上之數值是參考文獻[2] 中的數值。圖4-8和4-9為R 和xz R ，其所訂之數值為零。各圖的一致性相當良好。yz
因此由圖4-1~4-9可知，利用文獻[5]中所使用的方法可以準確的產生所要的紊流初 始條件。
圖4-1 平均速度 u 。 所定之數值， 模擬數值。
圖4-2 平均速度 v 。 所定之數值， 模擬數值。
v
y
圖4-3 平均速度 w 。 所定之數值， 模擬數值。
w
y
圖4-4 Reynolds stressR 對摩擦速度 uxx τ無因次化後數據 xx R uτ 。 所定 之數值， 模擬數值。 y δ xx R uτ
圖4-5 Reynolds stressR 對摩擦速度 uxx τ無因次化後數據 Rxx uτ 。 所定 之數值， 模擬數值。 y δ yy R uτ
圖4-6 Reynolds stressR 對摩擦速度 uxx τ無因次化後數據 xx R uτ 。 所定 之數值， 模擬數值。 y δ zz R uτ
圖4-7 Reynolds stressR 對摩擦速度 uxy τ無因次化後數據 Rxy2 uτ 。 所定之 數值， 模擬數值。 y δ 2 xy R uτ
圖4-8 R 。 所定之數值， 模擬數值。 xz y
δ
xz
圖4-9 R 。 所定之數值， 模擬數值。 yz y
δ
yz
4-2 channel flow
接下來為了驗證程式所模擬出的紊流流場正確性，本研究模擬了channel flow 的流場。其幾何外觀如圖2-1所示。 x 方向為stream wise的方向、 y 為垂直壁面方 向最後 z 為span wise的方向。 u 、 v 、 w 分別為其所對應的速度。
以下各項case A數據的取得為經過七個大渦流時間(large eddy turn over time)
τ Δ 後，再經過十個Δ 取平均。其中大渦流時間定義如下： τ uτ δ τ Δ = (4-6) 其中δ為管道的一半高度， uτ為摩擦速度。
Case B為case A將網格加密後。其case A、B兩者的物理條件均在表4-1中。 利用(4-2)式計算其 x 方向平均速度 u ，(4-5)式計算Reynaolds stressR 。並ij 將結果與文獻[2]中的比較。圖4-10為平均速度 u 對摩擦速度 uτ作無因次化以後，並 將其對 y+(y u y v τ ρ += )做圖。與理論值做比較，發現在靠近壁面附近的準確性較低。 且因為流量守恆，當壁面速度計算不正確時連帶影響整個速度曲線。因此平均速度 曲線與理論值差距較大。此為在計算中流場的viscous sublayer之內的網格較少， 是壁面計算不準確之最大原因。而圖4-11~4-13分別為R 、xx R 、yy R 無因次化後的zz 數據，並與Kim[2]中比較。發現在壁面處之數值與文獻[2]中的差距較大，且有震盪 的現象發生。離開壁面後與[2]之數值比較誤差不大。但在case B中，因為壁面網格 較密，因此震盪的現象消失。且其數值均與Kim[2]中的數值較接近。接著比較R 。xy 圖4-14與4-15為對不同數值做無因次化後的R ，並與[2]做比較。Case A中壁面處xy 的震盪現象依然存在，但在case B中震盪的現象並不存在，且以整體看來數值與[2] 相當接近。最後計算vorticity fluctuation並與文獻[2]比較。Vorticity定義如下： w= ∇× (4-7) V 其中V u i v j w k ∧ ∧ ∧ = + + ， i j k x y z ∧ ∧ ∧ ∂ ∂ ∂ ∇ = + + ∂ ∂ ∂ 。 因此先算出流場的平均vorticity，再求其vorticity fluctuation。
圖4-16~4-18為 x 、 y 、 z 方向無因次化後的vorticity fluctuation。比較後case A 壁面的不準確性依然出現。不過離開壁面後，整體的趨勢以及準確性都有不錯的表 現。Case B的準確性以及數值的穩定度皆比case A好很多。因此可知，在加密網格 後使得壁面的數值較穩定，且整體的準確性提高不少。
在本節中case A利用文獻[5]中的方法所產生的初始條件模擬channel flow。其 模擬出的流場性質分別為圖4-10~4-18。經過與文獻[2]比較後，準確性上有不錯的 表現。除此之外，本研究亦利用亂數產生初始條件來模擬channel flow，所模擬的 物理條件與case A相同。其流場模擬出的結果經平均後如圖4-19與4-20。圖4-19為 xx R 經無因次化後所得之數據。與圖4-11比對後可知，其order與文獻[2]及case A 所模擬出的數據相差甚大。在比對圖4-20與圖4-14，其結果也是與文獻[2]和case A 的數據落差很大。因此可知，在使用不同的初始條件進行模擬可得出不同的結果。 在(3-23)式中，隨著ε之改變會影響artificial viscous term A UΔ 在計算
時所佔的比例。當ε改變時計算之結果也會隨著改變。因此本研究模擬了當其他物 理條件皆不改變時，僅改變ε時所產生的變化。如圖4-21、4-22分別為R 及xx R 對xy 摩擦速度 uτ無因次化後的數據，並與文獻[2]做比較。其中物理條件為case A中的 條件，僅改變ε分別為0.03、0.5、1.0。由圖4-21、4-22可發現，當ε越大時其Reynolds stress也隨之變大。因此須不斷調整ε使其計算結果能接近正確值，此項為整個模 擬最艱難之處。
表4-1 流場物理條件。
Re Reτ uτ nx ny nz Δ x Δ y Δ z Δ t ε
case A 2700 180 0.086524 50 50 10 0.004 m 0.0013 m 0.0026 m 0.0005 s 0.03 Case B 2700 180 0.11248 80 100 20 0.002 m 0.0005 m 0.001 m 0.0005 s 0.1
