中 華 大 學 碩 士 論 文

中 華 大 學

碩 士 論 文

題目：Sanchez-Lacombe Equation of State 於兩 相流之初步探討

系 所 別：機械與航太工程研究所 學號姓名：M09508004 莊 詠 成 指導教授：牛 仰 堯 博 士

中華民國 九十七 年 七 月

中文摘要

本文的目的在建立一個計算二維可壓縮流氣、液混合數值模式 之程式。以二維保守形式之 Navier-Stokes 方程式做為統御方程 式，在空間離散中，則以混合型通量分離法 AUSMD 配合 MUSCL 方程 式以提高精度，再以具 TVD 效應之三階 Runge-Kutta 法處理時間離 散，使程式能在數值穩定與計算上求得較佳的平衡。另外，選擇 Sanchez-Lacombe 狀態方程式模擬氣、液混合相變化區。本文數值 驗證的例子，首先以一維震波管凝結現象解析來驗証數值模式之可 行性，隨後求取二維注油器相變化情形。

Abstract

In this study, the main purpose is to develop a compressible two-phase flow code with a saturation model. Here, a two-dimensional compressible mixture type Navier-Stokes equations is adopted as the governing equations. To keep good stability and efficiency at numerical approximations, the advective upwinding splitting scheme AUSMD with MUSCL type interpolation is used to evaluate the spatial differencing, also a three-order TVD type Runge-Kutta scheme is chosen to perform the time discretization. In addition, the Sanchez-Lacombe equation of state is studied and applied to model the liquid-vapor phase transition. In the numerical test cases, the validation of a one-dimensional shock-condensation tube and a two-dimension injection nozzle are performed to investigate the current numerical approaches.

致謝

首先感謝學生的指導教授 牛仰堯教授這兩年來的教導，對學 生在學業的學習上有很大的幫助，在程式撰寫方面，指導學生ㄧ步 ㄧ步的理解程式組織與架構，使學生能順利的完成論文，在此致上 萬分的謝意。

另外，也要特別感謝本實驗室的林於金學長、林宜勳學長、徐 唯淵學長以及邱勇憲同學、白凱元同學與許廷旭學弟，在這段時間 學業上的幫助，最後還要感謝家人在背後默默的支持，是學生能完 成論文的動力來源，感謝大家。

目錄

中文摘要………I 英文摘要………II 致謝………III 目錄………IV 圖與表目錄………VI 符號說明………VIII

第一章 緒論

1.1 前言………1 1.2 文獻回顧………6

第二章 數值模式

2.1 統御方程式………14 2.2 狀態方程式………15 2.3 數值通量分離法………25

第三章 數值結果

3.1 一維震波管凝結現象解析………29 3.2 二維注油器空蝕現象解析………33

第四章 總結

參考文獻

圖與表目錄

圖目錄：

圖 1 狀態方程式所得之壓力-密度對應圖……… 18

圖 2 狀態方程式所得之壓力-密度對應圖……… 19

圖 3 狀態方程式所得之密度-溫度對應圖……… 19

圖 4 Gibbs free energy-密度對應圖……… 20

圖 5 二氧化碳壓力-體積-溫度三相圖………20

圖 6 在時間 150μs 時的壓力與位置描述圖……… 31

圖 7 在時間 150μs 時的密度與位置描述圖……… 31

圖 8 在時間 150μs 時的蒸氣體積佔有率與位置描述圖………… 32

圖 9 在時間 150μs 時的速度與位置描述圖……… 32

圖 10 注油器之幾何外型比例結構圖……… 35

圖 11 注油器之計算網格……… 35

圖 12 維納收縮示意圖……… 36

圖 13 二維注油器流場局部放大之流線分佈圖……… 36

圖 14 注油器在時間 0.78S之體積佔有率局部放大圖……… 39

圖 15 注油器在時間 0.78S之壓力局部放大圖……… 39

圖 16 注油器在時間 0.78S之速度分佈圖……… 40

圖 17 注油器在時間 0.98S之體積佔有率分佈圖……… 40

圖 18 注油器在時間 0.98S之壓力分佈圖……… 41

圖 19 注油器在時間 0.98S之速度局部放大圖……… 41

圖 20 注油器在時間 1.28S之體積佔有率分佈圖……… 42

圖 21 注油器在時間 1.28S之壓力分佈圖……… 42

圖 22 注油器在時間 1.28S之速度分佈圖……… 43

圖 23 注油器在時間 1.755S之體積佔有率分佈圖………43

圖 24 注油器在時間 1.755S之壓力分佈圖………44

圖 25 注油器在時間 1.755S之速度分佈圖………44

圖 26 注油器在時間 2.00S之體積佔有率分佈圖……… 45

圖 27 注油器在時間 2.00S之壓力分佈圖……… 45

圖 28 注油器在時間 2.00S之速度分佈圖……… 46

圖 29 注油器在中心線之速度-入出口壓力差比較圖………46

圖 30 注油器之計算網格……… 47

圖 31 注油器細管後半段速度向量局部分布圖……… 48

圖 32 注油器細管前半段速度流線局部放大分布圖……… 48

圖 33 注油器體積佔有率局部分佈圖……… 49

圖 34 注油器體積佔有率局部放大分佈圖……… 49

表目錄：

符 號 說 明

流場混合聲速 a

氣體相之積體比 α_g 液體相之積體比 α_l 流體混合密度 ρ 氣相之密度 ρ_g 液相之密度 ρ_l 對比密度 ρ 參考密度 ρ^* 亞穩分域的對比密度 ρ _S 流體總能量 E 流體總內能 e

氣相之內能 e_g 液相之內能 e_l 流場總焓值 h

理想氣體的每單位質量焓 h_I 流場靜壓力 p 對比壓力 P

亞穩分域的對比壓力 PS

參考壓力 P^* 氣體常數 R 絕對溫度 K 溫度 T 對比溫度 T 參考溫度 T^* 亞穩分域的對比溫度 TS

流體分子重量 M_W 百萬分之一秒 μS X 方向之流體速度 u

Y 方向之流體速度 v

黏滯係數 ν 尺寸參數 r

第一章 序論

1.1 前言

在真實流場中，流體大部分以多相的形態存在。「相」是一種最普遍 的自然現象，在自然界中的物質基本上可定義為三種形態：固相、液 相和氣相。兩相流一般是指成分由單一種或多種元素所構成。根據不 同相態或元素狀態，兩相流可以分類為：氣-固兩相流、液-固兩相流、

及氣泡流等；當物質的相態多於二種，則稱之為多相流，其中氣泡流 廣泛存在於自然界和工程設備中，如油氣燃燒的熱傳與輸送現象以及 鍋爐等裝置中所產生的液-氣流動等。而氣-固兩相流、液-固兩相流 主要是應用於顆粒分離和除塵、水利工程中的泥沙運動和環境工程中 煙塵對空氣的污染、懸浮體的燃燒、含沙水流、血管流等。幾乎所有 的工程領域都會遇到一些在特定條件下多相流體運動所引起的複雜 問題。

