行政院國家科學委員會專題研究計畫成果報告

(1)

行政院國家科學委員會專題研究計畫成果報告

多相流場數值模擬技術開發(I) 研究成果報告(完整版)

計畫類別：個別型

計畫編號： NSC 95-2623-7-216-002-D

執行期間： 95 年 01 月 01 日至 95 年 12 月 31 日執行單位：中華大學機械工程學系

計畫主持人：牛仰堯共同主持人：賴明治

計畫參與人員：碩士班研究生-兼任助理：徐為淵、邱詠憲、莊詠成大學生-兼任助理：王凱弘、許廷旭

處理方式：本計畫涉及專利或其他智慧財產權，2 年後可公開查詢

中華民國 96 年 03 月 20 日

(2)

行政院國家科學委員會補助專題研究計畫期末報告

※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※

※ ※

※

多相流場數值模擬技術開發(I)

※

※ ※

※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※

計畫類別：█個別型計畫 □整合型計畫

計畫編號：NSC 95-2623-7-216-002-D

執行期間：95 年 01 月 01 日至 95 年 12 月 31 日

計畫主持人：牛仰堯

執行單位：中華大學

中華民國 95 年 12 月 31 日

(3)

摘要

在 AUSM 算則中，壓力-速度耦合(coupling)的目的是為了要修改通量的分離 (flux splitting)來達到空蝕氣泡流(cavitated bubbly flows)的準確模擬。如我們所知，

在多相流界面的解析中，壓力波在根據可壓縮的尤拉算則上仍是一個具有挑戰性的主題。現行的研究是計算在AUSM 算則中的數值消散項(diffusive terms)以及提出的一個強健方式來處理壓力-速度在 AUSM+上的耦合項。AUSM+修正的提出是在提高相界面的解析精度、非穩態震波管問題中的液體壓力波和與經證實的數據比較後的外流場空蝕流的計算結果。本計畫求取較為精確且廣泛的預測液-氣二相流體之相變化情形，研究發展混合守恆之二相流體相變化數值模式。以混合三種狀態方程式之近似守恆 Navier-Stokes 方程式做為統御方程，有別於傳統以二個統御方程式分別求取二種不同相位，本研究發展之混合模式配合體積佔有率及狀態方程式較能精確的模擬液-氣二相所佔之體積比，並拓展至三維二相流體之分析。

最後再以迎風算則之MUSCL 型式 AUSMD/V 法處理空間離散；具有 TVD 效應之三階Runge-Kutta 處理時間離散，使程式能在數值穩定與計算時間上求得較佳之平衡。

Keyword: Cavitation, Multi-Phase Flow, Upwind Splitting, AUSM

(4)

目錄

目錄 ……… II 符號說明 ……… III

一、緒論 ……… 5

1.1 前言 ……… 5

1.2 文獻回顧 ……… 5

二、數值模式 ……… 9

2.1 統御方程式 ……… 9

2.2 狀態方程式 ……… 11

2.3 數值方法 ……… 15

三、數值結果 ……… 19

3.1 一維震波管現象解析 ……… 19

3.2 二維水翼模擬……… 22

3.3 二維噴嘴模擬……… 23

3.4 二維鈍頭體模擬 ……… 23

四、結果討論 ……… 24

五、參考文獻 ……… 25

(5)

符號說明

a

流場混合聲速

αg

氣體相之積體比

αl

液體相之積體比

ρ

流體混合密度

ρg

氣體相之密度

ρl

液體相之密度

ρlsat

液體保持純相狀態之密度下限

E

流體總能量

e

流體總內能

eg

氣體相之內能

el

液體相之內能

(6)

一、緒論

1.1 前言

在真實流場中，流場往往以多相(Multi-Phase)的形態存在。常見如分離流 (Separated Flows)、氣泡流(Bubbly Flows)等。而在工業用管路系統中，多相流場則常出現於燃燒鍋爐、工廠排放管、火焰噴嘴等。當管流中液體受熱，或流道改變因而造成液體急速氣化時，稱之為空蝕(Cavitation)現象。

在計算可壓縮二相流體時有其困難之處，首先氣體的壓縮性要比液體來的強，事實上液體更為接近不可壓縮流。再者，空氣中之聲速與液體之聲速比約差四倍之多，這表示以傳統 Navier-Stokes 求解二相流體時必需解決極為複雜之數值方法。包含氣、液、混合三種狀態，數值方法中必需同時處理完全不同之聲速、

壓縮性、密度、能量、壓力等等。況且在流體力學中目前並無完美之相位轉換公式，這表示相位轉化時一定有相關性之數值震盪存在。二相流體在數值模擬上難度極高，為求降低相位轉換時發生的不正常震盪，每年都有無數學者發展出新的數值方法。在本文中將介紹求解超高速可壓縮流之數值方法，並預先以一維程式驗證其正確性。

1.2 文獻回顧

用來模擬完美氣流問題的高解析算則的發展在過去20 年間已經達到極大的進展。根據震波以及在單相流(single phase flows)中不連續的成功解析，這些成功的數值技巧都確實的應用在可壓縮多重流體(multifluid)和含有複雜相界面、材質不連續以及相轉變行為的多相流(multiphase flows)精確模擬。這些與上面提到多相轉變流類型的工程問題的應用都廣泛的在動力冷卻系統、燃料傳輸系統、水下空蝕現象以及液滴腐蝕問題中見到。如[1,2]所提到的，空蝕是一個在多相流中重要的現象。最為人所知的空蝕現象關係到三個觀點:由於流體內部的壓力低於飽和蒸汽壓導致流體內部氣泡的形成、成長和崩潰。

在大部份討論空蝕現象之文獻中，多偏向於研究低馬赫數與不可壓縮流

(7)

[3-5]。Delannoy.& Kueny.1990.[6]利用 TVD 法求解穩態 Barotropic 方程式，預測無黏性，不可壓縮流之空蝕現象，其結果大致符合實驗結果。

近幾年，Schmidt.et.[7] 發展出計算可壓縮流之空蝕數值模式，其程式中也另加一三階震波捕捉，以消除密度在相變化時所產生之擾動。Richard.和 Jean.1999[8]

則發展一計算高速可壓縮流(Very-High-Velocity Compressible Flow)之空蝕數值模式，由改變保守形式之尤拉方程，以達到消除相變化時所產生之數值誤差，並取得良好之成果。

