通用傳輸方程式

質量、動量與能量的守恆式可以用一具有物理性質 的傳輸方程式表 示。此控制體積的通用傳輸方程式示意圖如圖 2.1，其數學表示式：

c D s

s s s

d F d S F d S Q d Q d S

t  _

 

         



    

(2.1) 其中 t 表示時間，是控制體積，S 是控制表面，F_c  V 表示進出邊界 的對流通量，F_D則為進出邊界的擴散通量，^V 為流體速度，而Q_表示控 制體積內的源項，Q_s是邊界上的源項。然後藉由高斯散度定理，可將(2.1) 改寫為：

10

其中下標 1、2 分別代表兩個不同流體，



是流體 1 所占的體積分率。在

12

 

14

一壓力的不連續性，也就是P_v  P_l，故兩相之飽和溫度也不相等

16

2.8 邊界條件

(1) 進出口邊界條件

本研究中進出口採用壓力邊界，利用給定進出口不同壓力值形成壓力

差，此壓力差將會決定流量大小。而進出口速度計算採對流邊界條件，

c

0 t V

 

   



^，其中



為傳輸性質而V_c 為對流速度。

(2) 牆邊界條件

本研究中速度場的牆邊界條件採用無滑移(no-slip)邊界條件，即是(u=0, v=o)。而溫度場部分：定溫條件



^T=constant



^{或絕熱條件 T 0}

n

  

 

 

18

其中



_v 表示網格節點的 VOF 值，下標 i 表示網格節點周圍的網格，而

w

_i

20

格點的 VOF 即是根據式(3.7) 進行修正，修正後使主格點的 VOF 為１。

22

稱之為填充不足。此時主格點四個節點之 VOF 都大於 0.5，代表介面已 經通過主格點，但實際主格點的 VOF 卻小於 1。修正方法也是藉由權重 函數將下游網格的 VOF 往主格點補直到主格點之 VOF 為 1，式(3.6)為所 用的權重函數，並依照式(3.7)修正下游網格之 VOF。

3.2.3.4 枯竭不足(Under depleting, P⁰)

同樣在模擬 shear flow 時，如圖 3.4 介面尾端已經遠離主格點，理論 上主格點之 VOF 應該為 0，但實際主格點的 VOF 卻大於 0，稱此種情形 為枯竭不足。此時主格點之四個節點的 VOF 都小於 0.5，代表介面已經 離開主格點，但主格點的 VOF 大於 0。故將主格點多餘的 VOF 依權重 函數往下游網格分配，所用權重函數為式(3.6)而下游網格之 VOF 修正則 如式(3.8)。

3.2.4 VOF 計算流程

本研究使用 CISIT 的介面追蹤法求解 VOF 的計算流程為：介面位置 重建、求解 VOF 方程式 與進行 VOF 修正，其中每一時階將會進行兩次 的 VOF 修正步驟以確保計算範圍內每一個網格之 VOF 不會有殘留過大 或過小值。完成整個 VOF 計算流程後，除了介面所在的網格外，其餘網 格之 VOF 皆為 1 或 0 的均勻分佈。

3.3 動量方程式之離散化

24

 

_P

26

可得：

28

其中

30

2

32

法，最後整理為一線性代數式。

34

3.8 解題流程

在本章已經介紹 VOF 方程式、動量方程式與能量方程式之數值方法，

而雙流體流系統的解題流程如下：

1. 給定所有變數的初始值

2. 利用上個時階的體積通量求解 VOF 方程式得到新的體積分率



與重建 界面位置

3. 使用新的介面位置及上個時階所得溫度分佈求得界面處之熱傳量供相 變化計算使用，並且根據相變化造成的氣體體積變化修正此時階的 VOF 值

4. 求解動量方程式並且進行壓力修正以滿足連續方程式 5. 利用步驟 4 所得新的速度分佈求解能量方程式

6. 利用新的能量分佈求得溫度場分佈

7. 重複步驟 2 至步驟 6 至設定總計算時階數

36

第四章 結果與討論

4.1 相變化測試(一)，Stefan’s 1st problem

為了測試本研究所使用的相變化模型的準確性，將以兩個具有理論解

根據[21]，[22]此問題之 Neumann 解如式(4.5)及式(4.6)所示

 

