7.3. 以守恆變數法求解
7.3.2. 以遲滯密度法求解
7.3.2.1. 漸縮-漸擴噴嘴內流場
以中央差分與上風差分混合法執行三種不同背壓之一/二維漸縮-漸擴 噴嘴模擬測試，一維網格為100 CV，二維網格為 100×10 CV，測試結果分 述如下：
(1).一維漸縮-漸擴噴嘴模擬測試：
(A).在 Mref=1 時，k1必須小於 1 才不會發散且其值愈小(k1=0.675)，抑制 震波前後之振盪效果愈佳。
(B).當面上之密度近似法混合因子使用αρ=0.8 及速度項混合因子使用 0.6 時，可以得到較佳之結果，馬赫數分佈如圖7-223 所示。
(C).若加密網格至 200 個網格時，則可得到更佳的結果，時間間隔則需減 小為dt=1.E-5，馬赫數分佈如圖 7-224 所示。
(2).二維流場模擬之 K 值測試：
(Mref=1.、Δt=1E-4、ηu=1.0、ηv=1.0)
(A).Pb=0.870Po 時：ηp=0.2、αu=0.8、αv=1.0、αρ=0.8，k 由 1.0 調整至 0.6，計算結果詳如圖 7-225 至圖 7-229，當 k=0.6 時震波之捕捉效果 最好。
(B).Pb=0.769Po 時：ηp=0.2、αu=0.7、αv=1.0、αρ=0.8，k 由 0.9 調整至 0.5，計算結果詳如圖 7-230 至圖 7-234，當 k=0.6 時震波之捕捉效果 最好。
(C).Pb=0.645Po 時：ηp=0.1、αu=0.6、αv=1.0、αρ=0.8，k 由 0.8 調整至 0.6，計算結果詳如圖 7-235 至圖 7-238，當 k=0.6 時震波之捕捉效果 最好。
(D).由計算結果可知：
(a).當震波強度愈強時，CDS/UDS 之混合因子必須作調降，即 UDS 所 佔之比例增加，否則會造成發散現象。
(b).在本漸縮-漸擴噴嘴之流場模擬所使用之混合因子於 Pb 為 0.870 時 混合因子設為 0.8；0.769Po 時混合因子設為 0.7；0.640Po 時混合 因子設為0.6，可維持較佳之穩定度，而在相同測試條件狀況下，k 值通常取介於0.5~1.0 之間之值，在不發散之情況下，K 值愈小震 波之捕捉效果愈好。
(3).二維漸縮-漸擴噴嘴模擬測試：
(A).由上述之 K 值測試結果可以獲得最佳之計算情況，如圖 7-238 所示。
(B).在計算過程中發現有兩個主要參數會影響計算之收歛性，即壓力修正 方程之鬆弛因子及k 值，在使用此方法時要特別注意參數之調整特性。
由準一維及二維流場之計算結果顯示，不管低速或高速均與解析解相 當吻合，並對震波強度及位置之捕捉均有不錯的結果，證明本文之方法可 執行低、高速全速流場之計算。就數值振盪而言，遲滯密度法在一維流場 幾乎可得到無數值振盪之結果，而在二維流場亦可得到不錯的計算結果。
7.3.2.2. 流經下壁面圓弧之渠道內流場
使用與遲滯壓力法相同之計算網格，另在超音速部份由於流場中有複 雜之震波反射及交匯，使用本方法亦無法進行流場計算，尚待進一步探討 原因及謀求解決之道。
7.3.2.2.1. 次音速內流場之計算
入口馬赫數Min=0.5，分別搭配 SUPERBEE 與 Van Albada 限制函數進 行測試，次音速非黏性流場測試結果之壁面馬赫數分佈及馬赫數等值圖如 圖7-239 及圖 7-240 所示。不論是利用遲滯壓力法或遲滯密度法之守性變數 求解，所有參數包括Mref= 1、k=1 及 dm=2 等均不需作任何限制，且不論採 用何種通量限制函數，其所得之解及收歛情況(如圖 7-241)均相當。故遲滯 壓力法或遲滯密度之守性變數求解法對於次音速非黏性可壓縮流場有不錯 之結果。
7.3.2.2.2. 穿音速內流場之計算
入口馬赫數 Min=0.675，其邊界條件及網格系統與上述次音速流場相 同。面上對流速度及密度之近似值則使用 SUPERBEE 與 Van Albada 限制函 數差分法進行測試，詳細結果分述如下：
以SUPERBEE 限制函數法來進行測試，測試結果如下：
(1).當 Mref=1、 k=1 時，dm 值愈小數值振盪幅度愈小，如圖 7-242 所示。
(2).當 Mref及 dm 相同之情況下，k 值愈小則數值振盪幅度愈大，如圖 7-243 所示。
(3).當 k 及 dm 相同之情況下，Mref值愈小則數值振盪幅度愈小，如圖7-244 所示。
(4).由以上結果綜合研究顯示在穿音速流場 M=0.675，參數之設定 Mref=0.9，
k=1，dm=0.8 可得到較佳的結果，其與 Parameswaran [65]的計算結果相 近。其馬赫數等值圖及壁面馬赫數分佈如圖7-245 及圖 7-246 所示。
7.3.2.3. NACA 0012 翼型外流場
採用與前節原始變數法或特徵變數通量限制函數法相同的測試條件，
自由流馬赫數及攻角分別為(i) M∞=0.63、α=2°；(ii) M∞=0.80、α=1.25°；(iii) M∞=0.85、α=1.0°；(iv) M∞=1.2、α=0.0°。
以 O-型四邊形網格採用 Van Albada 及 SUPERBEE 限制函數進行測 試，在 M∞=0.63、α=2°時之參數設定為 Mref1=1、k=1、dm=2，SUPERBEE 通量限制子之計算結果如圖7-247 所示；在 M∞=0.80、α=1.25°時 Mref1=1、
k=0.5、dm=0.8，SUPERBEE 通量限制子之計算結果如圖 7-248 所示；在 M∞=0.85、α=1.0°時 Mref1=1、k=1、dm=1，Van Albada 通量限制子之計算結 果如圖7-249 所示；最後在 M∞=1.2、α=0.0°時 Mref1=1、k=1、dm=1，Van Albada 通量限制子之計算結果如圖7-250 所示。由計算結果顯示不論在次音速或超 音速流場均有不錯表現，尤其在穿音速流場對震波的捕捉有非常不錯的結 果。
第 8 章 結論
由前一章對各種流場以原始變數法、特徵變數限制函數法及守恆變數 法等求解類型之各種方式計算分析驗證結果，可以得到下列幾點結論：
1. 本文所發展之原始變數法之流場解子，使用 SIMPLE 算則之壓力修正法 結合正規化變數(NV)法或全變量消去法(TVD)之高階通量限制函數法，
應用重置變數無結構性網格，並採有限體積積分法來離散統御方程式，
藉由壓力修正方程式中之馬赫數控制因子，可自動隨流場特性而調整統 御方程式之型式，確實可求解低速至高速之全速流流場；在守恆變數法 之流場解子藉由遲滯壓力或遲滯密度來解決超音速場計算時之不穩定現 象，確實可求解非黏性低速至高速之全速流流場。
2. 本文發展之原始變數法解子，可依需求選用各種不同的對流項差分法，
包括CDS/UDS 混合法及 18 種高階限制函數法等；特徵變數限制函數法 解子，則有18 種限制函數可供選擇；守恆變數法之遲滯壓力或遲滯密度 解子，亦可選用與原始變數法相同的對流項差分法。可有效地滿足各種 不同的流場計算需求。
3. 原則上在次音速流場中，不論是採用原始變數法或守恆變數法求解時，
對流項之差分法使用 CDS/UDS 混合法或任何高階限制函數法，所得到 的結果均相當。而在穿音速或超音速具有震波的流場中，高階限制函數 法則具有抑制數值震盪的功效，在各種的限制函數中，SUPERBEE 的特 性具有較緊密(compact)的變化，雖準確度較高但數值振盪較劇烈；而 Van Albada 限制函數則具有較平滑的變化，故擴散性較大且數值振盪較輕 微，其他的情形則介於兩者之間。
4. 本文亦證明所發展的流場解子，有能力執行任意多邊形的無結構性網格 之計算。
5. 本文之方法採用高階之 TVD 或 NV 限制函數依流場梯度變化情形自動調
整中央差分法與一階準確上風差分之混合比例，來近似面上對流通量之 變數，從計算結果可以發現在震波前後均有抑制數值振盪的效果，雖有 些微的強度損失，但卻可得到令人滿意的震波解。
6. 使用特徵變數 MD1 之 Van Albada 限制函數法配合梯度修正限制法，在 震波之捕捉及計算效益之經濟考量下是較佳的選項。使用特徵變數MD2 具有較強的數值擴散現象，所以在Van Albada 限制法時幾乎可以完全抑 制震波前與後之數值振盪，對於加密網格則是不錯的選擇。即使是在 SUPERBEE 限制法時雖有輕微的數值振盪現象，仍不失為好的方法。
7. 以守恆變數法求解之遲滯密度法或遲滯壓力法而言，其在一維或二維之 噴嘴內流場、渠道及NACA 0012 非黏性流之次音速與穿音速流場計算，
對於震波位置與強度的捕捉均有不錯之結果，且對於震波前後的數值振 盪抑制效果亦相當良好，唯在超音速之複雜震波交匯之渠道內流場則面 臨困難，計算過程均產生發散，未來仍需再深入研究尋求解決之方法。
8. 由本文之計算結果可以發現面上密度之近似法，對程式的穩定性及解的 精確度有很大的關係。原則上若在沒有震波現象或梯度變化陡峭之不連 續流場，應採用完全的中央差分法來近似之；若在有震波之流場中，則 需混合少許的一階上風差分，則可消除震波不連續處之數值振盪，震波 的強度越強則一階上風差分的比例就要隨之增加。雖然亦採用 NVD 或 TVD 之限制函數加以調適，但在疊代計算過程中仍難以避免無預期的數 值劇烈變化，故即始使用限制函數差分法，密度的限制子仍需依流場特 性而加以限制，一階上風差分的比例愈高則愈穩定。
9. 本文發展之通量限制法可以很容易地使用各種 TVD 或 NV 限制函數，且 有兩種不同的變數基底來求得原始變數法之限制子，即原始變數或特徵 變數基底，兩種方法之優劣無一般性之區分，會隨不同流場之不同而異，
一般而言特徵變數基底法之收歛速度較快。
10. 梯度之計算方式對於震波強度及位置的捕捉有密切的關聯，若採用格點
平均法，則會有很強的擴散性；若採用相鄰網格內插法，則視近似差分 法之準確度而異，若在有震波區域之附近，二階以上之準確近似法容易 造成過與不及的預估值，有時甚至會導致解的數值振盪或疊代過程的發 散現象，故藉由二階準確近似之修正法及與相鄰網格之極值限制，不但 可維持高階準確，亦可抑制計算過程發散之效果。
11. 本文在網格之處理所發展之無結構性網格整合程式，具有二維或三維之 混合任意多邊形網格組合能力，有效地縮減複雜流場網格產生之時間，
並經由流場測試證明所整合網格之正確性及計算解子具有處理混合網 格之能力。
12. 未來仍需再進一步延伸至紊流模式及三維流場解子之發展，俾利應用於 各種飛行器外流場及引擎內流場相關議題之研究與分析。
第 9 章 參考文獻
[1] Jacob E. Fromm, ”A method for reducing dipersion in convective difference schemes”, Journal of Computational Physics, Vol. 3, pp. 176-189 (1968).
[2] Robert W. MacCormack, “The effect of viscosity in hypervelocity impact cratering”, AIAA Paper, 69-354 (1969).
[3] Jack R. Edward and Meng-Sing Liou, “Low-Diffusion Flux-Splitting Methods for Flows at All Speeds”, AIAA Journal, Vol. 36, No. 9, pp.
1610-1617 (1998).
[4] H. Luo, J. D. Baum and R. Lohner, “Extension of Harten-Lax-Van Leer scheme for flows at all speeds”, AIAA Journal, Vol. 43, No. 6, pp.
1160-1166 (2005).
[5] A. J. Chorin, “A Numerical Method for Solving Incompressible Viscous Flow Problems,” Jounal of Computational Physics, vol. 2, pp. 12-26 (1967).
[6] A. Rizzi and L. E. Eriksson, “Computation of Inviscid Incompressible Flow with Rotation,” J. Fluid Mech., vol. 153, no. 3, pp. 275-312 (1985).
[7] D. Choi and C. L. Merkle, “Application of Time-Iterative Schemes to Incompressible Flow,” AIAA Jounal, vol 23, pp. 1518-1524 (1985).
[7] D. Choi and C. L. Merkle, “Application of Time-Iterative Schemes to Incompressible Flow,” AIAA Jounal, vol 23, pp. 1518-1524 (1985).