蒸汽泡或是流體內的空蝕常常是以很高的頻率在崩壞並且產生動力壓迫去侵蝕微小的金屬元件壁面。因此，對空蝕現象有所認識並避免金屬壁面的侵蝕是必要的。在空蝕流的模擬中最主要的問題是由於多重聲音的現象[9,10]導致多尺寸 (multi-scale)時間存在於混合流中。在兩相混合中的音速和在個別的相中的音速相比可能會非常的低。因此，多相流具有頻繁的局部特徵，這個特徵在流體內伴隨著震波的存在可能會是次音速或甚至是超音速的，儘管流體的膨脹可能維持本質上的不可壓縮性。這些狀況顯示了一個巨大的挑戰，特別是在算則的數值穩定度上。目前傳統之氣-液二相數值模式往往無法有效處理高速可壓縮流在低壓區型成後所造成空蝕現象[11]，如水面下超高速投射物之空蝕現象[12]，在國防工業中高速魚雷、潛艇等都屬於此類之空蝕現象產生。當然發展此數值模式並非只能用在水面下高速投射物，類似的結構還有高速噴嘴。有些高壓注油器會以極高的速度將流體注入燃燒室內，在此類突縮流道中雖然壓力極高，然而一旦流體進入窄化流道後其產生之空蝕現象卻占了非常重要的地位。這種突然的氣化現象會引發非常嚴重的不良後果[13]，如何計算此類空蝕現象已成了流體力學中重要的課題。

為了求解保守形式之 Navier-Stokes 方程式，在空間離散上通常使用迎風

（upwind）法，細分又可分成通量差分分解法(FDS)與通量向量分解法(FVS)兩種。

Roe[14]與 Osher 等人[15]兩者即為 FDS 中較為著名者；而使用 FVS 中，Steger 和 Warming[16]與 Van Leer[17]兩者則最常被人引用。最近幾年，結合 FDS 法可精確

捕捉線性波的設計與具備 FVS 法之簡單本質特性，又不增加複雜性的混合

（hybrid）法被大量提出，Liou 於 1993 年發表的 AUSM[18-19]法為其中較為成熟

(8)

之一，近年來 AUSM 系列已發展出 AUSM+[20]與 AUSMDV[21]。在求解可壓縮多相流場時，由於相位轉換時任何狀態方程式都無法真正達到完美境界，造成 Navier-Stokes 方程式已非守恆形式，而且計算 eigensystem 有其困難之處，因此使用Roe 或者 Godunov 法變成極為麻煩。而使用 AUSM 來求解多相流場已被証實為非常有效率之方法，Chang[22]與 Niu[23]等使用 AUSM 法求解多相流場均取得不錯的結果。

通常統御方程式中時間項的積分必須以疊代(iteration)的方式隨時間趨進以接近解析解(analytical)，也就是所謂的時間趨近法[24]。Jameson et al.[25-26]發展出使用 Runge-Kutta 之多重步階法以處理通量平衡，常用於處理 Navier-Stokes 之 Runge-Kutta 法，可細分為二階、三階、四階。使用高精度之 Runge-Kutta 法雖然較為精準，但相對便耗時許多，於求解多相流場時更為明顯。為平衡計算時間與精度，本研究中將使用具TVD[27]效應之三階 Rung-Kutta 法。

另外，在求解高速二相流體時必需注意到的是，液-氣二相都必需是可壓縮性的，並且二者在計算過程中其純相狀態能獨立存在而互不干擾。首先考慮可用來計算數值模式則是液- 氣二相保守型式方程式 (two phase formulation of the conservation equation)，此傳統數值模式在二相中皆為可壓縮性，且有現成的文獻 [28-29]可參考。然而處理高速空蝕現象時，這類數值模式存在一些已知的缺點，

其一為在純氣、純液之狀態下此方程式並無確切的數學根據，其二則為無法精確處理二相界面，事實上根本無法有效計算液、氣之體積比，Rogue 等人[30]在 1998 年就曾明確的指出其缺點。基於上述原因，文獻[8]中改以混合模式求解高速空蝕現象。

多種數值方法在過去已經被提出來以精確的捕捉流體的界面，例如VOF法 [31,32]，等位函數法(level set method)[33-36]，向前追蹤法(the front tracking method)[37]以及界面捕捉法[38-41]。然而，可壓縮性的效果和相的轉變在他們的數值模式中都不予考慮。如同我們所知的，在低速空蝕流的數值計算中，低馬赫速液態流體通常都被認作是可壓縮的。然而，在液體壓力和密度場之間的耦合是非常弱的，特別是在接近不可壓縮的方式上。即使在流體內加入非常大的壓力梯

(9)

度，在密度上的變化也是非常微小的。

相對來說，在密度場中微小的數值誤差可能會導致壓力波估計上的劇烈跳動，以及相界面的不連續這些都是在液-氣空蝕流中重要的現象。此外，關於高速水下空蝕模式的模擬，Saurel和Cocchi [39] 不只考慮了壓縮性的效果，也考慮跟時間有關的飽和模式。並且，氣相是用理想氣體來做模擬，而液相是以Tait的狀態方程式來處理。他們使用含有人工消散項的高階中央有限體積的離散來計算組合過的狀態方程的混合真實流體模式去模擬相的轉變或是單相流的現象。傳統的守恆算則可能會造成在估計的擴散波裡溫度的不正確增加。為了要預防這種情形，

Saurel和Cocchi [39] 使用了一個非守恆形式的能量方程以及一個震波檢驗器來把算則回復到接近震波的守恆形式。這個求解技術的使用不易受到擴散波溫度的不正確提高而影響。然而，能量的守恆特性在次音速流裡將不再成立。為了克服Saurel 和Cocchi [39]方法的缺點，我們在此提出修正之前在多流體問題上使用的守恆上風分離法的研究[41]。

在我們先前的研究中[41]，在空氣對水1：1000的密度比下解析在兩相無空蝕流中的空氣-水的界面已展示過。我們建議在多流體組成的動力計算上使用加入 VOF模式的延展尤拉方程和用在界面捕捉的AUSM或AUSMDV算則為佳[42]。保守形式在模組方程式中已列入考慮。在之前的計算中，當液態流體的狀態方程式是由嚴謹的氣體方程式(或是Tammann’s模式)所表示時，AUSM將得到震盪的數值解。這種情形導致AUSM+算則無法把壓力波解析的很好。然而，AUSMDV和精確解以及在很多一維的兩相測試例子中所顯示的自由震盪物質界面相比證實了相當好的精度。不僅現行的多流體組成模式使用壓力守恆的更新來提高壓力在通過界面時的平衡，而且根據壓力和密度比的混合函數AUSMD/V也維持了在穿過不連續接觸時壓力的連續性。配合VOF模式的尤拉方程式變成一個可能的方法來達到界面和不連續材質的高解析精度。然而，AUSMDV雖能成功模擬不相融的兩相流卻還不能用在可融合流的例子上。因此，AUSMDV以及在[39-43]一個修正後的溫度相關飽和模式將在本篇壓力波的解析和相界面處來加以研究。測試的例子包含了一維的空蝕震波管、震波-凝結和二維內、外流場的暫態空蝕流。

