以守恆變數(ρu,ρv,ρw)求解之有限容積傳輸方程式可簡單表示如下：
( )
S Convection Source
transient Diffusion
d d q d
劃底線之部份以顯項來處理，將它放入源項中；另αφ代表混合因子(blending factor)，其值介於 0 與 1 之間(0=UDS；1=CDS)，或通量限制子(flux limiter)
( )r
同樣使用over-relaxed 法[81]增加計算的穩定性，可推得擴散項之離散式如 下：
2
CDS/UDS 之混合因子；若使用第四章所述之高階限制函數法，則以限制函 數ψ( )r 取代αφ，如下式所示：
為了使疊代過程較穩定，離散程式之鬆弛處理及線性方程求解方式，
與原始變數求解法相同。
5.4. 壓力、密度及速度偶合關係式
以守恆變數(ρu, ρv, ρw)求解動量方程，而壓力修正方程式與不可壓縮流 之推導方式相同，只是密度的變化不是單獨處理，而是與守恆速度(ρV )的 變化結合在一起，在求解可壓縮流時，為了使不可壓縮流橢圓形型式之壓 力修正方程產生在穿音速及超音速時之消散機制，必須使用遲滯密度來修 正動量方程式中面上的速度，或在動量方程式及壓力修正方程式中均使用 遲滯壓力。
5.4.1. 面上質量流率的處理
以不可壓縮流而言，質量流率之修正量係由面上速度之修正而求得，
說明如下：
**f [ (* * ) ]f
m = ρ V +V S′ i (5.29)
*
f f f
m′ =ρV S′i (5.30)
然而對可壓縮流而言，質量流率與垂直面上之速度分量Vn及密度均有 關，相對於不可壓縮流，密度的變化顯得非常大，所以質量流率之修正必 須包含速度與密度之修正量，若以守恆變數求解時，其壓力修正方程式與 不可壓縮流一樣，只是其速度之修正式已隱含密度之修正量，說明如下：
**f [( )* ( ) ]f f ( )*f f ( )f f
m = ρV + ρV ′ iS = ρV iS + ρV ′ iS (5.31) ( )
f f f
m′ ρ ′
⇒ = V iS (5.32)
面上的密度與速度修正量可表為：
(1). 若使用遲滯壓力法時：
** *
( )f ( )f ( )f f f
P f
V V V P D P
ρ ′ = ρ − ρ = −⎛⎜ΔΩA ⎞⎟ ∇ = − ∇′ ′
⎝ ⎠ (5.33a)
(2). 若使用遲滯密度法時：
** *
( )f ( )f ( )f f f
P f
V V V P D P
ρ ′ = ρ − ρ = −⎛⎜ΔΩA ⎞⎟ ∇ = − ∇′ ′
⎝ ⎠ (5.33b)
由於密度已和速度結合為守恆變數，因此在傳輸方程式中對流項的離 散式則必須將質量流率mf改為體積流率qf，通過面上的質量流率mf 之計算 方法與式(3.48a)相同，只是若使用遲滯壓力法時則式(3.48a)中之壓力必須改 用遲滯壓力，而通過面上的體積流率qf之計算方法表示如下：
(1). 若使用遲滯壓力法時：
[( ) ]
f f f f C C P f PC
m =ρ V Si −B p − p − ∇p iδ (3.34a)
f
f S f f
f
q d m
=
∫
V S V Si ≈ i = ρ (5.34b)(2). 若使用遲滯密度法時：
[( ) ]
f f f f C C P f PC
m =ρ V Si −B p − p − ∇p iδ (3.34c)
f
f S f f
f
q d m
=
∫
V S V Si ≈ i = ρ (5.34d)5.4.2. 壓力修正線性代數方程 將式(5.33)代入式(5.32)，可得
2
為相鄰網格係數AC之和，與不可壓縮流之型式相同，故必須藉由遲滯密度 或遲滯壓力之人工消散項之機制，方能模擬雙曲線型式之超音速流場特 性。式(5.37)之壓力修正線性方程式亦使用雙共軛法矩陣解子(BICG) [83]。
5.4.3. 邊界條件
邊界條件的處理與原始變數求解法相同皆可直接引用。
5.4.4. 求解程序
本文發展之守恆變數求解法則有兩種類型，即遲滯壓力法及遲滯密度 法，求解步驟說明如下：
1. 給初始入口速度Vin、面上質量流率m*f ，猜初始壓力場 p*P或p*P(遲滯壓力 法時)、溫度場TP*，而密度場ρ*P由狀態方程式求出。
2. 計算動量方程式之係數，並計算每一個網格加總之動量方程式殘值，利 用矩陣解子求出內點速度場(ρV)*P，並由邊界條件更新邊界點速度。
3. 計算面上速度(ρV)*f 及更新面上質量流率m*f 。
4. 計算壓力修正方程式之係數，並計算每一個網格加總之質量殘值，利用 矩陣解子求出內點壓力修正值p′，並更新邊界點之壓力修正量。
5. 修正網格中心點速度(ρV)*P*、面上質量流率m*f*及網格中心之壓力場p**P 或
**
pP (遲滯壓力法時)。
6. 由停滯溫度TO或總焓(HO)及速度場VP**[ (= ρV)** ρ*]更新溫度場TP**。 7. (1)若使用遲滯壓力法：在動量方程式及壓力修正方程式中所使用的壓力
為遲滯壓力p，在求解密度前須先將遲滯壓力 p轉換為真實壓力 p後，
再由狀態方程式更新密度場ρ*P*與體積流率q**f 。(2)若使用遲滯密度法：
由狀態方程式更新密度場ρ*P*，再利用真實密度ρ*P*求出遲滯密度ρ**P ，計 算通過面上的體積流率q**f 。
8. 若搭配使用限制函數法，則計算每一個控容面的限制因子ψ(r)，以作為 下一疊代之對流項差分式時使用。
9. 將求得的p**P 或p**P 當作新的p*P或p*P，重覆步驟 2-9，直到動量方程式之加 總殘值及壓力修正方程式之加總質量殘值滿足收歛條件(小於 1.e-3)為 止。
第 6 章 網格產生
無結構性網格在二維流場多採用 3 邊形或 4 邊形，在三維則為 4 面體 或 6 面體，當然亦有任意邊形或混合式網格，但其產生方式就較複雜。故 無結構性網格之計算解子(solver)均採用有限體積法(FVM)求解。
流場計算分析工作依程序可概分為三個部份，一是前處理即網格產 生、二是流場計算即計算解子、三是後處理即是流場性質顯示。無結構性 網格之產生有很多技巧及方式，然而，因其是流場計算分析工作之前處理，
產生之網格必須能依照一定的格式正確地輸入計算解子以進行流場計算，
雖然已具有簡單外形之 3 邊或 4 邊形網格之產生工具，鑑於對複雜外形網 格之產生須投入大量時間，以區塊網格方式產生後再合併處理，網格常需 配合計算需求而調整亦相當費時；且目前在市面上雖有許多發展相當成熟 且可快速產生網格之商業化工程分析軟體，可以快速產生二維或三維之複 雜外形無結構性網格，例如MSC.Patran、FLUENT、CATIA、SolidWork 或 其他Advancing Front 網格產生程式等，但因其格式均搭配其自身之計算解 子(solver)，對吾人自行發展計算程式所需之網格格式而言，無法直接使用，
然而其亦有一些問題，如網格之疏密控制或任意邊形之混合網格需求等。
為了解決網格產生問題，吾人發展一套網格產生處理程式，可以將不 同方式產生之不同形狀(3 邊(面)至 6 邊(面))、相鄰區塊之網格，整合為一組 滿足吾人所發展之計算解子所需網格，大大地縮短前處理時間。
建置完成網格整合產生處理程式之後，只需要獲得無結構性網格最基 本的元素，即「格點(vertices) x, y, z 座標」及「構成網格之格點編號」，就 可經由本套網格產生處理程式輕易且快速地獲得所需之網格資料，程式架 構詳參附錄B。
以下就簡單列舉三個案例說明：
(1).自由流流經一圓柱之計算域網格：
圓柱直徑為 1，遠場邊界為半徑 50 的外圓，在圓柱表面先佈置 8 層之
四邊形網格，分為上下兩個區塊來產生網格，外圍區塊則以Advancing fron 方法產生之三邊形網格，靠近圓柱周圍及尾流區之部份則作局部加密，俾 能清楚捕捉流場的特性。
當各別區塊之網格產生後，藉由點座標及網格的構成點編號所得到的 簡單網格資料，利用吾人自行發展的IGTP23D 網格整合程式進行網格的組 合，組合成一完整的計算域網格及產生計算解子所需之網格輸入資料，如 圖6-1 所示。
(2).自由流流經一三角柱之計算域網格：
三角柱為一邊長為 2 的正三角形，遠場邊界為半徑 50 的外圓，在三角 柱表面附近先佈置20 層之四邊形網格，分為 6 個區塊來產生網格，其中上 方之區塊則分別由2 個次區塊先行組成後再映像至下方並對稱於 y=0，外圍 區塊則以Advancing fron 方法產生之三邊形網格，同樣的在靠近圓柱周圍及 尾流區部份則作局部加密，俾能清楚捕捉流場的特性。
當各別區塊之網格產生後，藉由點座標及網格的構成點編號所得到的 簡單網格資料，利用吾人自行發展的IGTP23D 網格整合程式進行網格的組 合，組合成一完整的計算域網格及產生計算解子所需之網格資料，如圖6-2 所示。
(3).NACA 0012 外流場計算域網格：
在NACA 0012 翼型附近佈置 O-型四邊形網格，外圍遠場則佈置三邊形 粗疏網格，再利用吾人自行發展的 IGTP23D 網格整合程式進行網格的組 合，組合成一完整的計算域網格，如圖6-3 所示。
第 7 章 結果與討論
本文發展之流場解子(solver)共有三種求解類型，原始變數求解法(含 CDS/UDS 混合法及通量限制函數法)、特徵變數通量限制函數法、守恆變數 求解法(含遲滯密度及遲滯壓力法)等，每一種類型又可選用多種不同的對流 項差分法，為了測試是否可準確地預測自低速不可壓縮流到高速可壓縮流 之全速(all speed)流場，以下將就此三種求解類型分別進行非黏性流與黏性 流之流場測試。
首先以原始變數求解法求解非穩態之統御方程式，就一維與二維漸縮-漸擴噴嘴，來對震波之捕捉及對數值振盪之抑制效果進行測試；接著以二 維下壁面圓弧之渠道內流場，來進行次音速、穿音速及超音速流場之測試，
測試的方法包括CDS/UDS 混合差分法、原始變數之各種高階準確限制函數 差分法等，再從中選擇較佳之方法去進一步執行其他求解類型之測試，由 測試結果顯示SUPERBEE [106]限制函數具有較緊密(compact)的變化，雖準 確度較高但數值振盪較劇烈；而Van Albada [103]限制函數則具有較平滑的 變化，故擴散性較大且數值振盪較輕微，因此後續特徵變數通量限制法，
僅就SUPERBEE 與 Van Albada 此兩種代表性的限制函數進行測試；另以守 恆變數求解時使用遲滯壓力法及遲滯密度法等，亦搭配SUPERBEE 與 Van Albada 通量限制函數差分法進行流場計算測試。在擴散項或擬似擴散項均 使用二階準確之中央差分法。
本文以原始變數法求解之測試流場：黏性流包括低速空穴流、流經圓 柱之低速流外流場、流經NACA 0012 翼型外流場、雙喉部噴嘴內流場等；
非黏性流則有一/二維漸縮-漸擴噴嘴、流經下壁面圓弧之渠道內流場、流經 NACA 0012 翼型外流場等。以特徵變數通量限制法求解之測試流場：黏性 流包括低速空穴流、流經圓柱之低速流外流場、流經NACA 0012 翼型外流 場、雙喉部噴嘴內流場等；非黏性流則有流經NACA 0012 翼型外流場、流 經圓柱之高速流與流經三角柱之高速流外流場等。以守恆變數求解之測試
流場包括非黏性流之一/二維漸縮-漸擴噴嘴、流經下壁面圓弧之渠道內流 場、流經NACA 0012 翼型外流場等。
測試範圍包括內/外流場，黏性流體從低速不可壓縮流到高速可壓縮 流、非黏性流體則從次音速、穿音速到超音速流場，來充分證明本文發展 之流場解子具有全速流流場計算能力。三種求解類型之計算結果分述如后。