1、3 研究目的與方法

第一章 緒論

1、1 引言

空蝕現象(Cavitation)或稱為成穴現象，本質上為流場中一相當複 雜的機制，其所形成之特殊架構乃包含了氣相、液相與汽相等多相流 體(Multi-Phases)間交互作用的影響。一般而言，吾人所界定空蝕現象 的發生為，假設工作流體溶有氣體以及蒸汽之微小氣核(Nuclei)並懸浮 其中，當流場內某一處因為流道的突然擴張(Sudden Expansion)或者收 縮(Sudden Contraction)而呈現出局部的壓力下降至某一臨界值時（一般 為工作流體的蒸氣壓 P_v）。此刻，氣核將瞬間膨脹成為肉眼可見到之 氣泡(Cavities)，並且隨著流體的帶動而移往高壓區。氣泡則因為表面 受到壓力變化的影響會在流場中迅速的破裂，這種破裂通常伴隨著能 量的釋放並且引發流場的不穩定。

空蝕現象所導致的影響，其涵蓋的層面也幾乎囊括所有內、外流場 的領域；以外流場為例，高速快艇、潛艇、魚雷乃至於裝載有螺槳推 進器等流體機具在高速運轉的同時，都必須承受空蝕現象所導致的負 面影響，諸如，對機械表面造成侵蝕，整體性能的表現下降，誘導震 動以及噪音(Cavitation Noise) ..等等。而針對內流場而言，引擎中燃燒 室的注油孔(Injector Orifice)、氣渦輪機、工業用的燃燒爐、甚至噴墨 印表機的噴墨頭等工業產品，也因為空蝕現象的產生而有著燃燒不完 全，燃燒爐損毀、或著色位置錯誤等諸多不良後果。站在精密工業的 立場來看，實為一不容忽略的環節，應予以適當的防範。而本文則針 對內流場部分，希望藉由數值模擬的方法，並配合 ANSYS 計算軟體

來作一簡單的分析，藉由計算所得到之流場性質的分布情形，以期能 初步預測出空蝕現象所發生的位置與時機。

1、2 文獻回顧

空蝕現象的應用在處理相關工程流動課題上，一直是一重要又不 能被忽視的環節。根據 Knapp [1] 的文獻記載，空蝕現象的緣由可追 溯回 1894 年，驅逐艦在行進間螺槳運作至表面平均壓力為 11.75 psi 時，於螺槳表面產生一層近似真空的薄膜。這薄膜不僅影響到螺槳推 進的效率，也增加了不必要的油耗。此真空薄膜即便是空蝕現象的前 身。在外部流場的部分，機翼(Foil)和螺槳(Propeller)可以說是當今研究 的重點。以螺槳為例，螺槳葉片頂端因為渦度曳出所型成的漩渦空蝕，

在水力工程以及航太工程的領域範圍內，是個持續且不斷深入在探討 課題。其延伸的研究，有所謂的空蝕噪音現象。這現象就如同當直昇 機起飛或潛艇運作時，因螺旋槳與螺槳葉片頂端漩渦空蝕之交互作用 所產生的噪音，其噪音的特性如 Strsdberg [2] 所述。就以船舶的行進 為例，螺槳葉端的漩渦在整個流場中，是屬於相對較低壓的區域，一 但壓力遞減致工作流體的蒸氣壓時 (即空蝕現象發生)，大量的氣泡在 此低壓區域型成，而後隨著流場被帶開低壓區域。這些在低壓區域所 型成的氣泡一但脫離低壓區進入到一壓力相對較高的區域後，因為無 法持續的存在，氣泡將迅速破裂並且伴隨能量的釋放，除了造成流場 不穩定的擾動現象外也同時產生噪音。當然無可避免船體的震動與螺 葉之侵蝕，尤其是螺槳的尾端部分。因此，如何控制並延遲翼端漩渦 空蝕發生的時機，已成為首要探討項目，此基礎研究也成為船用流體 力學中一相當重要的課題。而針對機翼部分，因為空蝕現象通常發生

於翼後一或二倍弦長的範圍內，也由於流場可視量測技術的進步，LDV (Laser-Doppler)雷射都普勒儀試驗系統能允許我們更詳細的取得機翼 外部尾端漩渦的結構。近幾年來 Mazel. [3] 等人利用 LDV 量測超空蝕 翼周圍流場之速度分佈。接續 Mazel.，Brewer. [4] 等又利用動量原理 (Momentum Theory)進一步解出升力、阻力的參考係數。此外，根據 Scheiman. [5]，Stinebring. [6] 等研究報告指出，雷諾數、翼表面粗糙 程度、翼面攻角等因素也被認定為與空蝕的發生與否有直接或間接的 關聯存在。

而關於內流場部分，吾人關注的焦點則擺在燃油霧化或噴霧技術 的沿革為主，因為空蝕現象與液柱破裂和霧化現象之間具有一定程度 的相關聯性。換言之，空蝕現象的劇烈與否，會直接地影響到霧化的 品質[7] ，也由此間接決定出燃燒效率的好壞以及對環境污染的程度。

以引擎之注油器噴嘴為例，當燃油從注油噴嘴噴出的時候，噴流液柱 碎裂和噴霧之間的關聯相當複雜。D. T. Montgomery. [8]、W. Bergwerk [9] 算是最早初步完成在不同尺度的微小噴嘴內空蝕流體的研究。

Nurick [10] 的 研 究 發 現 指 出 ， 當 孔 口 (Orifice) 上 游 進 口 壓 力 達 到 25.8psia 時，孔口入口處將會產生一所謂的模糊區域(Fuzzy Region)，

且隨著壓力的增加，越往下游處移動，他認定其為空蝕現象開始發生 (Inception of Cavitation)。同時，他也指出此區域的發生與孔口處尖銳 的程度有一定比例上的關係。

1991 年 H. Hiroyasu. [11] 等人接續研究發現，當噴流的速度大於 某一臨界速度時，在孔口入口處會產生空蝕氣泡(Cavitation Bubble)，