(10)

二、數值模式

在本研究中，在空蝕流場之現象計算，在時間離散上採用具TVD 效應之三階 Runge-Kutta 法；空間離散則使用高階 MUSCL+Limiter Function 以提高精度；而在通量分離的部份，本研究以AUSMDV 作為主要之通量分離法。

2.1 統御方程式 混合模式

為了延續先前求解的混合模式結果，一個含有體積佔有率模式和黏滯項作用的二維二相混合模式用無因次的保守形式可表示為：

v v

F G

Q F G

t x y x y

⎛∂ ∂ ⎞

∂∂ +∂∂ +∂∂ =⎜⎝ ∂ + ∂ ⎟⎠ (1)

2

; = ; = 2

( ) ( )

u v

u p uv

Q u F uv G v p

v p p

E u E v

E

ρ ρ

ρ ρ ρ

ρ

ρ ρ

⎛ ⎞ ⎛ ⎞

⎛ ⎞ ⎜ + ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟ ⎜ ⎟

=⎜ ⎟ ⎜ ⎟ ⎜ + ⎟

⎜ ⎟ ⎜ + ⎟ ⎜ + ⎟

⎝ ⎠ ⎜⎝ ⎟⎠ ⎜⎝ ⎟⎠

黏滯通量向量F_ν，G 可在[36]中查到。並且，總能可以表示為： _ν

( )

∑

=

+ +

= ^m

i

i i i i

i e u v

E

1

2 2

2 1α ρ ρ

α (2) 此外，考慮一非保守體積佔有率方程為：

=0

∂ + ∂

∂

v y u x

t

α α

α

(3)

這裡α代表氣體的體積佔有率，其他的即為(1−

α

)。照例來說，

α

∈[0, 1]而ρ,u, E 分別為混合密度、混合速度和混合總能。每種流體均要滿足真實Tammann 狀態方程式：

−

∞

−

= e p

p ( γ 1 ) ρ γ

(4)

(11)

這裡每單位質量的內能表示為e，γ 表示比熱比而

p

_∞是已規定的類壓力常數，這些值可以用來描述在[33]裡說過的多個材質的性質在不同狀態下的影響。當 γ =1.4和

p

_∞ = 0時，真實氣體就變成完美氣體。此時狀態方程式的音速如下所示：

( ⁺

_∞

)

= P P

a ρ

γ

(5)

要開始計算多材質流體模組，像ρ_i , ,u_i p_i ,γ_i和

p

_i_,_∞這些原始變數，在網格上是已知的。像

ρ ρ ρ , u , v

和 E 這些混合變數是定義為對每個網格來說

ρ

_i

, u

_i 和

E

_i分之體積-重量的總合。混合速度 u 可以用以下式子求得：

∑ ∑

=

i i m i

i i i m

i

u

u u

ρ α

ρ α ρ

ρ

1

1 (6)

在計算混合壓力以前，

γ

對在壓力平衡的多流體成份的混合物的影響可以由以體

積佔有率為主的平均來得到

∑

= −

− =

m

i i

i

1 1

γ α

γ

(7)

很明顯地，當γ 為已知時，從第(5)式，未知的 p 和 p_∞可以根據從

⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟ −

⎠

⎜⎜ ⎞

⎝

⎛

= ∑ − ∑

=

m

i i

i m

i i

p p

1

1 γ 1

α γ

α

(8)

和

⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ + −

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

∑

−

∑

=

∞ ∞

m

i i

i m

i i

i

p

1 1

,

1 1

1

γ

α γ

γ

α

(9)

更新的γ 來得到。

此外，當網格是由同樣的γ_i 不同的 p 所組成時，混合壓力就可以由_i

(12)

1 m

i i i

p=

∑

₌α p 來求得。這也說明了如果流體的構成要素是在壓力平衡下，那麼 p 和u 在交界面上就是連續的。

由統御方程式(1)所得到的Jacobian matrix 為:

2 2 2

0 0 0

( 1) ( 2) ( 1) ( 1) 0

( 1) ( 1) ( 2) ( 1) 0

( 1) ( 1) ( 1) 0

( 1) ( 1) ( 1)

0

x y

x x y x x

y x y y y

x y

uU b U u u v

vU b v u U v

A a a a

U b Ub b Uu b Uv U

UY U

ξ ξ

ξ γ ξ γ ξ ξ γ ξ γ

ξ γ ξ ξ γ ξ γ ξ γ

γ ξ γ ξ γ γ

γ γ γ

ξ α ξ α

⎛ ⎞

⎜ − + − − − − − − ⎟

⎜ ⎟

⎜ − + − − − − − − ⎟

⎜ ⎟

=⎜⎜⎜⎜⎝− ⎛⎜⎝ − −+ + −⎞⎟⎠ ⎛⎜⎝ − + − −⎞⎟⎠ ⎛⎜⎝ − + − −⎞⎟⎠ ⎟⎟⎠

⎟⎟

在此 (p p ) ² γ a

ρ ^∞

+ = ;

2 2

( )

2 u v

+ b

=

其特徵值和特徵向量分別為:

1

2

2 2

3

2 2

4

5

x y

U U U a U a

U

ξ

λ λ

λ ξ ξ

λ

=

= + +

= − +

=

2 2 2 2

2 2 2 2 2 2 2 2

2 2 2 2

1 0 1 1 0

0

1 0

2 ( 1) 2 ( 1) 2 0

0 1

y x x

x x y x y

y y

x y x y

x x

x x y x y

a a

u u u

a a

v v v

X

v u

u v a u v Ua a u v Ua

ξ ξ ξ

ξ ξ ξ ξ ξ

ξ ξ

ξ ξ ξ ξ

ξ ξ

ξ γ ξ ξ γ ξ ξ

α α α

⎛ ⎞

⎜ ⎟