目前使用數值模擬是研究兩相流體相當普遍的重要方法之ㄧ，主 要目的是用來探討流場當中各相流體的壓力、運動速度、溫度、密度、

體積佔有率(固、液、氣相)、流動方式、相和相之間的交互作用與轉 變以及相與壁面間相互作用(如氣泡或顆粒的阻力與傳熱傳質)，觀察 當中的壓降、動量傳遞、熱傳、化學反應、相變化等，如管內的液- 氣兩相流，隨著兩相流體的速度變化，相應產生的氣泡流、霧狀流等 形態。

根據林於金[1]，兩相流的理論分析與單相流相比困難得多，描 述兩相流的數值方法與狀態方程式目前尚未發展完成。如今大部分學 者們採用的理論是兩種簡化模式的類型：第一種為「均相模式」，是 將兩相介質看成是一種混合得非常均勻的混合物，然後再假設處理單 相流體的概念和方法仍然適用於兩相流，並且合理的假設其物理性質 及傳遞性質；第二種為「分相模式」，將單相流的概念和方法分別用 於兩相系統的各個相，並且同時考慮兩相之間的相互作用。以上兩種 模式的應用仍然存在不少困難，主要影響的因素在於兩相流體間的質 量、動量和能量傳遞過程是一種非常複雜的傳輸現象，而且亦涉及到 一些基礎學科，如：相變化、流體力學、熱力學和熱傳學等，但是在 計算技術發展的推動下亦有相當進展。

相變化和臨界現象在物理學中是相當困難的領域。自1869年T.

Andrews [2]發現臨界點，1873年Van Der Waals [3]所提出的非理想氣 體狀態方程式以來，對相變化的實驗和理論研究已經有一百多年的歷 史 ， 從 早 期 的 傳 統 熱 力 學 、 統 計 力 學 ， 到 最 新 的 重 正 化 群 (renormalization)理論。相變化的成因是溫度和壓力，然而發生相變 化時，是一種形態的物質轉變成另一種形態的物質(例如液態水轉變 為水蒸氣、鐵磁與反鐵磁)，必然存在著物質相變化的質傳作用，同 時伴隨著潛熱釋放或吸取相變化所產生熱交換的過程。而臨界點是真 實流體在三相圖中一個特殊的點，為飽和液-氣共存的一個點，此時 兩相的狀態都是相同。超過臨界溫度以後，則無法分清是液態還是氣 態。目前也還沒有一種狀態方程式可以完美的描述相變化，這個領域 依然還有許多未解的謎。

在 1894 年，Reynolds[4]是最早開始以科學研究空蝕現象，他發 現管路中氣泡的生成與消失。之後，在 1902 年發生英國驅逐艦 Cobra 的螺旋葉片在使用數個小時後就發生嚴重的損壞，之後在水力機械上 也看到相同的現象。在當時認為螺旋葉片材料的剝落是由於海水腐蝕 所造成的，但是經過實驗證明在蒸餾水中運動的物體也會出現類似的 剝蝕，因而確認這種現象只是受到機械力衝擊的結果。但是根據之後 的研究分析，上述兩種因素對葉片的損壞都有關係，因此空蝕現象才

開始受到人們的重視與研究。

當流體加速通過物體表面的區域，造成其壓力因某種原因而突然 下降至與該區域液體溫度相應的氣化壓力以下時，使的部分液體在此 區域蒸發，形成流體中的「氣泡(Cavities)」（或稱空泡），此一過程 稱之為「空化現象」。氣泡隨著流體進入壓力較高的區域時，失去存 在的條件而突然有崩裂的傾向，原本氣泡周圍的流體運動使局部區域 的壓力驟增。如果流場中不斷形成、長大後的氣泡在固體壁面附近頻 頻崩裂，壁面就會遭受巨大壓力的反復衝擊，從而引起材料的疲勞破 損甚至表面剝蝕，此一過程稱之為「空化剝蝕現象」，簡稱「空蝕現 象(Cavitation)」。這種現象常常伴隨著能量的釋放並且引發流場的不 穩定。因此一個週期不規則的氣泡產生和成長程序就此發生，然而並 非所有空化現象都會造成物體的損壞，只有不穩定的空化現象才會，

如非穩態流動中出現的空化現象或者封閉氣泡的末端，才會引起空蝕 現象。因此，空蝕現象往往出現在物體的局部區域。空蝕現象與物體 受固體微粒或液滴衝擊而造成的損壞是不相同。而且空蝕發生的區域 都是極其複雜的兩相流體結構，在某些尺度下是極不穩定的。空蝕區 的結構是取決於物體表面的形狀和粗糙度的變化以及操作的條件等 等。空蝕的程度則以物體的重量的降低、體積的減少、表面粗糙度的

變化和穿孔數目作為基準。

空蝕現象是包含流體動力學、材料學和物理化學的複雜現象。空 蝕現象的產生常常會導致物體大量效能損失，因此對物體經常要做額 外的保養。在外流場時經常發生於水泵、水輪機和船舶螺旋槳的葉片 表面，以及高水頭泄水建築物的局部表面上，都必須承受因空蝕現象 所導致的不良影響，如對機械表面造成侵蝕，造成整體性能的表現下 降 。 而 內 流 場 則 經 常 發 生 於 引 擎 中 燃 燒 室 的 注 油 孔 (Injector Orifice)、氣渦輪機、工業用的燃燒爐、甚至噴墨印表機的噴墨頭等 工業產品，也因為空蝕現象的產生而有著燃燒不完全，燃燒爐損毀、

或著色位置發生錯誤等諸多不良後果，因此應想辦法作適當的防範和 減少空蝕現象造成的損壞。

空蝕氣泡(Cavitation Bubble)依型態可歸納為兩種：第一種為「蒸 氣空蝕(Vapor Cavitation)」，流場中的壓力降低到比液體的飽和蒸氣 壓更低，此時大量的液體可迅速氣化為蒸氣，充滿在低壓區中。第二 種為「空氣空蝕(Gaseous Cavitation)」，當流場中的壓力隨著流速加 大而降低至某一程度，其壓力尚未達到液體的飽和蒸氣壓時，然而此 時原本溶於液體中的部份氣體，因為壓力降低，導致其從液體當中被

釋放出來，因而有氣泡產生。

要估算空蝕現象，需要模擬從液相到氣相或者從氣相到液相的質 量傳遞率，並且了解空蝕過程的物理認知。兩相流體在數值模擬上難 度極高，在數值方法中必需同時處理氣、液、混合三種相態的聲速、

壓縮性、密度、能量、壓力等等，而且氣、液兩相之聲速相差約四倍，

加上相變化時存在的數值震盪之問題，以及傳統的 Navier-Stokes 統 御 方 程 式 在 求 解 兩 相 流 體 時 ， 所 遭 遇 到 的 複 雜 數 值 模 式 ， 然 而 Sanchez-Lacombe 狀態方程式在模擬相變化過程中並無法降低數值 的不正常震盪，因此加深了計算上的困難度，但仍舊以此一狀態方程 式做簡單的分析，藉由計算得到流場性質的分布情形，以期能初步預 測出空蝕現象所發生的位置與時機。