此時靠近孔口內的流場因為氣泡而產生紊流的現象，在結合液柱表面 的擾動後，噴流液柱碎裂的距離將縮短；反之，若噴流的速度過大，

此時流場跳離的效應將呈現(Hydraulic Filp)，不但拉長噴流液柱碎裂的

距離，同時霧化的效果也就變的不明顯。

以上的研究大多以低速的流體來進行，並不符合現今的引擎設 備，著鑑於此，N. Tamaki. & M. Shimizu. [12] 等人開始針對高速流體 流經不同 L/D 比之單孔噴嘴，進行一系列的實驗研究，以了解其中空 蝕現象與霧化機制之關聯性，所得到的結論為，當噴流的壓差過大時，

空蝕不存在的狀況下霧化不完全，而且分離線加長；反之，當時空蝕 現象出現時，則有助於促進霧化的現象，分離線也會跟著變短，並且 提出因為空蝕而引發噴嘴孔口內紊流是影響霧化最主要的因素。

在數值部分方面，Delannoy. & Kueny. [13] 等人，利用 TVD 數值 方法和穩態 Barotropic 方程式，預測無黏性之空蝕流(Cavitating Flow) 流經文氏管，其結果符合文氏管的定量行為，但是卻不能確切地指出 空蝕泡移動的頻率或者更進一步量出空蝕泡的長度，因為其假設密度 為壓力的一連續 sin 函數，即均為純相不可壓縮。也因為這樣的假設，

在穩態解的前提下，必須把汽相與液相之密度比調低，如此反而影響 到解的正確性。

最近幾年，Schmidt. [14-17] 等人發展出一套 2-D、暫態的數值模 式，特別是研究有關微小尺度又高速的噴嘴。不過，其一併把可壓縮 流體 (液體)、(汽體) 一起納入考量，此外，為了避免密度在相變化時 產生的擾動 (在此指的是液體和氣體交界面而言)，程式中還加入一三 階震波捕捉來消除微小的擾動並得到不錯的結果。

國內從事相關方面的研究，有成大卓英吉.[18]，楊授印.[19]等人，

由其報告中可得知當雷諾數增加至某一程度時，回流區會出現在凸縮 或凸張的地區，再經由簡單的因次分析或者利用計算軟體的模擬，可 以發現通常這些區域都是流場之相對低壓區。此低壓值與雷諾數均呈 現出一反比的型態。此外，這些區域也幾乎是空蝕現象所發生的區域。

1、3 研究目的與方法

真實情況下空蝕現象的成因相當複雜，且大多數的研究都傾向於 對物體結構的破壞上面，針對管路等內流場的研究則較少著墨。在利 用數值模擬來作初步預測的前提之下，乃分別採用 ANSYS 計算軟體 與 SIMPLER 運算法則來作個別不同的處理。以期對內流場空蝕現象 的起訖有一大致的了解。

在選用的數值方法上，基於 S. V. Patanker [20] 的 SIMPLER 算則 有成功結合紊流乃至燃燒模式，並得到相當不錯的結果，因此本研究 也相同採用 SIMPLER 作處理。就凸張管件部分比對實驗值並進一步 抓出次回流區，來驗證程式的正確性。而注油器部分則針對收縮比 1.87 及 5 之圓柱型、具有尖銳邊緣的管件，於 ANSYS 分析中加入k − ε 紊 流模式來判斷空蝕現象發生與否和雷諾數、擴張比之間的關聯。混合 突張收縮管部分則導入質量傳輸(Species Transport)的效果，藉以找出 有關沖蝕預測之方法。最後，再配合相關的結果作一比對與整合。

第二章 物理問題

2、1 物理現象簡介

就空蝕氣泡(Cavitation Bubble)而言，其本質上為一暫態(Transient) 且極度不穩定的物理現象。氣泡從起始終致破裂的整個過程，會隨著 不同的幾何外型與流場速度，而有不同的變化產生。就穩態、不改變 幾何外型的前提假設下，考量管流內空蝕現象發生的主要原因。一般 而言，為流場中出現突然擴張(Sudden Expansion)或者收縮(Sudden Contraction)的情形，或者流場呈現出大幅度的轉彎所引發。因為流體 從截面積較大處流經截面積較小處時，流體會分離並產生一維納收縮 (Vena Contracta)的現象。在維持質量守恆的條件下，流場面積縮小，

其收縮處之速度勢必要加快。又以白努力定律(Bernoulli Equation)來 看，同條流線上能量守恆。一但收縮端之速度加快，其壓力勢必要下 降，致使流場中產生一局部的相對低壓區，以符合上述定律。此時若 再加大工作流體的流速，此低壓區也會隨之增大且延伸，一但壓力下 降至工作流體的蒸氣壓時，氣泡將會出現在流場中而型成所謂的空蝕 流(Cavitation Flow)。

依照 Knapp 的區分，氣泡的型態可以大略歸納為兩類，一種是所 謂氣體性的空蝕(Gaseous Cavitation)，其產生氣泡的特徵，一般而言，

體積比較小而且均勻密佈於流場中較低壓的區域，另外一種則為蒸氣 性的空蝕(Vapor Cavitation)，其特徵上比氣體性的空蝕所型成的氣泡來 的大，且彼此之間性質的差異較大。再者，此兩種空蝕現象的成因也 不盡相同。對氣體性的空蝕而言，定義上是指針對工作流體內一含有

固定濃度的成分氣體，其蒸氣壓略低於工作流體。而蒸氣性的空蝕，

則是泛指由於流場中局部之壓降，低於工作流體的蒸氣壓時其生成之 氣泡流。

2、2 物理問題假設

由 上 所 述 ， 吾 人 可 以 明 瞭 空 蝕 現 象 其 實 包 含 至 少 二 相 (Multi-Phases)、且相與相之間有複雜的關聯與變化。若欲完全以數值 模擬的方法，來求解逼近真實情形的近似解時，則必須至少求解相關 的二相流偏微分方程組。而本研究中欲以單相流來作一初步的模擬，