t ² vt

38

漸與理論解保持相近距離。

4.2 相變化測試(二)，Stefan’s 2nd problem

接著進行另一個相變化測試，稱之為 Stefan’s 2nd 問題。此問題之示 維之 10mmx1mm 的均勻計算網格，共有 50x5、100x10 及 200x20 三組網 格大小，所用的計算時間步階為2 10 ^3 s。邊界條件設定：x=0mm 與

40

其中^

 

^{t 為介面位置，}為式(4.13)的根，_l為液相流體之熱傳係數，erfc 則 為 complementary error function。

首先以類似求解 Stefan 的方法測試採用何種混合熱擴散係數較為準 確，因此以三組網格進行三種熱擴散係數之測試，分別為調和平均、混合 型與液相流體之熱擴散係數。計算介面熱傳量所用

   n x

，計算所得結 果如圖 4.11 至圖 4.13 所示。可以發現到三組網格中，都以使用液相流體 之熱傳係數計算的結果最為接近理論解的介面位置。而從圖 4.14 至圖 4.16 之溫度分佈圖可以發現到使用液相流體之熱傳係數計算的溫度分佈 最為接近理論解的溫度分佈，這是因為氣相之流體均維持在飽和溫度，若 在單一流體模型中全部使用液相流體的熱傳性質可以改善在介面附近計 算格點因使用混合後熱傳性質造成的溫度不連續性。因此在本例子之特殊 情況下，完全使用液相流體的熱傳性質可以得到較接近理論解之溫度分 佈。

而圖 4.17 與圖 4.18 分別為三組網格中使用液相流體熱傳性質的介 面位置圖及介面速度對時間變化圖，可以發現以 200x20 這組網格的介面 位置最為接近理論解。從介面速度變化圖中可發現網格越密越可以表現出 在計算時間初期較高的介面速度。而隨著計算時間的增加，因溫度梯度的 減小使得介面速度也逐漸下降，其中以 100x10 與 200x20 這兩組網格之 結果較接近理論解。

在測試完熱擴散係數造成的影響後，吾人選定液相流體之熱擴散係數 搭配不同的

 n

進行測試，希望得知使用不同長度的

 n

計算介面處的熱傳 量對於相變化情形之影響。同樣使用了 50x5、100x10 及 200x20 三組網 格分別搭配一至三倍網格長度之

 n

進行計算。圖 4.19 為計算結果之介面 位置，可以發現在三組網格中隨著之

 n

增加介面位置移動變慢，其中三 組網格中、分別以 50x5 搭配 1.5 倍網格長之

 n

、100x10 搭配 1 倍網格 長的

 n

與 200x20 搭配 1 倍網格長的

 n

計算結果之介面位置最為接近理 論解。將其畫於圖 4.20。最後將計算結果整理於表三，其中以 200x20 使 用液體熱傳性質與搭配一倍網格長之

 n

計算結果的介面位置最為接近理 論解，誤差為 1.38%。

從圖 4.19 可以發現

 n

並非取越長越好，因為從圖 4.21 的溫度分佈圖 可以發現理論解之溫度分佈是往液相側擴散，若取過長的

 n

會使得求出 的溫度梯度較低進而減少相變化的質量，介面移動速度也就會較理論解慢。

因此需要選擇適當的

 n

避免影響計算相變化的準確性。之後研究將使用 1 至 2 倍網格長度作為

 n

的長度。

圖 4.22 為三組網格各自搭配最佳之

 n

介面速度對時間變化圖，其中 可以發現到 50x5 這組較粗的網格在計算初期無法準確表現出較大的溫度 梯度所以求得的介面速度較慢，也因此造成這組網格的介面位置始終落後

42

並且建議計算範圍接近其 most unstable Taylor 波長以便觀察氣泡生成，如 式(4.14)所示。

64x192(4 10 ^6 s)。邊界條件設定：左右兩側為對稱邊界條件，底部為定溫

44

變化圖，可以發現在計算時間 0.3 秒之前氣體體積的增加率較為趨緩，而 在 0.3 秒前後各組網格的氣體體積都開始快速增加，對應到氣泡生成過程 此時正是第一顆氣泡成形的時間。接著 64x192 這組網格約在 0.32 秒後氣 體總量增加速度開始減緩，而 16x48 與 32x96 這兩組網格與則約在 0.35 秒左右開始減緩。若對照圖 4.27 至圖 4.29 各網格之氣泡生成與速度場圖，