⎜ − + − ⎟

⎜ + + ⎟

⎜ ⎟

+ −

⎜ ⎟

=⎜ + + ⎟

⎜ ⎟

−

+ + +

⎜ + + + − ⎟

⎜ − + − + ⎟

⎜ ⎟

⎝ ⎠

(13)

2.2 狀態方程式

狀態方程式的選擇對於空蝕現象有著極重大的影響，本研究中所使用之狀態方程式，也關注於尾隨在水面下高速前進之物體所產生的空蝕現象。針對此結構，

我們必需特別注意尾隨在水面下投射物之後所產生的低壓空蝕區域及其熱力學行為。當流體之速度變的非常快時，液體的壓縮性則相對顯的非常重要，並且會引發強烈之震波(大於 50 Kbar)，另一方面，投射物後方形成的氣化現象所造成之極低壓區則是可以預見的。在純氣相與純液相之間產生如此劇烈之壓差，其結果將會導致高壓區與低壓區之間產生一液氣共存區，而在此共存區之中，液體與氣體將會是可壓縮性的。圖(2.1)顯示本研究中所使用之狀態方程式之壓度-溫度對應圖。

圖(2.1) 狀態方程式所得之壓力-溫度對應

A. 液體狀態

為求解可壓縮流，本文採用Cole.等人在 1948[44]所發展之壓力方程式，其表示如下：

₀ 1 ( )

( )

n

sat l sat

P K P T

T ρ ρ

⎡⎛ ⎞ ⎤

⎢ ⎥

= ⎜ ⎟ − +

⎢⎝ ⎠ ⎥

⎣ ⎦

(10)

(14)

其中：n = 7 ( For Water)，K = ₀ 3×10⁸ Pa

根據IAPWS-IF97 手冊[45]：

Specific internal energy ： =τγ_τ −πγ_π RT

P T u( , )

(11) Specific entropy ： =τγ_τ

RT P T h( , )

(12)

Speed of sound ：

ππ ττ

ππ π

π

γ γ τ

τγ γ

γ

− −

=