1.2 文獻回顧

計算流體力學(Computational Fluid Dynamics 簡稱 CFD) 的發 展大約起始於 1950 年，主要內容是由流體力學、數學、數值方法以 及電腦科技的結合而成，傳統上，CFD 方法目前已經被分成兩種主

要的類型：第一種為「密度為基本變量(Density Based)」法，用來解 決可壓縮流方程式並且適合應用在氣體動力學(Gasdynamic)的發 展；第二種為「壓力為基本變量(Pressure Based)」法，適合解決在 原始變量的不可壓縮流方程式。

為求解保守形式的 Navier-Stokes 方程式，在空間離散上選擇使 用迎風插分法(Upwind Differencing Method)，因為迎風插分法較能 準確模擬在不同相過程中可能存在的流體反應。一般在多相流場問題 包含分離運輸方程式(Separate Transport Equation)，在每個單獨的相 和模型分析中間相的質量、動量和能量。目前科學研究發展各種類型 的迎風算則提出針對兩相流模式上的求解，其中最常被採用的是混合 模式的方法。在混合模式之算則當中，常見的有通量差分分離法(Flux Difference Splitting Scheme 簡稱 FDS)和通量向量分離法(Flux Vector Splitting Scheme 簡稱 FVS)、Roe 類型的方法[5]以及修改典型的迎 風算則(如 Godunov 類型的算則)[6]等。典型運用這些算則的代表有 Romate[7]、Toumi[8]、Tiselj 和 Petelin[9]以及 Masella 等人[10]。

其中，FDS 是利用矩陣的計算，FVS 是運用純量的計算；因此，FVS 比 FDS 有效率。有許多學者欲結合 FVS 與 FDS，因為這種方式可同 時具有 FDS 準確捕捉線性波的設計與 FVS 的簡單本質特性，且又不

增加複雜性的混合法被大量提出，以 Niu[11]和 Edwards 等人[12]的 研究較為完善；Niu 以混合通量類型的通量分離法模擬多相流的數值 計算，而 Edwards 則以齊次平衡的兩相流模式模擬真實流場，並求解 的相變化發生情形。Niu 和 Edwards 等人數值方法的特性，在於使用 的算則在結構上是簡易的 Euler 方程式。Liou 所發表的 AUSM[13-14]

法為其中較為成熟的代表之一，探討與 Euler 方程式有關的多相流數 值解析。Evje 和 FJelde[15-16]考慮更以簡單的等溫兩相流模式，架 構出求解分離質量保守方程與混合動量方程式。

迎風算則在低速流時收斂慢而且答案較不精確，而低擴散通量分 離法，例如 AUSM＋、AUSMDV[17]和 LDFSS，基於波分解的特性 選擇傳統的迎風算則。這種方法顯示結合更多高準確技術與穩定性和 簡單的通量向量分離法。在求解可壓縮多相流場時，由於相變化時狀 態方程式都無法真正模擬並且達到完美的境界，造成 Navier-Stokes 方程式在計算上有其困難之處，而且在計算特徵系統亦有其困難的地 方，因此使用 Roe 或者 Godunov 法變成極為麻煩。而使用 AUSM 來 求解多相流已被証實為非常有效率之事，Chang[18]和 Niu 利用 AUSM 法求解多相流都取得不錯的成果。

早期的低擴散通量分離法是假設流體是一種理想氣體或者混合 物。低擴散通量分離法，對一般的狀態方程式，可以描述單相為液體、

氣體或者一種任ㄧ特定物質的超臨界的流體行為。Toumi[8]利用齊次 平衡模型模擬管子在沸騰的過程，Yabe[19]發展統一的模型模擬在固 體、液體以及氣體的流體動力學。

通常統御方程式中時間項的積分以疊代(Iteration)的方式隨時間 趨進以接近解析解(Analytical)，也就是所謂的時間趨近法[20]。

Jameson[21-22]發展出使用 Runge-Kutta 之多重步階法以處理通量平 衡，常用於處理 Navier-Stokes 之 Runge-Kutta 法，可以細分為二階、

三階、四階。使用高精度之 Runge-Kutta 法雖然較為精準，但是所需 時間將花費更多，尤其是在求解多相流場時將更為明顯。不過，為平 衡計算時間與精度，將使用具 TVD[23]效應之三階 Rung-Kutta 法。

根據 Knapp[4]的文獻記載，空蝕現象的發現始於 1894 年。1902 年英國的驅逐艦在行進間螺旋槳運轉時，在螺旋槳表面產生一層近似 真空的薄膜，這薄膜不僅影響到螺旋槳推進的效率，也增加了不必要 的油耗，此真空薄膜即為空蝕現象。Nurick[24]發現，當流體流入孔 口轉角處為圓形轉角時造成的空蝕現象，引發流體與孔壁分離導致流

場不穩定，使得下游提早發生霧化，而對於流體流入孔口轉角處為方 形 轉 角 時 造 成 的 空 蝕 現 象 ， 對 流 體 與 孔 壁 分 離 現 象 較 不 明 顯 。 Hiroyasu[25]接續研究發現，當噴流的速度大於臨界速度時，在孔洞 入口處會產生空蝕氣泡(Cavitation Bubble)，此時靠近孔口內的流 場，因為氣泡而產生紊流的現象在結合液柱表面的擾動後，噴流液柱 碎裂的距離將縮短；反之，若噴流的速度過大，此時將呈現流場跳離 的效應(Hydraulic Filp)，不但拉長噴流液柱碎裂的距離，同時霧化的 效果也會變的不明顯。Koivula[26]研究指出，噴嘴的孔口細管的長 度直徑(L/D)比和孔口轉角處的形狀對流體的流動特性有重要的影 響，當孔口轉角處較尖銳時，由於維納收縮現象(Vena contracta)的影 響，會造成流場截面積縮小導致流體的流速增加和壓力降低，因此當 壓力低於流體的飽和蒸氣壓時，空蝕現象便開始產生。

Schmidt[27-30]等人發展出計算二維、暫態可壓縮流的空蝕數值 模式，特別是研究有關高速微小尺度的噴嘴。在程式中也加入三階震 波捕捉，避免密度在液體和氣體交界面時所產生之擾動。傳統的氣- 液二相流數值模式，無法有效處理高速可壓縮流在低壓區所造成的空 蝕現象 Richard 和 Jean[31]則發展出一套計算高速可壓縮流、混合模 式之空蝕數值方法，其藉由改變保守形式之 Euler 方程式，以達到消

除相變化時所產生之數值誤差，並且有取得良好的成果。

Arcoumanis[32]發現，當高壓注油器(Injector Nozzle)以極高的 速度，使流體流入燃燒室內，由於注油器流道內壓力極高，因此一旦 有流體進入狹窄的流道後，造成局部的壓力低於飽和蒸氣壓，此時所 產生的空蝕現象占了非常重要的地位，因為這種突然產生的氣化現象 會影響油滴霧化的混合參數、流體噴射速度和油滴變形等，這些現象 會造成燃燒不完全使得未氣化的油滴在高熱下產生碳化，而漸漸的在 注油器和燃燒室發生嚴重積碳造成阻塞，使的油耗增加、產生爆震等 不良後果。