可以發現氣體總量體積開始減緩之時間點，都接近第一顆氣泡脫離底部的 時刻。

第一顆脫離的氣泡在脫離後受到表面張力影響開始被壓成月彎形狀，

而之後在其後面有第二顆小氣泡開始成型，慢慢的也會脫離底層。而在 32x96 這組網格可以發現到第二顆氣泡在脫離底部後，還會融入第一顆氣 泡中而形成一顆更大的氣泡。

若觀察圖 4.27 至 4.29 可發現類似的氣泡生成過程，首先由一開始給 定的微小擾動處開始有氣體增加，表示介面處存在一淨熱傳量足以克服相 變化所需的潛熱而開始產生相變化。隨著計算時間的增加在計算範圍的左 下角開始有越來越明顯的氣體生成現象並且造成速度場的變化，約在 0.20 至 0.25 秒間左下角開始有一明顯渦流產生，之後逐漸有一球型氣泡開始 在左側形成，而隨著氣泡逐漸成型時底部的氣體也形成薄層分佈情況。當 氣泡持續成長與上升時會形成一拉伸的氣柱。約在 0.35 至 0.40 秒間第一 顆氣泡脫離獨立成形。之後也會有第二顆小氣泡成形、拉伸最後脫離。

圖 4.30 是取計算範圍內左下角落網格中心的溫度與底部的壁溫 510 K 依式(4.16)求出此格點的

Nu

以表示底部之熱傳變化。



wall ⁰ sat



y ₀

Nu T

T T y





 

 

(4.16)若將 圖 4.30 與圖 4.27 至圖 4.29 對照來看，可以發現

Nu

的第一個極小值約 發生在 0.3 秒左右(16x48 為 0.32 秒)，此時正為氣泡從薄膜中生成逐漸成 形的時候，而隨著氣泡的上升

Nu

也隨著變大，約在 0.35 秒左右(16x48 為 0.38 秒)第一顆氣泡正要脫離時

Nu

也到達第一個極大值。而當氣泡脫離後

Nu

開始變小，以 32x96 這組網格為例，約在 0.45 秒時薄膜層又形成一氣 泡，此時

Nu

也來到第二個極小值的位置。之後隨著第二顆氣泡開始上升 同時

Nu

也開始上升，而

Nu

便在第二顆氣泡脫離前達到了第二個極大值 的位置。從以上現象可以發現當有氣泡在薄膜區形成時，

Nu

會處於相對 極小值；而當氣泡正要脫離薄層時，

Nu

會達到相對極大值。

46

第五章 結論

本研究利用一基於 VOF 的 CISIT 之介面追蹤法模擬一薄膜沸騰現 象。使用單一流體模型以一組質量、動量與能量守恆方程式的統御方程組 來描述此雙流體系統。由介面追蹤法所得之介面位置進一步求出介面處用 以驅動相變化產生之熱傳通量，並將相變化造成之質量變化量加入連續方 程式之源項，藉以滿足兩相間之跳躍條件。先以具理論解的相變化問題驗 證相變化之處理，可以發現溫度分佈之準確性與計算熱傳通量時所用之定 長

 n

影響相變化結果甚鉅。若選擇適當長度作為

 n

之設定，本研究之計 算結果(介面位置、介面速度)接近理論解。最後模擬了一薄膜沸騰現象，

除了可以觀察到氣泡的生成過程，此外同時計算底部某一定點的

Nu

，藉 此研究氣泡的生成情形與氣相內的熱傳性質變化的變化關係。可以發現每 當有氣泡在底部成形時，此時

Nu

會處於相對極小值；而當氣泡逐漸上升

除了可以觀察到氣泡的生成過程，此外同時計算底部某一定點的

Nu

，藉 此研究氣泡的生成情形與氣相內的熱傳性質變化的變化關係。可以發現每 當有氣泡在底部成形時，此時

Nu

會處於相對極小值；而當氣泡逐漸上升

在文檔中 利用一基於流體體積(VOF) 之介面追蹤法求解薄膜沸騰 (頁 20-0)