第三章 數值方法
3.4 流場求解流程及參數設定
3.4.3 參數設定
1. 第一個網格距導緣之水平距離為 1.0*10-2 弦長，J 方向距壁面的 間隔（spacing）為 1.0*10-3。
2. 流場分為兩個區塊，每個區塊網格數71*61。
3. c(弦長)為1，
µ
（黏度）為5.0*10-7，Re =ρ Uc / µ
，k(紊流動能) 為10-4，ε
(紊流動能消散率)為 10-3、及U∞(入流速度)為1，ρ
/ρ
v= 1000，r(弧長半徑及邊長)為 10，Cdest的選取範圍為1.0*10+4 ~ 1.0*10+5，Cprod的選取範圍為1.0*10+2 ~ 1.0*10+3，c （壓力修正 值）為0.1。
4.
u，v，p，k， ε
，vf 之疊代收斂值為 10-5。5. 壁函數(y+)之範圍可為 30 至 800 之間,但在本文計算中大部份選 擇在30 至 250 之間。
第四章 計算結果與討論
本文分別計算了四種翼形：NACA 0012、NACA 4412、Eppler 以及 NACA 66（mod），並加以分析，而以NACA 66(mod)與文獻[32]和 Senocak and Shyy[27]
做一比較。
4.1 空化區長度的定義
翼形空化發生的位置，發展的過程會隨著攻角的改變而不同，而對於使 用理論計算的研究，不同的空化係數也會產生不同的空化現象，由於本研究 使用二相流法去模擬空化的發生，因此重點將放在空化區的長度上，而不著 眼於空化區的大小，這是二項流法的特性所致，不像小板法可以較為明確的 定義出空化區的邊界。如何定義空化區的長度如下(以 NACA 66(mod)在
α =
4°，σ
= 0.91 為例)：圖 4.1 為計算所得之翼表面 Cp 圖，圖中可以發現空化區的壓力係數都 降為吾人所設定的 0.91 值，這是因為空化現象所造成空化區內的壓力皆變為 水的飽和蒸氣壓所致。
圖 4.1 NACA66(mod)攻角 4 度的 Cp 圖
接著，放大 A 點的區域來觀察，可以發現轉折點約莫落在 x / c = 0.4 的 位置，同時將放大區域的 Y 座標從 Cp 改為以 Volume Fraction 當成變數，
如圖 4.2
圖 4.2 放大 A 點的 Volume Fraction 圖
在圖 4.2 中，x / c = 0.4 可看出所對應的 Y 值約為 0.41，於是在圖 4.3 的 Volume Fraction 等值線圖裡定義 contour 值為 0.41，如下：
圖 4.3 Volume Fraction 等值線圖
圖 4.3 中弧長所對應的 X 座標值約為 0.38，即為吾人所定義的空化區長 度。
4.2 實驗結果比較
針對與實驗文獻[32]的驗證，吾人以攻角為 2.3 度以及 4 度的 NACA 66(mod)為主，比較結果及表格如下：
圖 4.4 攻角為 2.3 空化係數為 0.734 之比較圖
圖 4.5 攻角為 4 空化係數為 1.0 之比 較圖
σ 空化區長
度(實驗)
空化區長度(本 研究計算) 0.734 0.22 0.24 0.621 0.89 0.86 0.463 1.0 0.92
σ 空化區長
度(實驗)
空化區長度(本 研究計算) 1.0 0.34 0.35 0.775 0.8 0.77 0.518 1.0 0.86 表 4.1 攻角為 2.3 空化係數為 0.734 表 4.2 攻角為 4 空化係數為 1.0 之比
之比較表格 較表格
從圖和表格中可知，使用二相流法的模擬結果上，空化區長度跟實驗結 果比較大致是一樣的，惟獨在超空化現象(Super cavitation)的情況有較明顯 的差別，吾人認為應該是在翼尾端的網格點分布太少以致所計算出來的結果 有如此顯著的誤差
4.3 理論結果比較
由於本文的理論基礎是參考 Senocak and Shyy [27]，所以將最後所得的 結果與之做一比較分析，觀察彼此間的差異，如下：
圖4.6 攻角為 4.0 空化係數為 0.84 之比較圖
圖4.7 攻角為 4.0 空化係數為 0.91 之比較圖
關於圖 4.6 以及圖 4.7 的差別，吾人認為是數值方法以及網格密度的差 別，本文在速度、壓力、及其他方程式的疊代上都使用較低階的上風差分法 (Upwind scheme)，同時本文的網格密度比 Senocak and Shyy[27]在大小上 約小一階(order)，才使得在翼尾端的模擬上有明顯出入，然而在空化係數的 表現上，本文就有相當接近的結果了。
4.4 密度變化
觀察密度在空化區和沒有空化時的差別，本文分別就攻角上的變化以及
σ
上的變化來比較密度的差異，如下：4.4.1 攻角改變
由下圖可知，隨著攻角的變化，
ρ
m min就有越來越小的趨勢，然而由此圖 也發現，在x / c 等於 0 到 1 之間，密度都有變化，亦即有空泡的出現。
圖 4.8 NACA 0012 改變
α
的ρ
m圖4.4.2 空化係數改變
圖 4.9 NACA 0012 改變
σ
的ρ
m圖圖 4.9 說明空化係數越小則混合密度變化的範圍越大，總結圖 4.8 與 4.9 剛好證明了，攻角越大以及空化係數越小都表示越容易空化，則混合密度的 變化範圍就越大，反之，則變化範圍越小。
第五章 結論與未來展望 5.1 結論
本程式對於二維水翼空化問題的模擬，在與理論計算及實驗結果的比較 上，有合理的一致性，也有一些誤差需要再加以改進。對於不合理的假設，
程式也能反應不收斂和無法計算的結果，例如當壓力曲線已經呈現很平緩的 狀況時，即使再改變 Cdest與 Cprod，計算結果將不收斂，因此本程式將可提供 實驗者做快速的測試，以了解可能發生的空化現象，減少資源的浪費；而在 水翼性能方面的分析，可提供螺槳設計者對於空化的情況以及升阻力特性，
有一初步的概念。以下為翼形的簡單結論：
1. 隨著 Cdest與 Cprod的增加，可以發現 Cp的曲線越平緩，越接近空 化係數，所以這兩個參數值如何取決，對於最後的計算結果影響 甚大，但是模擬中也發現，若一開始疊代的值過大，將使得程式 無法收斂下來，因此在本文中，是以小的值慢慢循序漸進至大的 值。
2. 當二維水翼處於小攻角時，空泡多由導緣發生，而且空泡起始位 置由導緣逐漸後退，唯獨NACA 4412 正好相反。
3. 對於二維水翼，隨著空化係數的增加，則空化長度漸減，且阻力 減小，升力增加。
4. NACA 4412 與 Eppler 的比較中，在小攻角情況時，Eppler 翼形提 供較大的升阻比。
5. 密度和體積比率
α
l在空化區裡，已有明顯的改變，唯獨在尾緣的 地方有較長的擴散現象發生，這是需要改進的地方。5.2 未來展望
本程式對於空化初始已模擬成功，唯對於基本現象仍有些許的誤差，正 是未來需要再加以改進的地方，首先將網格點做一較完整的規劃如下：
圖 5.1 （上）本文之計算網格，（下）加密後之計算網格
圖 5.1（下）為在翼形的上表面等間距分佈 200 等分，並在 J 方向 0.1c 上將 格點集中，再均勻的散出去，為的是在 suction side 壓力降的地方能計算準
確，但是網格加密後，數值方法也必須加以改進，否則在流場參數的計算上 會有震盪的情形出現。
接著，二維翼形的計算發展完成之後，將進一步轉到三維翼形上面，
畢竟在實務上三維仍是主要的發展趨勢，利用在二維翼計算上得到的經驗和 心得，希望在三維翼形能有滿意的結果。
參考文獻
[1]
M. J. Lighthill, “A New Method of Two Dimensional Aerodynamic Design,”Technical Report ,Rand M2112, Arc, 1945
[2]
M. B. Giles, M. Drela, “Two-Dimensional Transonic Aerodynamic Design Method,” AIAA Journal, Vol.25, No.9, September 1987[3]
I. De Falco, R. Del Balio, A. Dell Cioppa, E. Tarantino, “Breeder Genetic Algorithms for Airfoil Design Optimisation ,” Evolutionary Computation , 1996, Proceeding of IEEE International, Conference, pp.71-75,1995[4]
辛敬業, “應用邊界元素法之二維翼面設計方法,” 第八屆造船及輪機工程 研討會論文集, pp.314-318,1995[5]
R. F. van den Dam, J. A. van Egmond, J. W. Slooff, “Optimization of Target Pressure Distributions,” Special Course on Inverse Method for Airfoil Design for Aeronautical and Turbomachinery Application, AGARD Report No.780 , Reference 3, Nov 1990[6]
J. A. van Egmond, “Numerical Optimization of Target Pressure Distributions for Subsonic and Transonic Airfoil Design (Inverse) and Optimization,”AGARD Conference Proceedings No. 463, Reference 17, March 1990.
[7]
I. De Falco, R. Del Balio, A. Della Cioppa and E. Tarantino, “A Parallel Genetic Algorithm for Transonic Airfoil Optimisation,” Proceedings of the 1995 IEEE International Conference on Evolutionary Computation. Part 1 (of 2), pp.429-434, Nov 29-Dec 1 1995[8]
Shigeru Obayashi, Susumu Takanashi, “Genetic Algorithm for Aerodynamic Inverse Optimization Problems,” Genetic Algorithms Engineering Systems:Innovations and Applications, 12-14 September 1995, Conference Publication No. 414, IEE, 1995
[9]
邱永裕, “以基因演算法與平行運算進行翼型優化,” 國立成功大學航太 工程學系, 碩士論文, 1999[10]
K. Deb, “ An efficient constraint handling method for genetic algorithms,”Computer methods in applied mechanics and engineering, Vol.186, pp.311-338, 2000
[11]
I. H. Abbott, A. E. von Doenhoff, “Theory of Wing Sections Including a Summary of Airfoil Data,” Dover publishing, Inc, 1959[12]
Kinnas, S.A. and Fine, N.E. “Theoretical Prediction of the Midchord and Face Unsteady Propeller Sheet Cavitation”, Proceedings of 5thInternational Conference on Numerical Ship Hydrodynamics, Hiroshima, Japan. (1989) [13]
Kinnas, S.A. “Leading-edge Corrections to The Linear Theory of PartiallyCavitating Hydrofoils”, Journal of Ship Research, Vol.35, No.1, March, pp.15-27. (1991)
[14]
Kinnas, S.A. “Inversion of The Source and Vorticity Equations for Supercavitation Hydrofoils”, Journal of Engineering Mathematics, pp.349-361. (1992)[15]
Fine, N.E., and Kinnas, S.A. “A Boundary Element Method for The Analysis of The Flow Around 3-D Cavitating Hydrofoils”, Journal of Ship Research, Vol.37, No.3, September, pp213-224. (1993)[16]
Kinnas, S.A., and Mazel, C.H. “Numerical Verus Experimental Cavitation Tunnel (A Supercavitating Hydrofoil Experiment)”, Trans. of ASME, Vol.115, pp.760-765. (1993)[17]
陳建宏,翁英傑,”二維全空化水翼的邊界元素法解析 ”第十屆中國暨 輪機工程研討會論文集, pp.204-211 , (1997)[18]
郭真祥,黃維信,何俊賢,”以源流分佈及高斯積分法則為基礎之高階邊 界 元 素 法 ”, 中 華 民 國 第 十 九 屆 全 國 力 學 會 議 論 文 集 ,pp.389-398,(1995)
[19]
郭真祥,田曉真,”空蝕水翼設計與性能解析(II)”, 第十一屆中國造船 暨輪機工程研討會暨國科會成果發表會,pp. 29-34,(1998)[20]
Chen, Y. and Heister, S.D. “A Numerical Treatment for Attached Cavitation”,Journal of Fluids Engineering, Vol. 116, pp. 613-618. (1994)
[21]
Deshpande, M., Feng, J. and Merkle, C.L. “Numerical Modeling of the Thermodynamic Effects of Caviation”, Journal of Fluids Engineering, Vol.119, pp. 420-427. (1997)
[22]
Delannoy, Y. and Kueny, J.L. “Cavity Flow Predictions based on the Euler Equation”, ASME Cavitation and Multi-Phase Flow Forum, 109, pp. 153-158.(1990)
[23]
Chen, Y. and Heister, S.D. “Modeling Hydrodynamic Nonequilibrium in Cavitating Flows”, Journal of Fluids Engineering, Vol. 118, pp.172-178,(1996)
[24]
Singhal, A.K., Vaidya, N., and Leonard, A.D., “Multi-dimensionalSimulation of Cavitating Flows Using a PDF Model for Phase Change,”
ASME Paper FEDSM973272, Proc. of 1997 ASME Fluids Engineering
Division Summer Meeting. (1997)
[25]
Merkle, C.L., Feng, J., and Buelow, P.E.O., “Computational Modeling of theDynamics of Sheet Cavitation,” Proc. 3
rdInternational Symposium on Cavitation, Grenoble, France. (1998)
[26]
Kunz, R.F., Boger, D.A., Stinebrung, D.R., Chyczewski, T.S., Lindau, J.W., and Gibeling, H.J., “Preconditioned Navier-Stokes Method for Two-phaseFlows with Application to Cavitation,” Comput Fluid, 29, pp.849-875. (2000) [27]
Senocak, I., and Shyy, W., “A Pressure-Based Method for TurbulentCavitating Flow,” J. Comput. Phys., 176, pp. 363-383. (2002 )
[28]
Boussinesq J., Essai Sur La Theorie Des Eaux Courantes, Vol.23, Mem.Pressentes Acad. Sci., pp46, (1877)