本實驗室關於空蝕現象的研究，以三維鈍頭體在不可壓縮流、空 蝕模式且不含相變化[33]和二維混合模式含相變化[34]，均有取得不 錯的成果。然而在狀態方程式的選擇上是分別在氣體、氣液混合、液 體各選用ㄧ個狀態方程式組合起來，成為ㄧ組可以描述流體在氣體、

氣液混合、液體時的線性狀態方程式，因此希望能有一個狀態方程式 能單獨描述流體在氣體、氣液混合、液體時的狀態。而在 2000 年時，

Edwards [12]嘗試使用 Sanchez-Lacombe 狀態方程式[35]以齊次平衡 兩相流模型模擬液-氣兩相過渡區，並使用低擴散通量分離法在液態

二氧化碳流過一個毛細管噴嘴模擬空蝕現象和液-氣凝結震波。而 Sanchez-Lacombe 狀態方程式是一個能單獨描述真實流體狀況的非 線性方程式，因此本文的研究目的將利用 Sanchez-Lacombe 狀態方程 式與本實驗室所發展的數值方法來驗證此一狀態方程式的穩定性與 準確性。

第二章 數值模式

在本研究中首先計算流場之空蝕現象，為方便求解非穩態流場，

本文在時間離散上採用具有 TVD 效應之三階 Runge-Kutta 法；空間 離散則使用高階 MUSCL+Limiter Function 以提高精度；而在一維和 二維程式使用 AUSMD 含 MUSCL 類型 minmod limiting interpolation 作為主要之通量分離法撰寫。

本章節將分為三個部份來敘述數值模式，分別為統御方程式、狀 態方程式與數值通量分離法，其中在狀態方程式的選擇方面，以往本 實驗室林益民、林於金、徐唯淵在選擇狀態方程式時，通常在液相、

氣相、兩相混合區時分別使用 stiffened 狀態方程式[36]或者 Cole [37]

所發展的狀態方程式、理想氣體狀態方程式、經驗公式，將真實流體 的狀態使用三個方程式描述，然而 Sanchez-Lacombe 狀態方程式[35]

則能單獨的以一個非線性方程式來描述真實流體氣相、液-氣混合相 與液相的熱力特性。

2.1 統御方程式 (Governing Equations)

在本研究中，流場為二維保守形式之 Navier-Stokes 方程式：

V V

F G

Q F G

t x y x y

∂ ∂

∂ ∂ ∂

+ + = +

∂ ∂ ∂ ∂ ∂ {2-1}

保守向量Q和無黏性通量F G, 定義如下：

u

Q v

E ρ ρ ρ

⎡ ⎤

⎢ ⎥

⎢ ⎥

=⎢ ⎥

⎢ ⎥

⎣ ⎦

；

2

( )

u u p

F uv

E p u ρ

ρ ρ

⎡ ⎤

⎢ + ⎥

⎢ ⎥

=⎢ ⎥

⎢ ⎥

⎣ + ⎦

； ₂

( )

v G vu

v p E p v ρ

ρ ρ

⎡ ⎤

⎢ ⎥

⎢ ⎥

=⎢ + ⎥

⎢ ⎥

⎣ + ⎦

黏性通量F G_v, _v定義如下：

0

xx V

xy

xx xy x

F

u v q

τ τ

τ τ

⎡ ⎤

⎢ ⎥

⎢ ⎥

= ⎢ ⎥

⎢ ⎥

+ −

⎢ ⎥

⎣ ⎦

；

0

yx V

yy

yx yy y

G

u v q

τ τ

τ τ

⎡ ⎤

⎢ ⎥

⎢ ⎥

= ⎢ ⎥

⎢ ⎥

+ −

⎢ ⎥

⎣ ⎦

其中：p、u、v、ρ、E和τ 分別為壓力、x 方向的速度分量，y 方向 的速度分量、密度、總能及剪應力。

為求解液-氣二相模式，各數值簡單定義如下：

l l g

g

ρ α ρ

α

ρ = +

( )

E=ρ e+2u

l l g

g

e e

e = α + α

1 = α

+ α

其中：下標l和g分別代表液體和氣體。

2.2 狀態方程式 （State equation）

Sanchez-Lacombe[35]狀態方程式是一種古典的流體晶格分子理 論，此狀態方程式的前身是統計熱力學的 Gibbs free energy，經由晶 格公式與 Sterlingr 近似法推導而成。此種狀態方程式類似於立方的 狀 態 方 程 式 ， 有 點 類 似 Van Der Waals [3] 狀 態 方 程 式 或 者 Peng–Robinson 狀態方程式，並可以為極性物質(polar substance)，例 如水和辛烷等，提供更好的預測臨界點的能力。所謂地極性物質是指 因為每種原子搶電子的能力都不相同(其中以氟最強、氧次之)，所以 造成分子表面電子雲分佈不均勻。當搶電子的力彼此互相平衡時，電

子雲就會均勻分佈於分子表面因而形成。

狀態方程式的選擇對於空蝕現象有著極為重大的影響，因為越能 真實描述流體的相變化性質，數值模式在模擬流體的空蝕現象發生地 時 機 、 範 圍 、 氣 泡 的 流 動 等 趨 勢 便 越 準 確 。 圖 1 到 圖 3 為 Sanchez-Lacombe 狀態方程式所繪製。圖 1 顯示一條二氧化碳在低於 臨界溫度的壓力-密度等溫線對應圖。圖 2 則是繪製數條二氧化碳在 不同溫度下的壓力-密度等溫線對應圖。圖 3 則是繪製數條二氧化碳 在不同壓力下的密度-溫度等壓線對應圖。圖 4 則是繪製二氧化碳的 Gibbs free energy-對比比容對應圖。圖 5 則沒有按造比例，大略繪製 的二氧化碳壓力-體積-溫度三相圖。

圖 1 在溫度為 283K的壓力-密度圖，可以清楚的看出純蒸氣的 區域，而且此區的壓力與密度變化曲線是近似於線性；純液體的區 域，則是由於較大的壓力變化引起密度的改變。在純蒸氣區與純液體 區之間產生強烈的壓差，其結果將會導致高壓區與低壓區之間產生一 個蒸氣和液體的兩相混合共存區，而在此兩相共存區之中，液體與氣 體將會是可壓縮性的。在兩相區內壓力的相應密度分別為點 A、B 和 C。A 和 C 點分別代表飽和蒸氣和飽和液態，而 B 點在物理上是無意

義的。而兩相區在 D 和 E 兩點有局部的壓力臨界值。利用飽和蒸氣 壓與臨界點壓力的確定蒸氣和液體的亞穩分域(spinodal)曲線，並把 兩相區分成蒸氣亞穩態(metastable)區、不穩定和液體亞穩態區。在 液體和蒸氣亞穩分域的點之間可以藉由系統平衡，將帶有蒸氣和液體 逸能(或者 Gibbs energies)取得相等的值，此時的壓力稱之為飽和蒸 氣壓P_v，可以藉由疊代 Sanchez-Lacombe 狀態方程式作為一個溫度函 數和透過齊次平衡方法模擬液-氣兩相變化區的計算得到。由圖 2 為 壓力-密度圖，可以發現當溫度越高，兩相共存區的存在將越來越不 明顯，到溫度 400 K 時，則可看到一條線性的直線，這是因為當流體 溫度超過二氧化碳的臨界溫度 304.25 K，此時流體會處於一種既不 同於氣態也不同於液態和固態的一種新流體狀態，在超臨界條件下，