2 2 2 2

) (

) , ( RT

P T

a (13)

R = specific gas constant (0.461526 KJ/Kg．K)

π γ

τ 、、

.. 等請參考文獻[45]。

B. 氣體狀態

水蒸汽之壓力可用理想氣體方程式求取：

RT

P = ρ

(14) IAPWS-IF97 手冊則定義其它之數值：

Specific internal energy ： ( , ) ( ₀ ) ( ₀ _τ)

π τ π

τ

τ γ π γ γ

γ

τ + − +

RT = P T

u (15)

Specific enthalpy ： ( , ) ( ₀ _τ)

τ

τ γ

γ

τ +

RT = P T

h (16)

Speed of sound：

) (

) 1

) ( 1

(

) ( 2

1 )

, (

2

2 2

2 2 2

τττ ππτ

πττ πτ

ππτ

πτ πτ

γ γ τ

τπγ γ πγ

π

γ π πγ

+

− + +

−

+

= + RT

P T

a (17)

C. 飽和狀態

在混合狀態下，氣–液二相之壓力、溫度均應相等：P_l =P_g =P_Sat，T_l =T_g =T_Sat；混合壓力則由IAPWS-IF97 手冊中求取。

(15)

4

5 . 0

2 4 )

( 2 1

)

( ⎥

⎦

⎢ ⎤

⎣

⎡

− +

= −

AC B

B

C MPa

T P_Sat

(18)

其中： A=θ² +n₁θ +n₂；B=n₃θ² +n₄θ +n₅；

c

=

θ

² +

n

₁

θ

+

n

₂

10 9

1

1 n

K T

n K

T

− +

θ =

(19)

1 2 5 16 43 110

sat 3 3 3 3 3 3

1 2 3 4 5 6

( ) 1

l c

T b b b b b b

ρ θ θ θ θ θ θ

ρ ^{= +} ⁺ ⁺ ⁺ ⁺ ⁺

(20)

2 4 8 18 37 71

sat 6 6 6 6 6 6

1 2 3 4 5 6

ln ^g ( )

c

T c c c c c c

ρ θ θ θ θ θ θ

ρ

⎡ ⎤

= + + + + +

⎢ ⎥

⎣ ⎦ (21) 在此

θ = − 1 T T /

_c

D. 狀態方程式求解過程

當守恆方程式的混合形式(1)都解決後，所有混合變數都可透過計算域: , ,ρ u 和 e來決定。然後在使用適當的狀態方程式來決定其他的熱力變數。求解的一開始是要先估計局部的溫度。然後用這個溫度分別來決定液體和氣體的飽和密度(20) 和(21)。如果混合密度ρ(從守恆方程式來決定)在此溫度下比液體的飽和密度高，

那麼這個流體就是純液體並且壓力可由(10)來計算，然後新的溫度由(22)得到。另外，混合密度也可能低於飽和氣體密度。在此壓力是由(14)決定，新的溫度是由(23) 而定。當混合密度在液體和氣體飽和密度之間時，體積佔有率由(24)和(25)來決定，

壓力由(26)，新的溫度由(27)來決定。如果新的溫度和一開始估算的溫度相等，那麼熱力的狀態就被唯一決定。否則就要使用最後的溫度來當一開始所估算的然後再使用疊代程序。此過程可以達到快速收斂且沒有困難性。

(16)

其中：

0 0

( ) ( )

l l l

e T

=

Cv T

−

T

+

e

(22)

在此 0 0

: : :

vl l

C e T

液體的等容比熱參考內能

參考溫度

0 0 0

( ) ( ) ( )

g g l

e T =Cv T−T +Lv T +e (23) 在此

Lv 蒸汽潛熱

:

, 1

2 , 1

g g l l

g g l l g l

E e u u

e e e

ρ α ρ α ρ

α α α α

= + = + ⋅

= + = +

(24) 在此 ρ: ^混和密度

sat( ) sat( )

g g

T

l l

T

ρ α ρ

= +

α ρ

(25)

在此

( , ) ( )

( ) ( ) ( , ) 1

l g

g l

l g

a T T

T T

T ρ ρ ρ

ρ ρ

α ρ α

= = −

−

= − ;

^sat

sat

: :

g

l

ρ ρ

氣體飽和密度液體飽和密度

1.5 3 3.5 4 7.5

sat 0 sat 1 2 3 4 5 6

ln( P / ) ( / P = T T

_c

)( a θ + a θ + a θ + a θ + a θ + a θ )

(26)

{

^sat ^lsat ⁰ ^sat ⁰

}

⁰

( ) _g _g ( ) _g _l ( ) _l ( ) _g _g ( ) ( ) / _l

e T

= ⎡⎣

α ρ T Cv

+

α ρ T Cv

⎤⎦

T T

− +

α ρ T Lv T ρ

+

e

⁽²⁷⁾

2.3 數值方法

此系統(1)是一個嚴謹的雙曲線方程式。這允許我們應用現有的雙曲線求解器像是近似黎曼解或者通量向量分離法，這些都成功的解決單相問題，去模擬多組成模式。在本研究中，我們使用 AUSM+算則配合壓力-速度耦合的消散項來做耦

(17)

合。而此系統方程式(1) 隨著時間的演進可以被一個叫預估 - 校正器 (predictor-corrector)法，配合固定的時間步階來做離散。可如下表示：

1/3

1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2

2 3 1/3 1/ 3

1/ 2 1/ 2 1/ 2 1/ 2 1/ 2 1/ 2

1 2/3 2/3

1/ 2 1/ 2 1/ 2 1/ 2

1 ( ( , ) ( , ))

2 1 ( ( , ) ( , ))

3 ( ( , )

n n n L R n L R

i i i i i i i i

n n n L R n L R

i i i i i i i i

n n n L R n

i i i i i i

Q Q t F Q Q F Q Q

x

Q Q t F Q Q F Q Q

x

Q Q t F Q Q F

x

+

+ + + − − −

+ + +

+ + + − − −

+ + +

+ + + −

= − Δ −

Δ

= − Δ −

Δ

= − Δ −

Δ ( Q

_i^L₋^{1/ 2}

, Q

_i^R₋^{1/ 2}

))

(28)

界面上的保守變數Q1^L+1/2,Q_i^R+1/2可以透過 minmod 限制子的 MUSCL 法來決定。在空間上的差分，界面的原始變數是透過 minmod 限制子達到二階或三階空間精度的 MUSCL 法來決定。本研究中，在二階或三階空間精度所計算出的解並沒有巨大的差異。所以在此只有展現三階空間精度解的結果。隨後，根據AUSM 形式算則

n

Fi₊₁_/₂的通量外差被用在離散質量、動量和能量的對流通量。此數值通量的離散將馬上在之後的部分做一個描述。

數值的通量分離

在計算多相流場時，為避免計算 eigensystem 之困難，本研究中使用 Y.Wada 和M.-S. Liou[21]所發表之 AUSMDV。此通量公式已經被他們在[41-43,46-50]中展示的現今用在解析震波的迎風算則和在兩相間高密度比的多組成問題的界面上所做的改進來證明。

AUSMDV 在求解多相流場已為主流通量分離法之一，首先定義界面上之質通量

( ρ u )

₁_/₂

= ρ

_L

u

_L⁺

+ ρ

_R

u

_R⁻ (29)

(18)

其中，L 與 R 分別為界面上之左與右。在這裡，速度分離(u ,_L⁺ u_R⁻)可如下表示：

) 2 1 2 (

4 )

( ² L L

L L

u u u

u a

a

u u +

−

⎥ +

⎦

⎢ ⎤

⎣

⎡ +

+ −

=

+ α α (30)

) 2 1 2 (

4 )

( ² R R

R R

u u u

u a

a

u u −

−

⎥ +

⎦

⎢ ⎤

⎣

⎡ −

− −

−

=

− α α (31)

其中

a

為界面上之共用音速：

( ) 2

1

2 /

1 aR aL

a = + (32)

於此，

a

_R_與

a

_L 有多種選擇。Wada 和 Loiu [54]則以壓力與密度的比來建立此二項數值：

R L

L

L p p

p

) / ( ) / (

) / ( 2

ρ ρ

α ρ

= + ，

R L

R

R p p

p

) / ( ) / (

) / ( 2

ρ ρ

α ρ

= + (33)

然後根據音速就可以定義左邊和右邊的馬赫數：

2

a

1

M

_L

= u

^L _，

2

a

1

M

_R =

u

^R (34)

然後接下來分出的馬赫數和壓力就可以用三個多項式的形式表示：

(

M ± M

)

= Μ^±

2 1

) 1

( (35)

⎪⎩

⎪⎨

⎧

Μ

−

±

= ±

Μ^± _±

,

, ) 1 (

) 1 (

2 2 2 )

4 (

M

otherwise M <1

(36) 和

⎪⎩

⎪⎨

⎧

Μ

−

±

= ±

Ρ ±

±

) 1 (

2 2 2

) 5

( (1 )

) 1 16 (

) 3 2 ( ) 1 4( 1

M

M M M

M ∓

otherwise M <1

(37)

在M 和 P 的下標數字是表示多項式的次數。在 AUSM 形式算則[51,52]裡界面的質量通量一般形式可定義為：

) (

)

( 12

12 12

12 a Lm Rm DP

u = ρ ⁺ +ρ ⁻ +

ρ (38)

(19)

在此

1 1 1

2 2 2

1( )

m^± = 2 m ± m (39) 根據方程式(38)，我們可有多個版本的質量通量分離。首先，AUSM 的原始質量通量界面形式可以在消散項D =0_p 的假設下回復，以及

) ( )

( ₍₄₎

) 4 ) ( 2(

1 AUSM M ML M MR

m ₊ = ⁺ + ⁻ (40)

AUSM+對可壓縮氣體流來說已經證實了是一個非常強健且精準的算則，然而，根據[51-53]，為了要消除由壓力在低速流下模擬的預估所引起的奇-偶退耦的錯誤，

AUSM+的質量通量是必須要用 pre-conditioning 形式的音速和壓力-速度耦合的消散項來做修正。新的質量通量在[51]是由 AUSMDV(P)來表示。如下所示：

) ( )

( ₍₁₎

) 1 ) ( 2(

1 AUSMDV M ML M MR

m = ⁺ + ⁻ (41)

和

{

}

( ) ( )

( )/( ) ( ) ( ) ( )/( )

P L R

L R

D M M M M

P P

M M M M

P P

δ δ

ρ ρ

δ δ

ρ ρ

+ −

⎡ ⎤

= ⎣ − ⎦

⎡ ⎤

×⎢⎣ − + ⎥⎦

⎡ ⎤

+⎣ + ⎦

⎡ ⎤

×⎢⎣ + + ⎥⎦

(42)

在此

(

⁽ ⁾ ⁽ ⁾

)

(M_L M₄ M_L M₁ M_L

M⁺ = ⁺ − ⁺

δ 以及 δM⁻⁽MR⁾=

(

M4⁻⁽MR⁾−M1⁻⁽MR⁾

)

在重新整理方程式(38)-(42)後，與壓力-速度耦合的 AUSMDV(P)質量通量就變成：

( ) { ( )

( ) ( ) ⁽ ⁾ ^}

1 1 (1) (1) (4) (1) (4) (1)

2 2

(4) (1) (4) (1)

L R

L R L R

u a

P P P P

P P M M M M P P

ρ ρ ρ

ρ ρ ρ ρ

+ − + + − −

+ + − −

= Μ + Μ + Μ −Μ +Μ −Μ

⎡ ⎛ ⎞⎤ ⎡ ⎛ ⎞⎤

×⎢ + ⎜ + ⎟⎥+ − − + ⎢ − ⎜ + ⎟⎥

⎝ ⎠ ⎝ ⎠

⎣ ⎦ ⎣ ⎦

(43)

(20)

然後界面的壓力通量可以寫成：

1

2 (5)( _L) _L (5)( _R) _R

p =P⁺ M p +P⁻ M p (44) 最後此界面AUSMDV(P)的數值通量可以簡單定義成：

¹₂

( )

¹

( )

¹

( )

¹₂

2 2

1 1

( ) ( ) ,

2 ^L ^R 2 ^R ^L

F = ρu ⎡⎣Φ +Φ −⎤⎦ ρu ρΦ − Φ +ρ P (45) 在此

[ ^{1, , ,} ^{u v H} ] ^{, P} [ ^{0, ,0,0} ^p ]

Φ = =

三、數值結果

3.1 一維震波管現象解析

(i)雙成分震波管

在第一個例子中，Abgrall[38]的研究中對含有兩種被界面分開的純液體的 Sod

的震波管問題已解決。我們考慮一震波管在不同的比熱比下，在左側充滿高密度的液體，而在右側則使用低密度液體。在兩邊界面的密度比都向上修正為10⁴。兩組的實驗數值資料可以表示如下：

⎟ ⎟

⎠

⎞

⎜ ⎜

⎝

⎛

=

⎟ ⎟

⎠

⎞

⎜ ⎜

⎝

⎛

∞

0 4 1

1 1 1

p p u

L

γ . ρ

and

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∞ 0

2 . 1

1 1 0001 . 0

p

R

p u

γ ρ

在整個管子中速度和壓力都假設為相等狀態。把解繪在無因次時間t=0.1 處。混合壓力、混合密度、速度以及空間佔有率的精確解在每個圖裡都表示為固體線。計算的解表現在100 個等間距的網格上。選擇 Δt/Δx=10⁻⁵來保持數值上的穩定度。使用AUSMDV 求解的混合密度、混合壓力和氣體空間佔有率、氣體速度的預估值都分別表示在圖一。圖中的比較顯示出精確解和計算解的一致性是令人滿意的。由

(21)

此可知在界面處正確的反應是記錄在密度和空間佔有率的圖表中。這也表示出即使密度比上升到1：10000，速度和壓力的平衡在此計算中仍維持著。

(ii)震波折射問題

隨著 Karni[32]和 Shyue [37]的研究，一個震波折射問題描述了一個含有空氣和氦氣的氣流被界面所分開，氣流被一入射的震波打擊然後開始運作的問題已被解決[32]。正如 Karni 所提到的，被入射的震波打中後，震波通過界面時會變的完全不穩定，如此一來任何錯誤的數值反應都會傾向於非物理力學的結果。在此例子中，在靠近材質界面的流體行為的AUSMDV 型算則的計算精度已經研究過了。除了含有單一不連續接觸的空氣/氦氣界面的問題之外，此例子也包含了一個非常狹窄的擴散波，一個傳送微弱的震波會比入射的震波要傳遞的快。同時，一個真實氣體常用來描述氦氣的狀態。這三種固定的狀態如下所示：

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∞ 0

4 . 1

5 . 1 3535 . 0

333 . 1

p p u γ ρ

，

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∞ 0.

4 . 1

0 . 1

0 0 . 1

p p u γ ρ

air, post-shock air, pre-shock

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∞ 0.25

0 . 1

0 1379 . 0

p p u γ ρ

helium, pre-shock

混合壓力、混合密度、速度和氣體空間佔有率包含精確解的結果如圖 2 所示。計

算的解表現在400 個等間距的網格上。選擇 Δt/Δx=10⁻⁵來保持數值上的穩定度。計算出的體積佔有率和解析的資料的估算相當吻合。根據相同尺寸的網格，在不連續接觸的圖中，AUSMDV 在速度和壓力的分佈圖中達到了和精確解非常一致的

(22)

解。此外，在靠近界面的非物理解都由於AUSMDV 算則的計算解完全地減低。

(iii)一維空蝕震波管

為了達到可用於計算水面下高速投射物體所引發之相變化現象，首要之事便是發展一可預測氣化波型範圍之數值模式，此數值模式可預先以一維程式測試其正確與有效性。數值解必需符合二個反方向之擴散波，如圖(2.2)所示。

圖(2.2) 空蝕模式測試條件與其擴散波示意圖

本例子[39]考慮一個管子最初在不變的壓力和密度下，充滿固定比熱的完美氣體。最初，流體佔據右半邊的管子從右邊開始移動，左半邊管子的流體從左邊開始移動。初始條件在本例子選用溫度是300K；密度 1kg/m³，左邊管子的速度是3000 m/s，反之右邊管子流體的速度就是-3000 m/s。本例子的解反映了兩個擴散波朝相對的方向移動，導致流體在中央區域的蒸發作用。

用 AUSMDV 計算多性質的結果以及它的新版本表示在圖三在 100μs 下，當兩個擴散波開始移動時，流體速度在穿過擴散波時會遞減直到中央區域到零為止，

這顯示出從初始值1000kg/m³到純蒸汽值約少於三個等級濃度大小時，越強的擴散波會導致越強的密度遞減。蒸汽在固定的溫度T=300k 下產生，這是因為在相轉變時局部平衡的假設所致。當流體變成純氣體時，擴散波的影響持續著。除此之外，

從圖三，壓力在純液體時突然下降，然後達到飽和壓力而且在靠近中央純氣相的區域維持不變。由於強擴散波的關係，所有的液體將被蒸發。穿過這些波，氣體的體積佔有率α 從管子的兩邊進入開始從 0 增加到 1。在本測試中，由 AUSMDV 計算的數值結果幾乎和舊版本所做的相等。

(23)

(iv)一維震波凝結管解析

在本例子[39]一維震波-凝結管中，流體顯示了相的轉變。在此，一個 1 公尺長的管子在左邊充滿了初始密度為1200kg/m³的高壓液體，而在右邊充滿了初始密度為500kg/m³的兩相(液-氣)混合液體。一開始，兩個流體都是在 300K 的溫度下。

兩相的混合是由液-氣各半的質量所組成。圖四和圖五分別顯示了由 AUSMDV 舊版和新版在150μs 下計算多性質的結果。這表示了當震波的增加進入到兩相混合右

邊的空間時，在中間區域扇形區的增加進入到左邊速度到達 0 的空間。在穿過震

波時，密度會從環境的500kg/m³增加到接近震波壓縮或初始1200kg/m³的密度。這也表示了氣體的體積佔有率，在震波之後兩相的混合會瞬間凝結。氣體的體積佔有率在通過震波時從α_g =0.5達到α_g =0。從速度資料得到的數值結果，值得注意的是波的增殖速度在純液體和兩相混合中是相當不同的。從圖四使用舊版的 AUSMDV 的計算結果圖來看，很明顯可以看到在接觸面上的壓力、密度、速度會有強烈的震盪錯誤。從舊的AUSMDV 所計算的體積佔有率圖來看，可以發現在純液相和混相之間出現模糊的界面。然而，所有使用新版AUSMDV 估算的結果都是令人滿意且沒有不必要的錯誤的。

3.2 二維水翼模擬

在展示本例子前，我們先定義空蝕係數：

K =1/2( ₀²)

0

U P P _v

ρ

− =

− ∞

− P P

P P _v

0 0

這邊P0、Pv、P_∞分別是總壓、飽和蒸汽壓和自由流(free-stream)壓。第二個例子符合了次音速的水流體在二維NACA 0012 水翼的測試情形。網格包含了 300 × 64 個節點，在水翼的表面和尾流可能發生空蝕的區域給予加密。初始條件假設雷諾數

(24)

為1×10⁵，溫度為300 K，以及進入固定速度 20 m/s。入口流體假設為次音速流，

且速度、溫度和外差的壓力都固定，在出口處的次音速流所固定的壓力、外差速度和溫度將在計算中被使用。圖六繪出空蝕氣泡開始在邊緣蔓延而且由於壓力低於蒸汽壓首先出現在尾流區，然後由於尾流的逆壓而向前增加到主要的邊緣，接著和主要邊緣周圍對流到下游的生長空蝕氣泡合併。在圖六擷取到水氣體積佔有率的演進。空蝕的發生只從空蝕係數為0.8 開始表示。

3.3 二維噴嘴模擬

在第三個例子中，接下來的參數在模擬中是視為變數來處理：the needle lift (H/D)、噴嘴外觀比例 (L/D)如圖四所示。這裡P_o、P_ν、P_∞分別為注射壓力、蒸汽壓力、末端壓力。由於注射壓力在我們的模擬中都屬高壓(150 bar)，因此在每個例子中的空蝕係數大約一致。數值的結果顯示噴嘴空蝕流在不同的壓降下所做的模擬，其中蒸汽-水密度比定在 1：1000、蒸汽密度為 0.0086 (kg/m³)、動黏度 9.18×10^-6 (kg/ms)、注射壓力 150 bar 以及 10 的 L/D 比。從模擬來看很明顯的在一開始的注射區，空蝕開始出現在噴嘴的入口稜角處，接著延伸到噴嘴的下游出口，然後當壓降時往回退縮並且潰散。噴嘴中氣相的體積佔有率的變化如圖五所示(T=300k、

U=10 m/s、1 Mpa 的壓降以及雷諾數 1×10⁵)。空蝕現象的出現在空蝕係數 1.252 時開始顯示。可以看到空蝕區隨壓降在一開始區域的增加而會有增大的趨勢。現在的空蝕模組也是從前的研究中所使用的，這促進了高壓降的噴嘴空蝕流結果。

3.4 二維鈍頭體模擬

最後的例子符合了由Rouse 和 McNown[55]所帶領的實驗，表現關於流體水在半球體上的流動情形。初始條件雷諾數都是

1 . 36 × 10

⁵，溫度為300 K，進入速度固定在 4.317 m/s，但是使用不同的空蝕係數(K)。網格使用 Edward 研究中[52]的

(25)

300 × 85 個節點，其他相關流體細部條件可以在[52,56-57]找到。空蝕係數對空蝕現象發生的影響已經模擬。如圖七所示，擷取到K=0.2 到 0.4 的例子中水氣的體積佔有率的變化且由密度場的計算所驗證。圖八驗證了和實驗值比對過的表面壓力分佈，也把空蝕係數的函數參數化。當空蝕係數從0.8 下降時，擴張區的壓力降到蒸汽壓的水準，導致氣相和空蝕氣泡的生成。和實驗的數值驗證也表現出隨著空蝕係數的增加數值的精度會跟著減低。由新版的AUSMDV 計算顯示當空蝕係數在 0.4 和 0.8 時和資料是相符合的，然而，空蝕係數在 0.2 和 0.3 時，更大氣泡區產生的表面壓力係數的估計並不是完全和實驗值相符，但是傾向令人滿意的方向來驗證。

四、結果討論

在本研究中，AUSMDV 的新型態算則已經展現了其在一維非穩態空蝕和震波 -凝結管計算上的強健和精確度，無論如何，修正後的 AUSMDV 在二維水翼和鈍頭體流中的空蝕現象達到了計算上的穩定。雖然紊流的影響在我們的計算中都不予考慮，但是在鈍頭體流的模擬和實驗值的數值比對上是令人滿意而且和 Ahuja 等人的研究[56-57]所顯示的表面壓力係數的估計是非常相近的。未來的研究目標需要去研讀含有多尺寸長度和時間的物理現象，例如紊流效應和在不同相間的討論等。未來開發多尺度的多相流的計算模式與技術將是刻不容緩的一個課題。

(26)

五、參考文獻

[1] Jiangong Qin,”Numerical Simulation of cavitating Flows by the Space-Time Conservation Element and Solution Element Method”, PhD Thesis, Wayne State University, Detroit, Michigan, 2000.

[2] David P. Schmidt, Christopher J. Rutland, and M. L. Corradini, “A Fully Compressible Two-Dimensional Model of Small High-speed Cavitating Nozzle”, Atomization and Sprays, Vol.9, pp255-275,1999

[3] Lemonnier, H., and Rowe, A., ‘Another Approach in Modeling Cavitating Flows’, Journal of Fluid Mechanics, Vol. 195, pp. 557-580, 1988

[4] Kubota, A., Kato, H., and Yamaguchi, H., ‘A New Modeling of Cavitating Flows:

A Numerical Study of Unsteady Cavitation on a Hydrofoil Section’ Journal of FluidMechanics, Vol. 240, pp. 59-96, 1992

[5] Kinnas, S. A., and Fine, N. F. ’A Numerical Nonlinear Analysis of the Flow Around Two and Three-Dimensional Partially Cavitating Hydrofoils’, Journal of Fluid Mechanics, Vol. 254, pp. 151-181, 1993

[6] Y.Delannoy and J.L.Kueny,’Two Phase Flow Approach in Unsteady Cavitation Modeling’,Cavitation and Multiphase Flow Forum,ASME FED, Vol. 98, pp.

153-158,1990

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

[8] 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.

[9] L. Van Wijngaarden, One-dimensional flow of liquids containing small gas bubbles, Annu. Rev. Fluid Mech.4, 369 (1972)

(27)

[10] H. B. Stewart, B. Wendroff, Two-phase flow: models and methods, Journal of Computational Physics 56 (1984) 363–409.

[11] Davies, G. A., Flanagan, O., Hillier, R., Hitchings, D., and Lord, S. J., ‘Hydraulic Shock Loading Due to Supercavitating Projectiles’ Proceedings of the 20th International Symposium on Shock Waves, edited by B. Sturtevant, J. E. Sheperd, and H. G. Hournuns, Vol. 2, World Scientific, NJ., pp. 1207-1212, 1995

[12] Thibaudier, C., and Tosello, R., ‘Interaction Jets de Charge Creuses-Eau’, Actes du 4eme Symposium International des Hautes Pressions Dynamiques, edited by AFP La Ferte, Saint Aubin, France, pp. 281-587, 1995

[13] 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

[14] Roe, P. L., ‘Approximate Riemann solvers, parameter vectors anddifference schemes,’ J. Comp. Phys., Vol.43, pp. 357-372. , 1981

[15] Osher, S. and Solomon, F., ‘Upwind difference schemes for hyperbolic systems of conservation laws,’ Mathematics and Computaters, Vol.38, pp. 339-374. , 1982 [16] Steger, J. L. and Warming, R. F., ‘Flux-vector splitting of the inviscid gas dynamics

equations with applications to finite difference methods,’ J. Comp. Phys., Vol.40, pp. 263-293. , 1981

[17] van Leer, B., ‘Flux-vector splitting for the Euler equation,’ LectureNote in Phys., Vol.170, pp. 507-512. , 1982

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

Phys.107, 23 ,1993

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

[20] M.-S. Liou, ’A Continuing Search for a Near-Perfect Numerical Flux Scheme, Part : AUSM+,’ ,NASA TM 106524, March 1994.

Ⅰ

(28)

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

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

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

[24] Hirsch, C., ‘Finite Volume Method and ConservativeDiscretizations,‘ Numerical Computation of Internal and External Flows, Vol. 1, Chap. 5, Wiley, Chichester, pp.

201-236. , 1988

[25] 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.

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

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

[28] Butler, P. B., Lambeck, M. F., and Krier, H., ‘Modeling of Shock Development and Transition to Detonation Initiated by Burning in Porous Propellant Beds’

Combusion and Flame, Vol. 46, pp. 75-93, 1982

[29] Saurel, R., Larini, M., and Loraud, J. C., ‘Ignition and Growth of a Detonation by a High Energy Plasma’ Shock Waves, Vol. 2, No. 1., pp. 91-102, 1992

[30] Rogue, X., Rodriguez, G., Haas, J. F., and Saurel, R., ‘Experimental and Numerical Investigation of the Shock-Induced Fluidization of a Particle Bed’, Shock Waves, Vol. 8, No. 1, pp. 29-45, 1998

[31] W. Hirt, B. D. Nichols, “Volume of fluid method for the dynamics of free boundaries”, Journal of Computational Physics 39 (1981) 201–225.

[32] Son., G., “Efficient implementation of a coupled level-set and volume-of-fluid

(29)

method for three-dimensional incompressible two-phase flows”, Numerical Heat Transfer, Part B; Fundamentals 43 (6), pp. 549-565, 2003

[33] M. Sussman, P. Smereka, S. Osher, “A level set approach for computing solutions of incompressible two-phase flows” , Journal of Computational Physics 114 (1994) 146–159

[34] Zheng, L.-L., Zhang, H., An adaptive level set method for moving-boundary problems： Application to droplet spreading and solidification, Numerical Heat Transfer, Part B; Fundamentals 37 (1-4), pp. 437-454, 2000

[35] Son., G., Hur, N., “A coupled level set and volume-of-fluid method for the buoyancy-driven motion of fluid particles”, Numerical Heat Transfer, Part B;

Fundamentals 42 (6), pp. 523-542, 2002

[36] Son., G., “A level set method for incompressible two-fluid flows with immersed solid boundaries”, Numerical Heat Transfer, Part B; Fundamentals 47 (5), pp.

473-489, 2005

[37] R. J. Leveque, K. M. Shyue, “Two-dimensional front tracking based on high resolution wave propagation methods “, Journal of Computational Physics 123 (1996) 354–368.

[38] R. Saurel, and R. Abgrall, A simple method for compressible multifluid flows, SIAM J. Sci. Comput. 21 (3) (1999) 1115–1145.

[39] Saurel, R., and Cocchi, J.P., “Numerical Study of Cavitation in the Wake of a Hypervelocity Underwater Projectile,” Journal of Propulsion and Power, Vol. 15, No. 4, 1999, pp. 513-522

[40] D. Pan, C. H. Chang, “The capturing of free surfaces in incompressible multi-fluid flows”, International Journal for Numerical Methods in Fluids 33 (2000) 203–222.

[41] Yang-Yao Niu, “Simple Conservative Flux Splitting for Multi-component Flow Calculations”, Numerical Heat Transfer, Part B, Vol.38 No. 2, pp.203-222, September, 2000

行政院國家科學委員會專題研究計畫 成果報告