則必須倚賴以下的假設條件，來對所求解的問題進行減化。本文所使 用的假設條件如下所述：

1. 假設所考慮的流體以蒸氣性的空蝕為主，所以排除溶於工作流體中 的氣體，或者是從入口處(外界)夾帶氣體進入流場中。簡化問題單 就以蒸氣性的空蝕來作分析。

2. 流體性質屬於牛頓流體(Newtonian Fluid)的範圍，即剪應力與剪形變 速率成正比。

3. 所考量的流場為不可壓縮(Incompressible)、穩態(Steady State)，紊流 (Turbulance)模式。

4. 假設所有超距力如重力、磁力等可忽略不記。

5. 工作流體之比熱 (Cp) 及黏滯係數 (µ) 為一定值。

6. 不考慮浮力的效應(Buoyancy Effect) 。

第三章 數值模式

3、1 統御方程式的建立與轉換

二維穩態 Navier-Stokes 連續方程式及動量方程式依據 S. V.

Patanker [20] 所述撰寫成 SIMPLE Algorithm 之標準型態：

連續方程式

0 ) 1 (

)

( =

∂ + ∂

∂

∂ r v

r u r

x

ρ ρ

(1)

X-方向動量分量

))) )

1 ( 3 ( ( 2 )

1 ( ) (

) 1 (

) (

x rv u

r r x

x r v r r x u x

x p r

r u u u r r r x u u

x u

∂ + ∂

∂

− ∂

∂ + ∂

∂

∂

∂ + ∂

∂

∂

∂ + ∂

∂

−∂

∂ =

− ∂

∂ ⋅ + ∂

∂

− ∂

∂ ⋅

∂

µ µ

µ

µ ρ

µ ρ

(2)

R-方向動量分量

2 2

))) )