流體的性質發生了極大的變化，其密度、介電常數、粘度、擴散係數、

電導率和溶劑化性能等都會發生改變。圖 3 為密度-溫度圖，可以知 道在一大氣壓下二氧化碳一直是維持在氣相狀態，在 25 atm 時，出 現了液-氣二相變化，當壓力超過臨界壓力 73 atm，流體的性質發生 了變化，所以在 100 和 500 atm 時，則無法明顯的區別出液-氣二相 變化。圖 4 為在壓力和溫度分別為 3.31 MPa、263K時的 Gibbs free energy-對比比容圖顯示出，G的兩個相對較小值分別代表為氣體和液 體，d G<0表示有功消耗，d G>0表示環境對系統供應功，d G=0則表

示功無損失亦無生成即為平衡狀態，此時即為流體的飽和狀態。圖 5 為二氧化碳的三相圖，曲線 1 為昇華和凝華曲線；曲線 2 為凝結和氣 化曲線；曲線 3 為凝固和熔解曲線，點 B 為三條曲線相交點，是固、

液、氣三相共存於相同溫度與壓力分別為 216K、5.2 atm 的一個點，

稱之為三相點；點 C 則為液、氣二相共存點，此時的溫度與壓力分別 為 304K、73 atm 將達到最高，稱之為臨界點。

圖 1 狀態方程式所得之壓力-密度對應圖

圖 2 狀態方程式所得之壓力-密度對應圖

圖 4 Gibbs free energy-對比比容對應圖

A. 單相狀態：

為求解可壓縮流，本文採用 I. C. Sanchez and R. H. Lacombe[35]

所發展之狀態方程式，其表示如下：

^{P= −ρ2−T⎡}⎢⎣^{ln 1}

( )

其中：

*

* * * *

P MW

T P

T P r

T P ρ RT

≡ ≡ =

T^*、P^*、ρ^*是參考值，T、P、ρ分別是對比溫度、壓力、密 度，r是尺寸參數，R是通用氣體常數，M_W流體分子重量，需 要注意的是0≤ <ρ 1。

表 1 不同物質的溫度、壓力、密度參考值

其它數值的參數定義：

Critical point density：

1

c 1 ρ = r

+ {2-4}

Critical point temperature：

T_c =2rρ_c² {2-5}

Critical point pressure：

( )

2 1

ln 1 1

c c c c c

P T

ρ ^⎡ ρ ^⎛ r^⎞ρ ^⎤

= − − ⎢⎣ − + −⎜⎝ ⎟⎠ ⎥⎦ {2-6}

Specific entropy：

^*

( )

*

ln 1

( , ) _I( ) P 2 1

h T h T T

ρ ρ ρ

ρ ρ

⎧ ⎡ − ⎤⎫

⎪ ⎢ ⎥⎪

= − ⎨⎪⎩ + ⎢⎣ + ⎥⎦⎬⎪⎭

{2-7}

Speed of sound：

² ( ) ( ) ( )

T T

T T

h P P h h

a h P

ρ ρ

ρ ρ

ρ

⎡ ⎤

− ⎣ − ⎦

= − {2-8}

Spinodal density：

2 ²_{S S} 2 T^S _S T^S 0

T r r

ρ +^⎛⎜ − − ^⎞⎟ρ + =

⎝ ⎠ {2-9}

Spinodal temperature：

2 (1 ) 1(1 )

s s

s

s s

T

r

ρ ρ

ρ ρ

= −

+ −

{2-10}

Spinodal pressure：

( ) ( ) ( )

( )

2 2

2 1 ln 1 1 1 1

1 1

S S S S S S

S

S S

P r

r

ρ ρ ρ ρ ρ ρ

ρ ρ

⎛ ⎞

− − + −⎜⎝ ⎟⎠ − +

= − + − ⎛ ⎞⎜ ⎟⎝ ⎠

{2-11}

Saturation points and coexistence line：

G(ρ_v, , )T P =G(ρ_l, , )T P {2-12}

其中：

P (¹ 1) ln(1 ) ¹ln

G T

ρ ρ r ρ

ρ ρ

⎡ ⎤

= − + + ⎢ − − + ⎥

⎣ ⎦

h_I 是 理 想 氣 體 的 每 單 位 質 量 焓 ， 下 標 S 是 指 亞 穩 分 域 (Spinodal)，亞穩分域的密度和壓力值都要大於零。

B. 齊次平衡兩相狀態

流體的密度是介於飽和值和平衡狀態方程式之間，則數值參數定 義為：

entropy：

ρ ρh( , )T =ρ_v( )T α ρ_v( , ) ( )T h T_v +ρ_l( )T α ρ_l( , ) ( )T h T_l {2-13}

vapor volume fraction：

( )

( , )

( ) ( )

l v

v l

T T

T T

α ρ ρ ρ

ρ ρ

= −

− {2-14}

liquid volume fraction：

α ρ_l( , ) 1T = −α_v {2-15}

Speed of sound ：

1_{2 2 2}

( ) ( )

v l

v v l l

a a T a T

α α

ρ ⁼ ρ ⁺ρ {2-16}

其中：下標l和v分別代表為液體和氣體。

2.3 數值通量分離法(Numerical Flux Splitting)

在計算多相流場時，為避免計算特徵系統之困難，本研究中使用 Y. Wada and M.-S. Liou[13]所發表之 AUSMD。

AUSMD 在求解多相流場已為主流通量分離法之一，首先定義介 面上之質量通量：

( ρ u )

= ρ

u

+ ρ

u

其中：L 與 R 分別為介面上之左邊與右邊。

速度分離u_L⁺和u_R⁻可表示如下：

) 2 1 2 (

4 )

( ² L L

L L

L L

u u u

u a

a

u u +

−

⎥ +

⎦

⎢ ⎤

⎣

⎡ +

+ −

+ =α α {2-18}

) 2 1 2 (

4 )

( ² R R

R R

R R

u u u

u a

a

u u −

−

⎥ +

⎦

⎢ ⎤

⎣

⎡ −

− −

−

− =α α {2-19}

其中：

( )

2 1

2 /

1 aR aL

a = + {2-20}

a為交界面上之共用音速。

於此，α_L與α_R有多種選擇。Wada and Loiu 則利用壓力與密度的 比率來建立此二項數值：

2

2

L L

L R

R R

L R

p

p p

p

p p

α ρ

ρ ρ

α ρ

ρ ρ

⎧ ⎛ ⎞

⎪ ⎜ ⎟

⎪ = ⎝ ⎠

⎪ ⎛ ⎞ +⎛ ⎞

⎪ ⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

⎪⎨

⎪ ⎛ ⎞

⎪ = ⎜ ⎟⎝ ⎠

⎪ ⎛ ⎞ ⎛ ⎞

⎪ ⎜ ⎟ +⎜ ⎟

⎪ ⎝ ⎠ ⎝ ⎠

⎩

{2-21}

壓力通量則與AUSM⁺所定義的相同：

^{1/ 2}

1/ 2 L L

R R

M u a M u

a

⎧ =

⎪⎪⎨

⎪ =

⎪⎩

{2-22}

壓力之分離如下：

( 1) (22 )

, 1

4

( )

2 ,

P M M

if M

P P M M

otherwise M

±

⎧ ± >

= ⎨⎪⎪

⎪ ±

⎪⎩

∓

{2-23}

介面上之壓力則簡單的使用平均數來取得：

最後，介面上之通量可表示如下：

_{1/2 1/2} ( )_1/2 ( ) _1/2 2

) 1 (

) 2(

1 u u P

F = ρ φ_L +φ_R − ρ φ_R −φ_L + {2-25}

其中： ^{φ =}

[

]

[

]

第三章 數值結果與討論

3.1 一維震波管

在第一個例子中，將考慮一維震波管的凝結現象。一開始先假設 一根長達一公尺的直管在左半邊填滿了初始密度為 1200 kg/ m³的高 壓液態水；並且在右半邊則填滿了初始密度為 500 kg/ m³的液-氣二 相混合區。兩種流體最初的溫度都為 300 K。二相混合區由質量比 分別是由 50 %的液體和 50 %的氣體所組成。數值驗證則是在 201 個 格點數上進行。此外為了保持計算過程的穩定性，選擇了 CFL 為 0.5。

在圖 6 到圖 9 顯示時間為 150μs 的壓力、密度、氣體體積佔有率和 速度性質結果圖。

隨著時間步進，高壓液態區開始往低壓液-氣二相混合區移動，

造成速度由 0 m/s增加到約 250 m/s，此時在交界面處因為高壓與低 壓相遇造成壓力和密度產生劇烈的變化而形成震波，震波也同時由直 管左半邊的高壓液態區向右半邊之低壓混合區傳遞，相應產生的擴散 波也由右半邊混合區向左半邊的液態區前進。因為高壓往低壓區移動

的關係壓縮液-氣二相混合區，使的密度由 500 kg/ m³漸漸增加至密 度約 1050kg/ m³，由於震波的傳遞造成原本液-氣二相混合區逐漸冷 凝形成液體，而冷凝現象的產生造成溫度有下降的趨勢。

圖 6 顯示了壓力的變化情形，當震波傳遞到在右半邊液-氣二相 混合區而擴散波傳遞進入左半邊時的壓力圖，由圖中可見壓力由高壓 區逐漸往低壓區推進，而冷凝現象的出現造成交界面之壓力開始降 低。圖 7 顯示了密度的變化情形，由於震波與擴散波的影響，交界面 的 密 度 分 佈 也 因 為 冷 凝 現 象 的 產 生 而 逐 漸 逼 近 密 度 値 為 1050

/ m3

kg 。圖 8 顯示氣體體積佔有率的變化情形。當震波通過液-氣二相 混合區後，混合區內的氣體凝結成液態水。氣體體積佔有率也從

5 .

=0

αg 開始，在震波通過後達到α_g =0。圖 9 描述速度分佈的情形。

可以觀察到震波擴散波的傳遞速度在純液體和液-氣二相混合區裡是 不相同的，而且由於左右二半邊因為壓力差的關係，高壓區往低壓區 移動造成純液體往液-氣二相混合區逐漸推進，使得速度由初始值 0

s

m/ 開始增加，最終達到約 250 m/s左右。在交界面的壓力波和不連 續處的解析，則由於使用 AUSMD 含 MUSCL 類型 minmod limiting interpolation 為主要之通量分離法的結果，而沒有看見虛假的數值震 盪誤差。

圖 6 在時間 150μs 時壓力與位置的描述圖

圖 7 在時間 150μs 時密度與位置的描述圖

圖 8 在時間 150μs 時蒸氣體積佔有率與位置的描述圖

圖 9 在時間 150μs 時速度與位置的描述圖

3.2 二維注油器

當流場截面積產生變化時（如突縮管、突張管等），流體為了保 持質量、動量、通量守恆便會產生變化，如壓力與密度變化、產生回 流區等現象。第二個例子將以二維注油器（假設為圓柱型、孔口轉角 處具有尖銳邊緣）來模擬空蝕現象。圖 10 和圖 11 分別顯示注油器數 值模擬所使用之幾何外形尺寸比例結構圖和網格外型圖，總網格數為 150 X 100，而細管長度與細管直徑比則為 10：1，粗管與細管之管 徑比為 8：1。

為了可以簡單預測空蝕現象發生的區堿與範圍，假設二維注油器 的入口速度為 10 m /s、液態二氧化碳流體密度為 1200 kg/ m³和氣體 體積佔有率為 0 %、初始溫度為 300K、入口壓力為 10 MPa、壁面 為無滑移邊界速度為 0 m/s，為了使數值計算穩定，因此將細管流道 長度加長，由於速度會因流體黏滯性的影響，導致細管中心流速最 大，越接近管壁，則流速會漸漸降低直到為 0 m/s。

當流體進入細段處時，在轉角處會開始產生局部的低壓區，當流 體壓力降至飽和蒸氣壓以下時便產生空蝕現象。以質量、動量守恆之

觀點來說，當流體進入細段處時，造成流場截面積縮減與回流區的產 生，導致細管中心流速、動壓增加和靜壓減少，若此時壓力低於流體 的蒸氣壓時，空蝕現象就有可能發生。

在注油器中，隨著時間步進流體突然由粗管尖銳處進入管徑大幅 度縮小的細管前端時，由於幾何外型上的尖銳轉角處，會迫使連續的 流體分離成兩部份，而且在入口處產生維納收縮現象(Vana Contracta) 如圖 12，同時在轉角處開始產生回流區如圖 13 所示，所謂的維納收 縮現象為流體流過具有尖銳轉角的細管前端時，其流體直徑會小於細 管直徑，這種流體直徑收縮的現象稱之。維納收縮現象和孔口間的距 離與流體流動的狀態（雷諾數）、孔口直徑與流場截面積的比值有關。

由文獻[38]中可得知回流區的大小取決於幾何外型的管徑比以及流 體速度的大小。當雷諾數越大時表示流體速度越快，因而回流區的長 度也會相對增長，但是當流場形成紊流時，主回流區的長度將維持定 值。

圖 10 注油器之幾何外型比例結構圖

圖 11 注油器之計算網格

圖 12 維納收縮示意圖

圖 13 注油器流場局部放大之流線分佈圖

圖 14 到圖 28 分別顯示在初始溫度為 300 K、時間分別為 0.78、

0.98、1.28、1.755、2.00 S的氣體之體積佔有率、壓力與速度分佈 示意圖，圖 29 為注油器在無因次化之中心線出口速度-入出口壓力差 與 Chaves[39]實驗值比較圖。

從圖 14 到圖 16 可知，當流體在 0.78 S時從粗管進入細管時，

可以看到細管入口處的壓力和速度變化劇烈，而空蝕現象也在尖銳轉 角處出現，此時空蝕區域的壓力低於蒸氣壓，速度則是因為流體流入 尖銳轉角處，造成流場不穩定迫使流體結構重新調整並分離成兩部份 產生維納收縮現象，導致流體由於黏滯力的關係在靠近細管壁面時速 度較低、中心處速度較高。從圖 17 到圖 19 可知，0.98 S時因為回 流區開始出現引起流體與細管壁面發生較大的分離現象，結果造成流 體能量的損失，並導致流場截面積縮減、流場中心速度增加，空蝕區 和低壓區範圍開始漸漸增加，當離開維納收縮區後，因為壓力回升至 蒸氣壓以上，造成部分流體由氣相變為液相，根據 Bernoulli 定律可 知，當流場截面積縮小導致流速加快，則壓力會相對降低，因此流體 進入注油器細管而導致速度加快、壓力降低的現象，符合 Bernoulli 定律。圖 14 到圖 28 中所示的壓力與速度圖亦符合此一現象。從圖

20 到圖 22 可知，1.28 S回流區增至最大範圍約從 0.4 mm 成長至 1.2 mm，氣體體積佔有率在 90％以上佔大多數、70％以下數量開始逐漸 增多並且氣泡會隨著時間步進而漸漸往出口移動，低壓區範圍約從 1.2 mm 成長至 3 mm，從圖 23 到圖 25 可知，1.755 S為回流區消失 的前ㄧ刻，此時氣體體積佔有率在 90％以上維持在約 2.8 mm，70％

以下的範圍劇增，低壓區範圍幾乎佔滿細管的前半段，從圖 26 到圖 28 可知，在 2.00 S時細管壁面佔滿了氣泡而往出口移動的氣泡最終 會消失，低壓區的範圍幾乎維持不變，回流區消失。由這些圖可得知，

第一點為空蝕現象集中的區域除了維納收縮區之外，大部分集中在流 體從粗管進入細管的尖銳轉角處和細管的管路壁面；第二點為低壓區 的範圍幾乎與產生空蝕現象的區域同步；第三點為回流區的形成對空 蝕區與低壓區範圍有ㄧ定的影響力，而且速度在流場中心時會最快。

圖 29 為 Chaves 利用定義噴嘴入出口在不同壓力下之實驗，並且在空 蝕現象出現後量測出口速度之無因次化後數據比較圖，由圖中可看出 數 值 解 的 趨 勢 近 似 實 驗 值 而 速 度 略 有 些 偏 低 ， 這 可 能 是 因 為 Sanchez-Lacombe 狀態方程式無法真實模擬相變化過程的熱力性質 導致數值計算上產生誤差。文獻[40]拍攝空蝕現象的照片與本文所使 用之數值方法與狀態方程式，所得之數值解與物理現象相似。

圖 14 注油器在時間 0.78S之體積佔有率分佈放大圖

圖 15 注油器在時間 0.78S之壓力分佈放大圖

圖 16 注油器在時間 0.78S之速度分佈圖

圖 17 注油器在時間 0.98S之體積佔有率分佈圖

圖 18 注油器在時間 0.98S之壓力分佈圖

圖 19 注油器在時間 0.98S之速度分佈放大圖

圖 20 注油器在時間 1.28S之體積佔有率分佈圖

圖 21 注油器在時間 1.28S之壓力分佈圖

圖 22 注油器在時間 1.28S之速度分佈圖

圖 23 注油器在時間 1.755S之體積佔有率分佈圖

圖 24 注油器在時間 1.755S之壓力分佈圖

圖 25 注油器在時間 1.755S之速度分佈圖

圖 26 注油器在時間 2.00S之體積佔有率分佈圖

圖 27 注油器在時間 2.00S之壓力分佈圖

圖 28 注油器在時間 2.00S之速度分佈圖

為了使流場達到完全發展流，因此將細管流道長度加長為細管直 徑的 25 倍，網格數為 225 X 100 個，如圖 30 所示。當注油器下游部 分的速度向量分布達到ㄧ致時即為完全發展流，如圖 31 所示。圖 32 為細管前半段的速度流線分布圖，可以看到回流區形成的範圍，由於 細管長度拉長關係使得回流區範圍變大。圖 33 則為此時的氣體體積 佔有率局部分佈圖，在圖中可以看到整個細管幾乎佈滿了氣泡。圖 34 為氣體體積佔有率局部放大分佈圖，從圖中可發現空蝕現象氣化 較劇烈的部分主要集中維納收縮區和細管壁面處。

圖 30 注油器之計算網格

圖 31 注油器細管後半段速度向量局部分布圖

圖 33 注油器體積佔有率局部分佈圖

第四章 總結

在本文所發展的數值模式是利用 Sanchez-Lacombe 狀態方程式 描述液相、混合相、氣相的熱力性質，藉此模擬流場性質的分布情形，

分析空蝕現象所發生的位置與時機。根據先前所做模擬之過程與結 果，分為下列幾點作總結：

1.在第ㄧ個例子ㄧ維震波管中，在交界面的壓力波和不連續處的解 析，則由於使用 AUSMD 含 MUSCL 類型 minmod limiting interpolation 為主要之通量分離法的結果，而沒有看見虛假的數值震盪誤差。

2.在第二個例子二維注油器中，當流體進入窄化的流道時，根據 Bernoulli 定律會導致流場速度加快和壓力降低，當壓力低於蒸氣壓 時，相變化與空蝕現象也會相應產生。而空蝕現象出現的區域大部分 集中在流體從粗管進入細管的尖銳轉角處、維納收縮區與細管的管路 壁面。

3.將第二個例子延伸，把二維注油器細管的長度與半徑比由原來的

乎佈滿整個細管。

4.Sanchez-Lacombe 狀態方程式配合液-氣平衡公式，用以預測真實 流體液、氣和兩相混合區的性質，搭配 AUSMD 數值方法所模擬出兩 種的例子其結果大致上符合物理現象。

5.Sanchez-Lacombe 狀態方程式是一個非線性方程式，因此在求解流 體的熱力性質時是不穩定，所以本文模擬成功的算例很少，在尋求更 適當的狀態方程式上有待努力。

6.未來期望能將此數值方法做進一步的改善並延伸到更多維的數值 應用上，以及將表面張力、阻力項等因素加入數值的模擬，使模擬能 更接近真實的流體特性。

參考文獻

[1] 林於金, 雙方程雙相流模式之數值計算, 中華大學機械與航太工程研究所 碩士論文, 1996.

[2] T. Andrews, On the continuity of the gaseous and liquid states of matter, Philosophical Transactions of the Royal Society of London, Vol. 159, (1869), pp.

575-590.

[3] Gordon J. Van Wylen Richard E. Sonntag Claus Borgnakke, Fundamentals of Classical Thermodynamics, Fourth Edition 1995 pp.463-468.

[4] R. T. Knapp, J. W. Daily, and F. G. Hammitt, “Cavitation,” McGraw-Hill Book Company, 1970.

[5] P.L. Roe, The use of Riemann problem in finite difference schemes, Lect. Notes Phys. 141 (1980) 354–359.

[6] A. Harten, P.D. Lax, B. Van Leer, On upstream differencing and Godunov schemes for hyperbolic conservation laws, SIAM Rev.25 (1981) 35–61.

[7] J. E. Romate, An approximate Riemann solver for a two-phase flow model with numerically given slip relation, Comput. Fluids, 27 (1998) 455–477.

[8] Toumi, I., “A Weak Formulation of Roe’s Approximate Riemann Solver,”

Journal of Computational Physics, Vol. 102, No. 10, 1992, pp. 360–373.

[9] Tiselj, S. Petelin, Modelling of two-phase flow with second-order accurate sche me, Journal of Computational Physics, 136 (1997) 503–521.

[10] J. M. Masella, I. Faille, T. Gallouet, On an approximate Godunov scheme, Int. J.

Comput. Fluid. Dyn. 12 (1999) 133–149.

[11] Yang-Yao Niu, ‘Simple Conservative Flux Splitting For Multi-Component Flow Calculations‘, Numerical Heat Transfer, Part B, 38:203-222, 2000

[12] J. R. Edwards, R. K. Franklin, M. -S. Liou, Low-diffusion flux-splitting methods for real fluid flows with phase transition, AIAA J.38 (2000) 1624–1633.

[13] M. -S. Liou, and C. J. Steffen, Jr., ‘ A New Flux Splitting Scheme ’ J. Comput.

Phys. 107, 23 ,1993

[14] M. -S. Liou, ‘ On a New Class of Flux Splittings ’, Lecture Notes in Physics 414, 115, 1993

[15] S. Evje, K.K. Fjelde, Hybrid flux-splitting schemes for a two-phase flow model, Journal of Computational Physics,175 (2002) 674–701.

[16] S. Evje, K.K. Fjelde, On a rough AUSM scheme for a one-dimensional two-

phase flow model, Comput. Fluids 32 (2003) 1497–1530.

[17] Wada, Y., and Liou, M.-S., “A Flux Splitting Scheme with High-Resolution and Robustness for Discontinuities,” AIAA Paper 94-0083, Jan. 1994.

[18] Chih-Hao Chang and M.-S. Liou,’ A New Approach to the Simulation of Compressible Multifluid Flows with AUSM+ Scheme’, AIAA 2003-4107

[19] Yabe, T., Xiao, F., and Zhang, Y., “Strategy for Uni. ed Solution of Solid, Liquid, Gas, and Plasmas,” AIAA Paper 99-3509, June 1999.

[20] Hirsch, C., Finite Volume Method and Conservative Discretizations, Numerical Computation of Internal and External Flows, Vol. 1, Chap. 5, Wiley, Chichester, pp. 201-236., 1988.

[21] Jameson, A., Schmidt, W., and Turkel, E., Numerical simulation of Euler equations by finite volume methods using Rung-Kutta time stepping schemes, AIAA Paper 81-1259, AIAA 5th Computational Fluid Dynamics Conference.

1981.

[22] Jameson, A. Sucesses and challenges in computational aerodynamics, AIAA Paper 87-1184,Proc. AIAA 8th Computational Fluid Dynamics Conference.1987.

[23] Yee, H. C., A Class of High-Resolution Explicit and Implicit Shock -Capturing Methods, NASA TM-101088, February 1989.

[24] Nurick, w. H.,“Orifice Cavitation and Its Effect on Spray Mixing”,Fluids Engineering, Vol. 98, no. 4, 1976, pp.681-687.

[25] H. Hiroyasu, M Arai, and M. Shimizu, “Break-up Length of a Liquid Jet Internal Flow in a Nozzle,” ICLASS-91 Gaithersburg, MD, U.S.A., pp.

275-282, 1991.

[26] Koivula, T., “On Cavitation in Fluid Power, ” Proc. of FPNI-PHD Symp.

Hamburg , Vol.1, pp. 371-382, 2000

[27] David P. Schmidt, Christopher J. Rutland and M. L. Corradini, “A Numerical Study of Cavitating Flow Through Various Nozzle Shapes,” Engine Research Center, University of Wisconsin, Madison. 971597.

[28] David P. Schmidt, M. L. Corradini, “Analytical Prediction of the Exit Flow of Cavitating Orifices,” Engine Research Center, 1500 Engineering Dr. University of Wisconsin, Madison. WI 53706.

[29] David P. Schmidt, M. L. Corradini, “One-Dimensional Analysis of Cavitating Orifices,” Engine Research Center, 1500 Engineering Dr. University of Wisconsin, Madison. WI 53706.

[30] David P. Schmidt, Tzay-Fu Su, Kayhan H. Goney, P. V. Farrel, M. L. Corradini, “Detection of Cavitation in Fuel Injector Nozzles,” Engine Research Center, University of Wisconsin, Madison. WI 53706.

[31] Richard Saurel and Jean Pierre Cocchi,’Numerical Study of Cavitation in the

Wake of a Hypervelocity Underwater Projectile.’, J. of Propulsion And Power.

Vol.15, No.4, July-August 1999.

[32] Arcoumanis, C., Gavaises, M., and French, B, ‘Effect of Fuel Injection Processes on the Structure of Diesel Sprays,’ SAE Transactions, Vol. 103, No. 3, pp. 1025 -1064, 1997

[33] Yang-Yao Niu and Wei-Yuan Hsu, “A Preliminary Calculation of Three- Dimensional Unsteady Underwater Cavitating Flows near Incompressible Limit “, Communications in Computational Physics, Vol. 4, No. 4, pp. 894- 910, OCT, 2008 (SCI)

[34] Yang-Yao Niu ,Yee-Ming Lin and Yu-Chin Lin ,”A Simple and Robust Advection Upwind Flux Splitting to Simulate Transient Cavitated Water-Vapor Flows”, Numerical Heat Transfer, Vol. 51(7), pp. 679-696, 2007 (SCI)

[35] Sanchez, I. C., and Lacombe, R. H., “An ElementaryMolecular Theory of Classical Fluids: Pure Fluids,” Journal of Physical Chemistry, Vol. 80, No. 21, 1976, pp. 2352–2362.

[36] F. Harlow, A. Amsden, Fluid dynamics, Technical Report LA-4700, Las Alamos National Laboratory (1971).

[37] Cole, R. H., Underwater Explosions, Princetion Univ. Press, Princeton, NJ, 1948 [38] 洪健元 , 內流場空蝕現象之初步數值, 中華大學機械與航太工程研究所

碩士論文, 1991.

[39] H. Chaves et al., Experimental Study of Cavitation in the Nozzle Hole of Diesel Injectors Using Transparent Nozzles, SAE Paper 950290, pp. 199-211, 1995.

[40] J.Philip Drummond, M. Yousuff Hussaini, and Thomas A. Zang, “Spectral Methods for Modeling Supersonic Chemically Reacting Flowfields”,AIAA J., Vol.24, No.9, pp.1461~1467 , 1986