1 ( 3 ( ( 2 )

1 ( ) (

) 1 (

) (

r v x

rv u r r r

r r v r r r u x

r p r

r v v v r r r x v v

x u

µ µ

µ µ

µ ρ

µ ρ

∂ − +∂

∂

− ∂

∂ + ∂

∂

∂

∂ + ∂

∂

∂

∂ + ∂

∂

−∂

∂ =

− ∂

∂ ⋅ + ∂

∂

− ∂

∂ ⋅

∂

(3)

為了方便計算時不用考慮單位，將變數集中於雷諾數與尤拉數。統御 方程式之無因次化，其變數如下

∞

∞

∞

∞

∞

=

=

=

=

=

=

=

P P P

u R

r r

u

v v

u u R x x

0 0

ρ ρ ρ

µ µ µ

(4)

再重新整理統御方程式分別如下

連續方程式之無因次化

0 ) 1 (

)

( =

∂ + ∂

∂

∂ r v

r u r

x ρ ρ (5)

X-方向動量分量之無因次化

)]

1 ( ) 3 (

[ 2 Re

2

Re )]

( 2 Re )

( 2 1[

x r v r r x v x

u

x E p r

r u u v r r

x r v

u u x r

r ^u

∂

∂

∂ + ∂

∂ + ∂

∂

− ∂ +

∂

− ∂

∂ =

− ∂

∂ + ∂

∂

− ∂

∂

∂

µ

ρ ρ

(6)

R-方向動量分量之無因次化

] 2

) 1 (

1 ) 3 (

[ 2 Re

2

Re )]

( 2 Re )

( 2 1[

r 2

v r

u r

r v v r

r r x r u

r E p r

r v v v r r

x r u

v u x r

r ^u

− ∂

∂

∂

∂ + ∂

∂ + ∂

∂ +

− ∂ +

∂

− ∂

∂ =

− ∂

∂ + ∂

∂

− ∂

∂

∂

µ µ

µ

ρ ρ

(7)

其中無因次變數的意義：

1. Reynolds Number：(雷諾數)

∞

∞

= ∞

µ

0ρ Re 2u R

此為流體慣性力與黏性力之比，當雷諾數數值越大表示流體流速越 大，慣性力也相對增加。當雷諾數數值越小表示流體流速越小，黏性 力相對較大。其中 R₀ 為參考半徑。

2. Euler’s Number：(尤拉數)

∞2

∞

= ∞

u E _u P

ρ

此為壓力差和流體流動的壓力(static pressure)之比。

統御方程式之轉換

為了使計算方法能夠適用於所有複雜幾何外型之流場，物理平面 (Physical Domain X-Y)轉換至計算平面(Computational Domain ξ -η) 的過程將不容被忽視。亦即重合曲面座標系統與幾何外型，以便將複 雜之物理域映射到一矩形的計算域上。如此一來，在計算域中就能省 去對邊界幾何形狀上的考量，進一步避免瑣碎的內差而導致解的準確 性不足。除此之外，在邊界條件的處理上，也不需要因為更改幾何外

型而作大幅度的修改。更增添應用上的便利性。本文以垂直速度分量 (Contravariant Velocity Components)作為轉換變數，其轉換後通式如下：

S g

g

v g g

g g u

g g

∂ = Γ ∂

−

∂ + ∂

∂ Γ ∂

∂ −

∂

)

1 ( ) 1 (

22

11

η φ

φ η ρ

ξ φ φ

ξ ρ

φ

φ

(8)

帶入適當的變數得下列統御方程式：

轉換後之連續方程式

0 ) (

)

( =

∂ + ∂

∂

∂ g u g ρv

ρ η

ξ (9)

轉換後之 u 方向動量方程式

)

1 ( )

1 ( ) 1 (

22

11

u eff

eff

p S E g u

g

vu g g

g u g uu g g

∂ +

− ∂

∂ =

− ∂

∂ + ∂

∂

− ∂

∂

∂

ξ µ η

η ρ µ ξ

ξ ρ

(10)

轉換後之 v 方向動量方程式

)

1 ( )

1 ( ) 1 (

22

11

v eff

eff

p S F u

g y g

u vy g g

u g y

g u uy g g

∂ +

− ∂

∂ =

− ∂

∂ + ∂

∂

− ∂

∂

∂

η µ η

η ρ µ ξ

ξ ρ

ξ

ξ ξ

ξ

(11)

而其中E = X_ξ² +Y_ξ²、F = X_η² +Y_η²，S _u 、S _v 分別為 source term 在

u、v 方向上的分量。 g 為轉換後之 Jacobian 值、Γ 為轉換後之擴散 係數。

3、2 交錯格點配置

對於計算流體力學而言，在進入數值演算之前需先將計算空間離 算成數個有限的控制體積，接著利用每個控制體積內分別守恆 (質 量、動量、能量) 三大定律以及其交介面通量，配合欲求解的流場變 數形成非線性的插分式，然後才聯立求解。一般而言，數值運算會將 所欲求解的變數 (p、u、v) 全部計算於相同的網格點上，但這種處理 方式容易造成速度場震盪(Wavy Velocity)以及鋸齒狀的壓力場分布 (Checkerboard Pressure)。基於上訴理由，S. V. Patanker 建議採用交錯 網格(Staggered Grid)配置，把包括壓力、密度等純量變數儲存在主格 點，而速度向量則獨立置於與其垂直之控制體積的介面上，其使用上 的優點在於壓力因配置在速度之控制體積的控制面上，在數值觀點上 相鄰兩格點的壓力差正好是格點間速度分量的驅動力(Driven Force)。

如此一來不但避開鋸齒狀的壓力分布，且另外一方面速度項 (u ,v) 置 於邊界上，不需運用內插便可以直接求得，這使得速度項將更符和於 實際流場的真實狀況。

3、3 SIMPLER 數值運算模式

SIMPLE 演算法最大的特點在於解動量差分方程式的過程中，忽 略鄰近點速度修正項的影響。雖然這種做法並不會影響最終的收斂 值，但是維持等式兩邊平衡的條件下，忽略鄰近點速度修正項的值，

也幾乎是放大了壓力修正項。所以，為了改進壓力場表現不加的缺點，

SIMPLER (Semi-Implicit Pressure Linked Equation Reverse) 演算法於 是出現。只針對壓力修正方程式無法確切反映在壓力場上，作一變更，

將其利用來修正速度場，再由解出來的速度分布反推解得壓力場。也 就是 SIMPLER 中壓力場與速度場是相互協調的。所以，在計算的過 程當中，一次的遞迴將多了一壓力修正的方程式，理論上，計算時間 會有所增加，不過就整體而言遞迴次數卻比 SIMPLE 來的少，所以收 斂速度也較為快 [21]。

3、4 k − ε 紊流模式與處理

ANSYS 所外掛的紊流模式為將動量方程式中雷諾應力項(Reynolds Stress Term)以微分方程式或代數方程式來作表達。在剔除暫態項後，

ε

k −

model 動量方程中紊流動能 k 以及消散項ε ( Turbulence Kinetic Energy；Dissipation Rate) 分別如下所示：

ρε σ µ

µ

ρ ₋

∂

∂













∂

∂

∂ + ∂



 





∂

∂

∂

= ∂

∂

∂

j i i

j j i t i

k t i i

i

x u x u x u x

k x

x k u ) (

(12)

C k x u x u x u C k

x x

x u

j i i

j j i t

i k t i i

i

2 2 1

)

( ε µ ε ρε

σ µ ε

ρ −

∂

 ∂













∂

∂

∂ + ∂



 





∂

∂

∂

= ∂

∂

∂

(13)

其中 µ_t = C_µ ρ ^kε²

ρ

與 U 為 流 體 之 平 均 密 度 與 速 度 ， 而

µ

_eff 為 有 效 黏 滯 性 (Effective Viscosity)，其定義為

µ

eff

= µ + µ

t (14)

依 經 驗 所 需 要 的 常 數 值 分 別 為 C_µ =0.09_、 C₁ =1.44 、 C₂ =1.92 0

.

=1

σk 、σ_ε =1.3。

3、5

質量傳輸守恆方程式

假設管壁上有些微的離子進行解離。針對 SIMPLER 程式而言，加 諸於上的質傳方程式(Mass-Transport Based)如下：

1 0 )

( =





∂

− ∂

∂ + ∂

∂

− ∂

∂

∂

r D Y r vY r r

r x D Y x uY

i i i

i i

i ρ ρ ρ

ρ (15)

其中 Y _i 表某一成分之質量分率， D _i 為其質量的擴散係數。

3、6 邊界條件設定

對於 SIMPLER 程式而言，比照實驗給予一近似的邊界條件，其 給定條件如下：

進口條件：入口速度呈現拋物線型分布，v=0，Y_i = 0 出口條件：考慮質量守恆，其餘皆以 = 0

∂

∂ r

ϕ 處理，其中 ϕ=u v , ,ρ

對稱軸：v=0，其餘皆以 = 0

∂

∂ r

φ 處理，其中φ = u ,Y_i

邊牆條件：假設邊牆無滑動(No Slip Condition)，即 u , v 大小值為零 Species 邊界條件：Y_i = Y_i,0

而對於 ANSYS 而言，所採用的計算元素為 2D-fluid 141、3D-fluid 142，其基本的元素性質參照使用手冊 [22]，其邊界給定如下：

進口條件：入口速度呈現拋物線型分布，v=0 出口條件：考慮質量守恆，其餘皆以 ^{= 0}

∂

∂ r

ϕ 處理，其中 ϕ=u v , ,ρ

邊牆條件：假設邊牆無滑動，即 u , v 大小值為零

3、7 求解之程序與收斂標準

SIMPLER 演算法，其解法的程序如下所述：

1. 先將壓力梯度項(Pressure Gradient Term)從動量方程式中剔除，並 由此動量方程式計算出虛假速度 u、v。

2. 聯立連續方程式和兩組動量方程式，並藉由推導出來的壓力方程式 求出每一格點的壓力值。

3. 解包含壓力梯度項的動量方程式，已獲得交錯速度(Starred Velocity)

u∗、

v

^∗ 。

4. 利用解得的交錯速度 u^∗和

v

^∗ 推得壓力修正方程式 p′(Pressure Correction Equation)。

5. p′ 只用來修正速度向量，但不修正壓力項。

6. 判別是否滿足收斂標準(Convergence Criterion)，若不滿足，則將所 求的變數值當成新的猜測值，並且重複一新的循環計算過程。

以上程序均為利用改變截面方向前後掃過的 Line by Line TDMA 來解 矩陣化的代數方程式。最後將所解出的變數從計算域(Computational

Domain)轉換至實際域上(Physical Domain)。

收斂標準

理論上，就 SIMPLER 而言，經由疊代方法求解統御方程式時，等 式兩邊值不可能完全相等，此時必定存有一誤差在。此時吾人設定一 殘存值(Residue)來終止疊代。

∑ ⁻

−

= a a b

R

_{es p}

φ

_{p nb}

φ

_nb (16)

採 用 的 收 斂 標 準 應 視 不 同 物 理 性 質 而 有 所 不 同 ， 本 文 所 有 利 用 SIMPLER 之計算條件均以下降 12 個 order 為主。

而 ANSYS 則考量計算方法之不同 (ANSYS 為積分，SIMPLER 為 微分)，其收斂標準與殘存值的定義參照 [22]，本文所有利用 ANSYS 之計算條件均以下降 8 個 order 與 SIMPLER 作一近似比對。

第四章 結果與討論

4、1 測試突張管件之流場

本文先以突張管之穩態流場作一簡單測試。突張管件的研究中，

分別在層流、紊流等不同的條件下，與實驗所得到流場之物理現象來 個別驗證 SIMPLER 程式和 ANSYS 計算結果的正確性和差異程度。首 先，針對突張管格點部分，為了滿足格點獨立(Grid Independency)使所 求得的解不至於影響答案，以及在增進計算效率的考量之下，格點數 乃採用 150 ×23、300 ×46來比較格點獨立解影響雷諾數與迴流區長 度，其計算網格示意圖如圖 4.1、4.2 所示，SIMPLER 程式部份所用的 網格如圖 4.3 所示。比較格點獨立影響雷諾數與迴流區長度的關係圖 如圖 4.4 所示，從圖 4.4 上可以清楚的了解粗格點與細格點影響迴流區 大小的差異並不大。格點數採用150×23其解應屬合理範圍，故可當作 為代表解。圖 4.5 為突張比 1.5 時雷諾數與迴流區長度的關係圖，由圖 上可得知，當流場為層流時，其雷諾數R _e 與迴流區長度X _r 的關係會 隨著雷諾數的上升而增長。一般而言，雷諾數越大表示工作流體的流 速越快，也由於流速的加快，相較之下流體的黏滯力道也越趨於薄弱 (即流體慣性力 Internia Force 大於黏滯力 Viscosity Force)，經與實驗值 比對之後，可以發現到此一趨勢相當。接續上述的理論，吾人繼續增 加突張比致 1.94 時，其雷諾數與回流區長度的關係圖。如圖 4.5 所示，

演算後出來的結果則驗證了 Durst 的說法。Durst 認為在平面突張管 中，隨著雷諾數的漸增，於突張處附近壁面兩側會出現不對稱的迴流 區。其示意圖如下 4.6 所示，比較的結果如上圖 4.5 所示。相較於兩種

不同的突張比，吾人可以得到一規律，當突張比加大的時候，在同樣 管長的前提下，次迴流區的現象較明顯。其原因可能為，當管路截面 積變大時，流體分離的現象較為明顯，且於突張管後段將產生一低壓 的迴流區(Recirculation Zone)。從能量的觀點出發，迴流區的現象可解 釋為能量的滯留，由流線圖 4.7 上可以清楚的看出背階處的主迴流區 呈現一順時針旋轉，而於上壁面將會產生一逆時針旋轉的次迴流區來 維持能量守恆。

不過經與實驗值比對其無因次化的背階長度之後，有些微的差異 存在。就 SIMPLER 程式而言，較 ANSYS 的計算結果更貼近於實驗值。

雖然兩著的趨勢相同，不過也都會隨著雷諾數的上升而加大誤差。此 差異的推測可能由於，當流場處於高雷諾數的狀態時，其本質上可能 已經脫離層流的模式，亦即有可能為暫態或者是紊流模式。從能量的 觀點出發，高速流體在流動的時候，其動量會傳遞至邊界層內，也因 此所得迴流區長度較實驗數據為短，而 ANSYS 的誤差則有可能為積、

微分的處理模式不同而導致數值上的誤差。(ANSYS 採用有限元素 法，SIMPLER 採用有限體積法)。

經由上述兩種不同擴張比之突張管的流場分析中，吾人可以比較 出 ANSYS 有限元素法、SIMPLER 有限體積法，其計算出答案的精確 程度與結果。雖然在準度上 ANSYS 不若 SIMPLER 來的精準，不過在 趨勢預測上仍可提供一近似的解答以便參考。此外，ANSYS 商業軟體 也提供使用者自行加入算則處理相關問題，但此以超出本研究之範 圍，需由後人再接續處理。

圖 4.1 突張管 ANSYS 之計算網格示意圖 1 (150×23₎

圖 4.2 突張管 ANSYS 之計算網格示意圖 2 (300 ×46₎

圖 4.3 突張管 SIMPLER 之計算網格示意圖(150×23)

0 200 400 600

Re NO.

0 2 4 6 8 10

Xr/h

Reasonable grid distribution---150*23 Refine grid distribution---300*46

圖 4.4 格點獨立解影響雷諾數與主迴流區長度圖

圖 4.5 突張比 1.5 雷諾數與主迴流區長度之關係圖

圖 4.6 突張比 1.94 雷諾數與主(次)迴流區長度之關係圖

圖 4.7 突張管流場於(Re = 600)之示意圖

4、2 定義空蝕係數(Cavitation Number)

在進入

空蝕

分析之前，吾人需先定義一空蝕係數以簡化分析的程 序。

我們先行定義一點P_O 為上游(Upstream)的一點，其表示遠端一速 度接近於零的地方(Far Field)。P_c 位於維納收縮段(Contraction)的中間 點， P _e 則位於收縮段出口處(Exit of Orifice)，如圖 4.81 所示。我們 須假設當空蝕發生時會影響流體流動,且其對於孔口的截面積有固定 的比率，這比率定義為C_C (Coefficient of Contraction)，流體流過P_c 點 之截面積為 A_C ，而孔口之截面積為 A，所以

A

C _C = A ^C (17)

假設流體密度固定且沒有相的質傳(Inter-Phase Mass Transfer)，所 以流量定義為：

C

C V

A

m_& = ρ ⋅ ⋅ (18)

流體流經收縮段時，剪力可以忽略，假設由P_O 到P_c 之間為一個 勢流場(Potential Flow) (如圖4.8a)，因兩點之間的壓力差，可以有效的 穩定流場與降低亂流的產生(Suppress Turbulence and Stabilize Flow)。

又在 P_c 分別存在液體與蒸汽，所以我們假設P_c 壓力等於蒸汽壓，所 以由P_O 到 P_c 之間的 Bernoulli's equation :

2

2 1

C V

O P V

P = + ρ (19)

由(17)~(19)可推得：

) (

2 _{o V}

C P P

AC

m_& = ρ − (20)

又定義一分離係數 C_d (The Coefficient of Discharge)：

e o

d P P

V

C = −

2

2 1 ρ

(21)

所以由上式，吾人可得：

e o

V o

d P P

P C P

−

= − (22)

此時空蝕參數定義為 K =

e o

V o

P P

P P

−

− (23)

因此在高壓噴流機構中，K值會趨近於 1，當入口與出口壓力均 降低時，K值會變的非常大。當入口壓力減少時或出口壓力增加時，K 值也會隨之而變大。

4、3 注油器之流場與空蝕分析

經由突張管的基礎驗證並確定 SIMPLER 程式與 ANSYS 的正確性 與可行性之後，吾人接續探討粗細段收縮比為 1.87 以及 5 的兩種不同 注油器內其流場分布與壓力下降的情形 (注油器均假設為圓柱型、具 有尖銳邊緣)，並期簡單預測出空蝕現象所發生的區域。流場示意圖如

圖 4.8b 所示。大略區分為層流與紊流兩種模式分別進行探討。其計算 網格示意圖如圖 4.9 所示。圖 4.10 為 SIMPLER 之計算網格示意圖。

不論粗細段比為何，其流場的分布狀況則相當類似，均在同一側邊之 轉角處、與收縮段會出現迴流現象。而其大小則取決於幾何外型的收 縮比，以及其工作流速的大小。如圖 4.11 表示注油器其噴嘴收縮處雷 諾數與迴流區長度的關係圖，由圖上可以看出此迴流區長度隨著雷諾 數的增加而變長。且迴流區的起始位置也緩慢往後移動，其流場之流 線示意圖及放大圖如圖 4.12、4.13 所示。理論上，工作流體從寬邊尖 銳處入口流進窄邊前端時，由於幾何外型上的尖銳轉角，迫使連續的 流體分離成兩部份，而導致所謂維納收縮現象(Vena Contracta)產生於 噴嘴入口處，此即為上述迴流區的成因。而此迴流區只在雷諾數約大 於 300 時才會出現。再者以動量守恆的觀點而言，迴流區的出現也導 致流場截面積的縮減，一但截面積減少，流體流經維納收縮處後端時，

流速勢必增加，更加速壓力的下降程度。若此時壓力降至工作流體的 蒸氣壓時，空蝕現象就有可能發生。而此壓力相對極小值出現的位置，

也正好是流場中首先出現低壓的地區。不過在 1.87 的計算例中其收縮 比不大，雖然收縮段壓力值會隨著雷諾數增加而變小，如圖 4.14 所示，

但並不足以降至工作流體之蒸氣壓以下(3.5674 *10⁵ Pa)。在相同的操 作環境下，當收縮比加大致 5 時，雷諾數增加至約 8000 時，才會有空 蝕現象出現。如圖 4.15 所示。由以上研究結果發現到突縮比越小的管 路，雷諾數必須放到越大才會有空蝕現象發生的可能，由圖 4.16 空蝕 係數與雷諾數的關係可以了解，雷諾數在 8000 以下幾乎沒有空蝕現象 發生 (因為壓差未達標準)。所以，在提高雷諾數的前提下，吾人必須 利用紊流模式以得到更精確的答案。然而，在現有的 SIMPLER 程式 不足之情況下，僅能以層流模式來逼近紊流的解，並配合 ANSYS 來

作驗證。此外，針對 SIMPLER 程式中加入質量傳輸方程而言，探討 的焦點則擺在加入質傳效應後對迴流區大小尺度有無影響。如圖 4.17~19，由圖中可知，當雷諾數越高時，相對被工作流體帶離壁面的 離子也越多。不過，從現階段的研究中發現質傳效應的加諸與否對於 迴流區的大小影響並不大。圖 4.20~21 為突張比為 5 之注油器其壓力 分別於雷諾數 4000、8000 之分佈圖。從其分佈結果，吾人可以觀察出 一趨勢，即當雷諾數越大，其低壓區域也會隨之增大。

圖 4.8a 定義維納收縮示意圖 Reattachment

Vena Contracta Expansion and

Reattachment Zone

Contraction Zone Entrace

o c e

Exit

圖 4.8b 注油器之流場示意圖

(摘錄自 Fundamentals of Fluid Mechanics [25] )

圖 4.9 注油器之 ANSYS 整體格點示意圖(150×43₎

圖 4.10 注油器之 SIMPLER 格點示意圖(150 ×23)

圖 4.11 突張比 1.87 注油器噴嘴收縮處雷諾數與回流區長度之關係圖

圖 4.12 突張比 5 注油器於(Re = 800)流場之流線示意圖

圖 4.13 突張比 5 注油器於(Re = 800)流場之流線放大圖

圖 4.14 突張比 1.87 注油器於噴嘴附近之壓力與雷諾數關係圖

圖 4.15 突張比 5 注油器於噴嘴附近之壓力與雷諾數關係圖

圖 4.16 突張比 5 注油器雷諾數與空蝕係數之關係圖

圖 4.17 突張比 1.87 注油器於(Re = 400)之質傳分佈

圖 4.18 突張比 1.87 注油器於(Re = 600)之質傳分佈

圖 4.19 突張比 1.87 注油器於(Re = 800)之質傳分佈

圖 4.20 突張比 5 注油器於(Re = 4000)之壓力分佈

圖 4.21 突張比 5 注油器於(Re = 8000)之壓力分佈

4、4 混合突張收縮管之流場與空蝕分析

經由以上兩個例子，吾人認為空蝕現象發生的條件取決於壓力 差，而壓力差異的大或小則決定在幾何外型與其工作流體的速度上。

而發生空蝕的位置均位於幾何外型變化劇烈處附近 (出口或收口處)。

針對混合型管件而言，要發生空蝕現象，由於出口處多一寬管的效應，

其壓力降不如注油器於同位置來的大 (假設注油器與混合管件突張處 前端之長度相似)。所以，預明瞭混合管件發生空蝕的時機與機制，必 須在高雷諾數下以紊流模式來進行分析。其幾何外型示意圖與計算網 格示意圖如下圖 4.22、4.23 所示 (混合突張收縮管均假設為圓柱型、

具有尖銳邊緣)。大體上，流場的分布狀況相當類似前述突張管與注油 器兩者之合併。吾人以突張處的流場分布來作一探討，由圖 4.24、4.25、

4.26、4.27 為雷諾數分別在 2000、4000、8000、12000 下縮張比 5 之 混合管件尾端部分的流線示意圖可以看出，當 ANSYS 加掛紊流模式

ε

−

k model 後，其尾端會呈現出一對稱的分布，且此迴流區的長度隨 著雷諾數增加而加大。而由圖 4.28 流線示意圖上吾人可以了解。當縮 張比加大致 5 時，雖然管徑變窄，不過，因為後端多一突張管所致，

其壓差不如注油器管件來的大。因此，於收縮段壁面處因流體黏滯性 與摩擦力之影響，所造成負的壓力梯度不足以抵擋使噴流前進之壓 力，流場分離所導致之維納收縮現象反而不明顯，如圖 4.29 收縮段流 線放大圖所示。

不過在縮張比 2 的計算例中由於縮張比不大，雖然收縮段壓力值 會隨著雷諾數增加而變小，但並不足以降至工作流體之蒸氣壓以下 (3.5674 *10⁵ Pa)。在相同的操作環境下，當縮張比 2 加大致 5 時，雷

諾數增加至約 12000 時，才會有發生空蝕現象的機會。如圖 4.30、4.31 所示。由以上研究結果發現到突縮比越小的管路，雷諾數必須放到越 大才會有空蝕現象發生的可能，由圖 4.31 空蝕係數與雷諾數的關係可 以了解，雷諾數在 12000 以下幾乎沒有空蝕現象發生 (因為壓差未達 發生條件)。

在加入質量傳輸方程式後，如圖 4.32~34，由圖中可知，當雷諾數 越高時，質傳效應也越明顯，即相對被工作流體帶離壁面的離子越多。

跟注油器情形類似，質傳效應的加諸對於迴流區的大小影響並不明 顯。圖 4.35 為縮張比 5 之混合管於 Re = 12000 之收縮段壓力分佈圖，

其分佈情形與注油器相仿，但以達發生空蝕現象的標準。

至於縮張比 5 之混合管 3D 部份，其格點分佈情形取近似 2D，如 圖 4.36~38 所示，其流場於 Re = 12000 之分佈如圖 4.39~4.42。從其分 佈情形來看，大體上整體趨勢與 2D 類似。不過針對細部迴流區大小，

來作一比較 2D 與 3D。吾人發現在尺度上，3D 較 2D 來的小一些，如 圖 4.43 二維與三維後段迴流區大小所示。此即所謂的三維尺度效應。

由本研究可得知，當雷諾數偏高時，因為多了一維的考量，粒子隨著 流線的變化會有旋轉的現象產生。流場後端的分佈情形也較為紊亂。

圖 4.22 縮張比 5 混合管之 ANSYS 計算網格示意圖(215×30)

( 2 D ) 0 9 J u l 2 0 0 1 

0 2 0 4 0 6 0 8 0

X 0

1 2 3 4 5 6 7 8 9 1 0

Y

( 2 D ) 0 9 J u l 2 0 0 1 

圖 4.23 縮張比 2 混合管之 SIMPLER 計算網格示意圖(215×30)

圖 4.24 縮張比 5 混合管於(Re = 2000)突張處流線放大圖

圖 4.25 縮張比 5 混合管於(Re = 4000)突張處流線放大圖

圖 4.26 縮張比 5 混合管於(Re = 8000)突張處流線放大圖

圖 4.27 縮張比 5 混合管於(Re = 12000)突張處流線放大圖

圖 4.28 縮張比 5 混合管紊流之收縮處於(Re = 12000)流線放大圖

圖 4.29 縮張比 5 混合管紊流收縮處於(Re = 12000)流線再放大圖

圖 4.30 縮張比 5 混合管雷諾數與壓力之關係圖

圖 4.31 縮張比 5 混合管雷諾數與空蝕係數之關係圖

圖 4.32 縮張比 2 混合管於(Re = 1000)之質傳分佈

圖 4.33 縮張比 2 混合管於(Re = 2000)之質傳分佈

圖 4.34 縮張比 2 混合管於(Re = 12000)之質傳分佈

圖 4.35 縮張比 2 混合管於(Re = 12000)之壓力分佈

圖 4.36 縮張比 5 混合管 ANSYS 三維計算網格示意圖(215×30×20)

圖 4.37 縮張比 5 混合管之 ANSYS 三維計算網格前段部份

圖 4.38 縮張比 5 混合管之 ANSYS 三維計算網格後段部份

圖 4.39 縮張比 5 混合管之 ANSYS 三維整體流場於(Re = 12000)分佈

圖 4.40 縮張比 5 混合管之 ANSYS 三維前段流場於(Re = 12000)分佈

圖 4.41 縮張比 5 混合管之 ANSYS 三維後段流場於(Re = 12000)分佈

圖 4.42 縮張比 5 混合管之 ANSYS 三維前段流場於(Re = 12000)分佈 放大圖

圖 4.43 縮張比 5 混合管於(Re = 12000)突張處二維與三維流線放大比 較圖

第五章 結論

經由上述的研究結果，吾人可以大致歸納出以下幾點結論：

1. 管流結構中流場最先出現相對低壓的位置，大致都發生在流場開始 變換大方向的位置附近。依據白努力定律，伴隨流場速度加快的結 果，即是壓力降的增加。倘若壓力降至工作流體的蒸氣壓時，空蝕 現象也就因此引發，故可以進一步推測此處為空蝕現象最先發生的 位置。

2. 當流場的速度越快時，於收縮處或突張處所產生的回流區將相對應 變長。一但流場進入紊流的範圍時，主回流區的長度將維持一定值。

3. 在突張比或收縮比過小的情況下，空蝕現象不易發生。以 1.87 開口 之注油器為例，必須將流場制定在紊流模式(於高雷諾數下)才會抓 出空蝕現象的可能。

4. 質量傳輸守恆方程式的加諸與否，對本計算例而言，當雷諾數越高 時，質傳效應將越明顯。此趨勢不論是應用在注油器或者混合管件 均呈現出一致性，然而從現階段的研究中發現質傳效應對回流區長 短的影響並不大。

5. 當空蝕現象出現時，其空蝕係數大致為趨近 1 的數.

6. 當粗細段內徑固定時，其細段的長度越短，越容易發生空蝕現象.

參考文獻

1. R. T. Knapp, J. W. Daily, and F. G. Hammitt, “Cavitation,”

McGraw-Hill Book Company, 1970.

2. M. Strsdberg, “Propeller Cavitation Noise After 35 Years of Study, ” ASME Symposium “Noise and Fluids Engineering,” pp. 88-99, 1977.

3. S. A. Kinnas and C. H. Mazel, “Nurmerical Versus Experimentail Cavitation Tunnel,”Journal of Fluids Engineering, Vol. 115,pp. 760-765, 1993.

4. S. A. Kinnas and Wesley H. Brewer, “Expriment and Viscous Flow Analysis on a Partially Cavitating Hydrofoil,” Journal of Ship Research, Vol. 41, No.3, pp. 161-171, 1997.

5. J. Scheiman, J. L. Megrail, J. P. Shivers, “Exploratory Investigation of Factors Affecting the Wing Tip Vortex,” NASA TMX-2516, 1972.

6. D. R. Stinebring, K. L. Farrel, M. L. Billet, “The Structure of a Three Dimentional Tip Vortex at High Reynolds Number,” Journal of Fluids Engineering, ASME, Vol. 113, pp. 496-503, 1991.

7. 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.

8. D. T. Montgomery, M. Chen, C. T. Chang, P. V. Farrell, and R. D.

Reitz, “Effect of Injector Nozzle Hole Size and Characteristics and the Performance of Heavy Duty D.I. Diesel Engine,” SAE Paper 962002, 1996.

9. W. Bergwerk, “Flow Pattern in Diesel Nozzle Spary Holes,”

Proceedings of the Institute of Mechanical Engineers, Vol. 173, 1959.

10. W. H. Nurick, “Orifice Cavitation and Its Effect on Spary Mixing,”

Journal of Fluids Engineering, pp. 681-687, 1976.

11. 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.

12. N. Tamaki, K. Nishida, M. Shimizu, and H. Hiroyasu, “Effect of the Internal Flow in a Nozzle Hole on the Breakup Processes of Liquid Jet,” ILASS AMERICAS 96, 9^th Annual Conference on Liquid Atomization and Spary System, pp. 255-259, 1996.

13. 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.

14. 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.

15. 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.

16. David P. Schmidt, M. L. Corradini, “One-Dimensional Analysis of Cavitating Orifices,” Engine Research Center, 1500 Engineering Dr.

University of Wisconsin, Madison. WI 53706.

17. 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, 1996.

18. 卓英吉, “以數值方法探討注油器管內流場之空蝕現象,” 國立成功 大學航空太空工程研究所碩士論文, 1997.

19. 楊授印, “以數值方法探討噴流機構內流場空蝕現象,” 國立成功大 學航空太空工程研究所碩士論文, 1998.

20. S. V. Patankar, “Numerical Heat Transfer and Fluid Flow,” McGraw-Hill, New York, 1980.

21. G. D. Raithby and G. E. Schneider, “Numerical solution of problems in incompressible fluid flow; treatment of the velocity-pressure coupling,”

Numerical Heat Transfer, Vol. 2, pp. 417-440, 1979.

22. ANSYS, CFD FLOTRAN Analysis Guide, Release 5.5, 1998.

23. B. F. Armaly, F. Durst, J. C. F. Pereira and B. Schonung, “Experimental and theoretical investigation of backward-facing step flow,” J. Fluid Mech, Vol. 127, pp. 473-496, 1983.

24. M. K. Denham, “The Development of a Laser Anemometer for Recirculating Fluid Flow Measurments,” Ph.D. Thesis, 1974.

25. B. R. Munson, D. F. Young and T. H. Okiishi, “Fundamentals of Fluid Mechanics,” 2^th, John Wiley & Sons Inc, 1994.