國 立 交 通 大 學
機 械 工 程 學 系
碩 士 論 文
利用一基於流體體積(VOF) 之介面追蹤法
求解薄膜沸騰
Use of a VOF based interface tracking method for film
boiling flow calculations
研 究 生：賴 胤 男
指 導 教 授：崔 燕 勇 博士
利用一基於流體體積(VOF) 之介面追蹤法求解薄膜沸騰
Use of a VOF based interface tracking method for film
boiling flow calculations
研 究 生：賴胤男 Student：Yin-Nan Lai 指導教授：崔燕勇 Advisor：Yeng-Yung Tsui 國 立 交 通 大 學 機械工程學系 碩士論文 A Thesis
Submitted to Department of Mechanical Engineering Collage of Engineering
National Chiao Tung University in Partial Fulfillment of the Requirements
for the degree of Master of Science
In
Mechanical Engineering July 2011
Hsinchu, Taiwan, Republic of China
利用一基於流體體積(VOF) 之介面追蹤法求解薄膜沸騰
研究生：賴胤男 指導教授：崔燕勇 博士
國立交通大學機械工程學系(研究所)碩士班
摘要
本研究使用單一流體模型搭配 CISIT 介面追蹤法求解雙相流之熱傳 與質傳現象。其中所用的離散方法為有限體積法 (finite volume method)， 而 CISIT 介面追蹤法則是利用求解流體體積分率的傳輸方程式來追蹤介 面之位置。利用介面處的質量與能量跳越條件求出相變化所造成的質量變 化，並將其加入連續方程式中的源項作為修正。將動量方程式與連續方程 式進行耦合，所得速度場將滿足介面處的質量跳躍函數。為了驗證相變化 模型的準確性先求解具有理論解之相變化問題，發現溫度場之準確性與介 面溫度梯度之計算都會影響模擬相變化的準確性。最後模擬一薄膜沸騰的 氣泡生成，可觀察到氣泡的生成與成長過程中對於熱傳及質傳變化的影 響。
Use of a VOF based interface tracking method for film boiling flow
calculations
Student：Yin-Nan Lai Advisor：Dr. Yeng-Yung Tsui
Department of Mechanical Engineering National Chiao Tung University
ABSTRACT
A single-fluid model with a VOF based interface tracking method (CISIT) is used for simulation of two-fluid flows with heat and mass transfer in this study. The finite volume method is used as the discretiztion method. The motion of the interface is tracking by the solution of a transport equation for the volume fraction. By using the interface position obtained from the interface tracking method the heat flux which drives the phase change on interface can be calculated. The mass change across the phase boundary is computed by taking the mass and energy jump conditions at the interface. It is added as a source term in the continuity equation. This continuity equation is used in the pressure-velocity coupling algorithm. After the velocity field obtained from solving the momentum equation the jump condition for the continuity is satisfied. To assess the accuracy of this mode two test problems for which analytical solutions are available, are considered. It is shown that the precision of the temperature profile and the length used for calculating temperature gradients at the interface would affect the accuracy of simulation. Finally, a simulation of film boiling on a horizontal plate is presented. It shows that how the bubble growth affects the heat and mass transfer of the flow.
誌謝
在這兩年的碩士生活裡，首先要感謝崔燕勇教授的悉心指導，讓 我學習如何去解決所遇到的問題與挑戰。並且崔老師在我遇到問題時 適時給予指導，讓我有個解決問題的方向指引。再來最要感謝林仕文 學長，學長不論在課業研究或日常生活方面都給予我許多教導與關心， 多虧學長的熱心協助我才能完成現階段的研究進度。也要感謝胡育昌 學長提供在 Fortran 語言的編寫經驗與技巧，幫我找出自己程式編寫 的錯誤。 此外也感謝陳信宏學長、郭大慶學長讓我很快的融入 307 實驗室 這個大家庭中。感謝林子翔與黃義政兩位同學在這兩年的陪伴與鼓勵， 你們認真的研究態度是激勵我不至於放棄的最大動力，謝謝你們！也 要感謝這兩年認識的同學與學弟妹們，與你們一起修課與玩樂是我難 忘的回憶。 接著要感謝立丘小集的張義雄老師，感謝老師讓我在學校領域外 學到許多寶貴且重要的處事觀念與態度。也要感謝李俊賦、陳二常、 張維仁和謝易珊等大學同學們在這幾年間的相互鼓勵與打氣。最後感 謝我的爸爸、媽媽與哥哥，謝謝你們一路的支持與鼓勵，讓我能得到 這兩年寶貴的成長與學習經驗。目錄
摘要 ... i ABSTRACT ... ii 誌謝 ... iii 目錄 ... iv 表目錄 ... vi 圖目錄 ... vii 符號表 ... ix 第一章 緒論 ... 1 1.1 前言 ... 1 1.2 文獻回顧 ... 2 1.2.1 求解雙流體流之介面問題 ... 2 1.2.1.1 Lagrangian 型式 ... 2 1.2.1.2 Eulerian 型式 ... 2(a) Front tracking method ... 3
(b) Level set method ... 3
(c) Marker and cell method ... 3
(d) Volume-of-fluid method ... 3 1.2.2 伴隨熱傳及質傳現象之介面問題 ... 4 1.2.2.1 僅考慮熱傳現象，無質傳(相變化)現象 ... 4 1.2.2.2 考慮熱傳現象與質傳(相變化)現象 ... 5 1.3 研究內容 ... 8 第二章 數學模式 ... 9 2.1 前言 ... 9 2.2 通用傳輸方程式 ... 9 2.3 質量及動量守恆方程式 ... 10 2.4 能量守恆方程式 ... 11 2.5 VOF 方程式 ... 13 2.6 相變化模型 ... 14 2.7 表面張力 ... 16 2.8 邊界條件 ... 17 第三章 數值方法 ... 18 3.1 前言 ... 18 3.2 VOF 方程式 ... 18
v 3.2.2 VOF 方程式之離散 ... 19 3.2.3 VOF 修正 ... 20 3.2.3.1 過度填充(Over filling, P1) ... 20 3.2.3.2 過度枯竭(Over depleting, P0) ... 21 3.2.3.3 填充不足( Under filling, P 1) ... 21 3.2.3.4 枯竭不足(Under depleting, P 0) ... 22 3.2.4 VOF 計算流程 ... 22 3.3 動量方程式之離散化 ... 23 3.3.1 非穩態項(unsteady term) ... 23 3.3.2 對流項(convection term) ... 23 3.3.3 擴散項(diffusion term) ... 24 3.3.4 源項(source term) ... 24 3.3.5 線性代數式之整理 ... 25 3.4 連續及動量方程式耦合處理 ... 26 3.5 求解動量方程式所用之流體性質 ... 31 3.6 能量方程式之離散化 ... 32 3.6.1 前言 ... 32 3.6.2 單一均質模型 ... 32 3.6.3 線性代數式之整理 ... 33 3.6.4 介面處之網格溫度處理 ... 33 3.7 介面處熱傳量之計算 ... 33 3.8 解題流程 ... 35 第四章 結果與討論 ... 36 4.1 相變化測試(一)，Stefan’s 1st problem ... 36 4.2 相變化測試(二)，Stefan’s 2nd problem ... 39 4.3 薄膜沸騰(film boiling)... 42 第五章 結論 ... 46 參考文獻 ... 47
表目錄 表一、一大氣壓之飽和純水之流體性質 ... 50 表二、Stenfan 問題之計算結果比較表 ... 50 表三、Sucking 問題之計算結果比較表 ... 50 表四、薄膜沸騰測試所用之流體性質( 5 1.0135 10 sat P Pa) ... 51
vii 圖目錄 圖 2.1 通用傳輸方程式示意圖 ... 52 圖 2.2 流體分布及曲率符號之示意圖 ... 52 圖 3.1 控制體積及其鄰格點之關係示意圖 ... 53 圖 3.2 對計算主格點之體積通量守恆示意圖 ... 53 圖 3.3 一介面前端跨越過一計算網格之示意圖 ... 54 圖 3.4 一介面尾端跨越過一計算網格之示意圖 ... 54 圖 3.5 計算主格點與鄰格點和計算面之示意圖 ... 55 圖 3.6 介面處熱傳量計算之示意圖 ... 56 圖 3.7 介面處熱傳量計算之示意圖 ... 56 圖 4.1 Stefan 1st 問題之示意圖 ... 57 圖 4.2 網格(50x5) 計算結果之介面位置圖 ... 58 圖 4.3 網格(100x10) 計算結果之介面位置圖 ... 58 圖 4.4 網格(200x20) 計算結果之介面位置圖 ... 58 圖 4.5 網格(50x5) 計算時間 t=10s 之溫度分佈圖 ... 59 圖 4.6 網格(100x10) 計算時間 t=10s 之溫度分佈圖 ... 59 圖 4.7 網格(200x20) 計算時間 t=10s 之溫度分佈圖 ... 59 圖 4.8 Stefan 問題計算結果之介面位置圖 ... 60 圖 4.9 Stefan 問題計算結果之介面速度對時間變化圖 ... 60 圖4.10 Stefan’s 2nd 問題之示意圖 ... 61 圖 4.11 網格(50x5) 使用不同混合熱傳係數計算結果之介面位置圖 ... 62 圖 4.12 網格(100x10) 使用不同混合熱傳係數計算結果之介面位置圖 ... 62 圖 4.13 網格(200x20) 使用不同混合熱傳係數計算結果之介面位置圖 ... 62 圖 4.14 網格(50x5)不同混合熱傳係數計算時間 t=2s 之溫度分佈圖 ... 63 圖 4.15 網格(100x10)不同混合熱傳係數計算時間 t=2s 之溫度分佈圖 ... 63 圖 4.16 網格(200x20)不同混合熱傳係數計算時間 t=2s 之溫度分佈圖 ... 63 圖 4.17 Sucking 問題計算結果之介面位置圖 ... 64 圖 4.18 Sucking 問題計算結果之介面速度對時間變化圖 ... 64 圖 4.19 液相流體之熱傳係數搭配不同n計算結果之介面位置圖 ... 65 圖 4.20 液相流體之熱傳係數搭配最佳n計算結果之介面位置圖 ... 65 圖 4.21 液相流體之熱傳係數搭配最佳n計算結果之溫度分佈圖 ... 66 圖 4.22 液相流體之熱傳係數搭配最佳n之介面速度對時間變化圖 ... 66 圖 4.23 薄膜沸騰模擬所給的初始介面分佈圖 ... 67 圖 4.24 薄膜沸騰模擬在不同網格之第一顆氣泡升起情況 ... 67 圖 4.25 薄膜沸騰模擬在計算時間 t=0.35s 不同網格之溫度分佈 ... 68 圖 4.26 氣相體積對初始氣相體積之比值對時間之變化圖 ... 68
圖 4.27 網格 16x48 之氣泡生成與速度場分佈圖 ... 69
圖 4.28 網格 32x96 之氣泡生成與速度場分佈圖 ... 70
圖 4.29 網格 64x192 之氣泡生成與速度場分佈圖 ... 71
ix
符號表
符號 定義 v c 等容比熱 fi C 通過 fi面之體積通量 e 內能 C F 面上的對流通量 D F 面上的擴散通量 k 熱傳導係數 m 質量流率 f m 介面上的質量流率 p 流體壓力 f S 面的法向量 w f S 液體浸潤截面之法向量 t 時間 v 流體速度 希臘符號 流體體積分率 Pnb P 到 nb 距離向量 曲率半徑 黏滯係數 密度 剪應力 表面張力第一章 緒論
1.1 前言 有關雙流體流的自由介面問題已經被廣泛研究，因為它不論是在於科 學或工業應用方面都極為重要。例如對於生物科技、船舶工程或石化工業、 鑄模工業等領域，如何準確預估自由介面的位置都是重要的研究議題。在 生物科技方面希望能藉由探討自由介面的運動，研究血管內物質傳輸現象。 而在氣液雙相流研究中，模擬沸騰氣泡的生成也是重要議題，首先必須準 確預測氣液介面位置。在石化工業方面希望模擬自由表面的運動，來幫助 設計出安全可靠的燃料載器。而鑄模工業希望以數值模擬作為開發新模具 的先行設計，若能準確有效的預測液態金屬在模具中的注入情況，以降低 過去試誤法所需較高的開模時間與模具成本。這些研究都希望獲得準確的 介面分佈，並且同時保持介面的鮮明性。 此外雙流體流問題常伴隨著熱傳現象，例如在模擬鑄模過程時，若同 時考慮熱傳現象，便可得知模具內的溫度變化，避免過高的局部溫度造成 模具的破壞，也同時能掌握金屬固化情形。而在模擬沸騰(boiling)現象時， 因氣液間溫度差造成的熱傳現象更是研究重點，因為若要準確的預測相變 化的情況，必須先準確地掌握介面處的熱傳量及質傳量。接著介紹有關求2 1.2 文獻回顧 1.2.1 求解雙流體流之介面問題 一般求解雙流體流的介面問題的數值方法可以分為 Lagranian 型式 (重建網格) 及 Eulerian 型式(固定網格) 兩種類型。 1.2.1.1 Lagrangian 型式 此種方法是在每個時階都會根據新的介面位置重建新的網格，確保介 面隨時都是準確預測並且保持其鮮明外型。而下個時階的介面位置是透過 計算已知的速度場而得到。當得到下個時階的自由表面邊界(free surface boundary)時，便開始重建網格並且更新流體的性質分佈。因為其邊界條件 是根據新時階之介面位置來設定，所以它能非常準確的預測介面移動。但 是當流場變化太過複雜時，如劇烈的變形或翻騰後造成的破碎表面，都會 降低其預測介面的準確性[1][2]。 1.2.1.2 Eulerian 型式 Eulerian 型式採用的是固定式的網格，所以不會有變形網格造成誤差， 但是其最主要的缺點就是會有數值擴散現象。這個現象會造成界面擴散至 數個網格上，如此介面不再明確。因此有許多學者提出不同方法希望降低 此類數值擴散現象。
(a) Front tracking method
此種方法是在計算開始前，將很多不具質量的追蹤質點均勻分布在初
始的介面上。而流體計算是在一 Eulerian 系統的進行，而追蹤質點相當位 於一 Lagrangian 系統，要即時將 Eulerian 系統的資料傳送到 Lagrangian 系 統是很複雜，而且當推廣到三維系統，追蹤質點的數量將會大增並造成計 算量與所需記憶體急遽增加[3]。
(b) Level set method
Level set 方法是以依連續函數(level set function)定義計算格點至界面 的最短距離，如此介面上的 level set function 值永遠都是零。然後經求解 level set 之傳輸方程式可以得到介面所在位置。此方法所求得介面為連續 分佈，但是此方法並未能保持質量守恆，當介面有劇烈變形時，會因為質 量不守恆而產生誤差[4][5]。
(c) Marker and cell method
MAC 方法是將許多不具質量的標記點均勻散佈在某一流體的區域， 而這些質點隨著計算過程變化而標定出此流體的介面分布。此方法優點是 能提供準確的介面定義，但隨著計算模型的複雜及計算時間的增加，而這 些紀錄質點將會消耗大量計算時間及計算儲存量[6]。 (d) Volume-of-fluid method VOF 方法是將某一流體的體積分率訂為 1 作為一指標函數(indicator function) ，而另一流體的體積分率則為 0，體積分率則介於 1~0 的網格即
4 到介面位置，並且 VOF 變化可滿足質量守恆。通常 VOF 方法可大致分作 介面捕捉法及介面重建法兩大類。 所謂介面捕捉法是利用不同的離散法去直接描述介面處的不連續性， 例如以一階上風法最為穩定但會造成數值擴散的現象使得介面不再明確； 而下風法雖不穩定，但具有反擴散(anti-diffusion)的特性可保持介面的形狀。 因此有許多學者提出了混合型式的介面捕捉法，結合高解析(易擴散)與壓 縮性(較不穩定)的離散法來準確預測介面運動並將介面擴散情形保持在數 個計算網格內[7]、[8]、[9]。 而介面重建法則是利用 VOF 的資訊在網格內重建介面的位置，如最 早 Noh 及 Woodward 學者所提出的 SLIC[10]，他們預測的介面僅有水平 或垂直兩種方向；之後 Young 學者提出 PLIC 能提供更為準確的介面位 置，但也更為複雜[11]。本研究中計算 VOF 方程式便是屬於介面重建法， 由崔燕勇及林仕文所提出的 CISIT (conservative interpolation scheme for interface tracking )[12]，將會在 3.2 節中介紹此方法。
1.2.2 伴隨熱傳及質傳現象之介面問題
1.2.2.1 僅考慮熱傳現象，無質傳(相變化)現象
1995 年，Pericleous et al.[13] 發展在一個部分為液體所填滿的空穴中，
及液體在空穴中所佔有的體積可以視為一在介面擁有極強烈不連續性的 性質的單一流體，假設無相變化。經由求解一個具兩個值的 scalar advection equation 而得到兩個流體各自所佔有的體積。此研究應用 Van Leer TVD 差 分方法在 scalar advection equation 上，同時在動量及焓方程式中應用通量 修正。最後求解兩個問題，首先為一在封閉並具牆上熱傳的空穴中的液體 柱的倒榻過程。第二個則是模擬液態鋁金屬填滿模具及瞬時的冷卻現象。
2002 年，Davidson and Rudman[14]發展一新的 VOF 方法計算一變動
的兩流體介面處的熱傳或質傳現象。假設無相變化，並且假設此兩相介面 無熱阻存在，所以介面上任一點皆可達熱平衡。在處理能量方程式時，是 個別以單位體的熱(q)去描述兩相的能量方程式，並且以一熱源相表示經過 兩相介面的熱傳量。求得兩相各自的 q 後再加總起來，進而求得溫度分佈。 在數值方法方面，處理對流項時使用體積通量來校正，而在兩相介面上使 用溫度及擴散熱通量的連續條件。此研究先測試一固定圓球體的熱擴散現 象(無對流效應)，其計算結果符合解析解。再進行(a)一環狀氣泡的成長與 其熱傳現象；(b)一成長的液滴(droplet)的質傳現象。 1.2.2.2 考慮熱傳現象與質傳(相變化)現象
2004 年，Esmaeeli and Tryggvason [15]以前端追蹤法(front tracking) 及
6
條件(jump conditions) 而將兩流體的質量、動量與能量方程式組合併成一 組統御方程組。其中在介面處因為跳躍條件造成兩相間動量及能量方程式 之差異，各增加一 delta function 形式之源項(source term)作為修正。此研 究進行了一些具理論解或與可與其他數值結果比較之測試，最後成功模擬 一水平平板上之薄膜沸騰(film boiling) 現象。
2007 年，Yang et al.[16]對於一 R141B 冷媒在盤管(coiled tube) 中的沸
騰現象分別以數值模擬及實驗量測方式進行研究。在數值模擬方面，使用 商業軟體 Fluent 6.0 中的 VOF 雙相流模型搭配蒸發模型的 UDF(user defined function)。並且使用 模型作為紊流模型，而在求解動量方程 式時使用 PISO(pressure-implicit with splitting of operators) 將速度及壓力 場進行耦合，在離散動量及能量方程式時使用二階的上風法(upwind)。將 數值模擬及實驗量測結果進行比較，由數值模擬預測之氣泡分布與實驗觀 測相當接近，也發現到溫度分布受到相分布之影響，例如高溫區皆發生於 氣相區。
2009 年，Sander and Weigand [17]使用 DNS(direct numerical simulation)
搭配 VOF 方法，模擬內含有液態冷卻金屬的引擎氣閥柱(hollow valve stem) 在引擎運轉的條件時，因受到加速度而產生的雙相熱傳現象。研究幾何形 狀、加速度或者液態冷卻金屬的填充程度等影響因子對於冷卻效果的影響。 模擬結果發現以閥內完全填滿液態冷卻金屬的冷卻效果最佳，但是過高或
過低的填充率在某些特定轉速時容易產生溫度急遽上升的 shaker-effect， 易造成機件損壞，因此較佳填充率為 30~50％。由於此研究模型十分複雜， 採用較多簡化模型，此研究結果僅作為定性分析較為適當。
2010 年，Haelssig et al. [18]採用單一流體(single-fluid)模型以一組包
含質量、動量、能量及物種方程式組來模擬一多流體組成之氣液流在介面 處的質傳與熱傳現象。因為作者將物種與能量在氣液介面處的跳躍條件 (jump condition) 都包含在統御方程式之源項(source term)，如此一來此方 法便可適用於多相物質的相變化情況。此研究使用商業軟體 Fluent 6.3.26 搭配 PLIC(Piecewise Linear Interface Calculation) 求解 VOF 方程式。而動 量方程式則以 PISO 求解。動量、能量及物種方程式採用二階的上風法進 行離散。作者利用此方法首先模擬乙醇(液態)-水(氣態)兩項流體在長 30mm 高 6mm 之二維管道中的質傳及熱傳現象。模擬結果的兩相間質傳 現象與實驗結果接近。
8 1.3 研究內容
本研究希望以有限體積法(Finite Volume Method) 搭配一新發展的 VOF (Volume of fluid) 型式的介面重建法求解雙流體流之熱傳與質傳現 象。 在模擬雙流體流時，採用單一流體(Single fluid) 的模型進行模擬，也 就是以一單流體代表原本的兩相異流體，此單一流體擁有連續的速度場及 溫度場分佈，而此單一流體之流體性質則是以流體體積分率的方式進行混 合。最後可用一組質量、動量及能量之統御方程式組描述此單一流體的模 型。 兩流體間的質傳現象(相變化)，先利用介面重建法中所得到的介面位 置計算出介面處的熱傳量，再利用介面處的質量與能量的跳躍條件(jump condition) 求出兩相間的質傳量。並且將相變化造成的質量變化量以一源 項型式加入連續方程式中，在求解動量方程式時將同時滿足連續方程式， 藉此模擬兩流體間的相變化情形。為了驗證相變化模型之準確性，本研究 先測試兩個具有理論解之相變化問題，最後模擬一薄膜沸騰之氣泡生成過 程。
第二章 數學模式
2.1 前言 本章將會描述計算模擬流場的數學模式。本研究中將兩個不相融的相 異流體視為一均質、連續的流體進行計算，但此兩種流體仍具有各自的流 體性質，在兩流體相接介面處保有不連續的流體性質。使用一步階函數型 式的體積分率表示兩流體的交接介面。此體積分率會影響流體性質分佈而 且定義出兩流體的介面分佈。在描述此雙流體(two-fluid)系統的數學模式 時使用以下基本假設：暫態、不可壓縮的黏性流、二維並且考慮重力影響 及表面張力。 2.2 通用傳輸方程式 質量、動量與能量的守恆式可以用一具有物理性質 的傳輸方程式表 示。此控制體積的通用傳輸方程式示意圖如圖 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 是控制表面，Fc V 表示進出邊界 的對流通量，FD則為進出邊界的擴散通量，V 為流體速度，而Q表示控 制體積內的源項，Qs是邊界上的源項。然後藉由高斯散度定理，可將(2.1) 改寫為：10 C D S d F d F d Q d Q d t
(2.2) 若考慮控制體積內一點，其微分型式的傳輸方程式可表示為： C D S F F Q Q t (2.3) 2.3 質量及動量守恆方程式 質量及動量守恆方程式，也可使用通用傳輸方程式表示。表示質量守 恆方程式時，物理性質
，並且省略擴散項及源項。而動量守恆方程 式中，物理性質
V ，而系統為非穩態且不可壓縮流，在處理擴散項使 用黏性牛頓流體假設，在源項考慮重力及表面張力。則質量和動量守恆方 程式可分別表示為： 0 V t (2.4)
V VV P V g f t (2.5) 其中
為密度，V 為速度，P為壓力，
為黏滯係數，g為重力加速度而 f 為表面張力。 因為此雙流體系統之流體被視為一均質流體進行分析，此均質流體具有連 續的速度分佈與溫度分佈。故系統的動量方程式可以用(2.5) 式表之，其 中密度及黏滯係數則以流體體積分率
進行混合如(2.6)和(2.7)式所示：
1 1 2 (2.6)
1 1 2 (2.7)其中下標 1、2 分別代表兩個不同流體，
是流體 1 所占的體積分率。在 單一流體的區域使用流體本身的流體性質，而介面附近則是使用體積分率 進行混合。所以整個計算範圍內所有性質是呈現片段連續(piecewise continuous)的分佈。 2.4 能量守恆方程式 考慮系統的能量方程式時，在低速流情況下，將忽略黏滯力所造成的 能量消耗項，而能量方程式可用流體內能型式表示如式(2.8)。因吾人希望 將兩相各自的能量方程式以類似處理動量方程式的方式進行合併。在擴散 項部份因為變數不一致無法直接將混和之流體性質提出，因此做了一些化 簡的動作。
v e k Ve e t c (2.8) 內能型式之能量方程式 如式(2.8) 所示，以內能型式表達能量方程式中的擴散項，再將兩流 體各自的能量方程式以體積分率進行混合：
1 2 1 1 2 2 1 1 2 2 1 2 1 2 (1 ) (1 ) (1 ) v v e e V e e t k k e e c c (2.9) 定義：12
1 1 2 2 1 2 1 1 e e e (2.10) 根據式(2.10) 整理式(2.9)，並且將擴散項化簡可得式(2.11)
1 2 1 2 1 2 1 2 1 2 (1 ) (1 ) v v v v e V e t k k e e c c k k e c c (2.11) 定義 1 2 1 (1 ) 2 v v v k k k c c c (2.12) 而內能型式之能量方程式可改寫成：
v k e V e e t c (2.13) 定義e1及e2如式(2.15), (2.16)
1 1 v 1 o e c T T (2.14)
2 2 v 2 o e c T T (2.15) 其中To為參考溫度，若將式(2.14)和(2.15) 代回式(2.10)，可得：
1 2 1 1 2 2 1 2 1 1 v o v o c T T c T T e (2.16) 因為本研究將雙流體系統視為一均質流體，故其溫度場為連續分佈，也就 是T T1 T2，因此可將式(2.16) 之溫度提出如式(2.17)。
1
1 1 2 2 1 2 2 1 1 o v v v v e T c c T c c (2.17) 求得 e 後，根據式(2.17) 可求得溫度分佈。 2.5 VOF 方程式 在本研究中使用 VOF(volume-of-fluid)方法，其體積分率定義如式(2.18) 1 Volume of fluid Volume of control volume (2.18) 1 , 1 0 , 2 0 1 , for fluid for fluid
for transition region
(2.19) 從[19]可知對雙流體流而言，體積分率的物質微分(material derivative)為 零： 0 D V Dt t
(2.20) 若將式 (2.20)在等號兩邊加上
V
，則可整理得到式(2.21)：
V V V t (2.21) 若考慮系統為不可壓縮流，可將式(2.21)改寫成：0
V
t
(2.22)14 2.6 相變化模型 在雙相流中考慮熱傳現象時，若兩相流體處於不同溫度時便會開始有 熱傳現象產生，當其一相達到此時系統壓力的飽和溫度時，並且界面處的 熱傳淨通量高於此時系統壓力的相變化潛熱時，便會在介面處開始產生相 變化現象。 相變化代表有部分質量的流體轉換流體性質，此時原本在介面就存在 的流體性質的不連續性需要進一步進行修正以滿足介面的跳躍條件 (jump conditions)，這些跳躍條件包括了質量、動量與能量，如式(2.23) 至 (2.25) [14]。
l ul uf n v uv uf n mf (2.23)
I f v l v l v l m u u n p p n n (2.24) f fg f v l v l T T m h q k k n n (2.25) 其中ul 、uv、uf 分別表示液相流體、氣相流體與介面速度，mf 為介面處 相變化之質量流率，而
、p、
依序表示剪應力、壓力與表面張力項，qf 表示介面之淨熱傳通量。 在本研究中假設介面溫度為系統壓力的飽和溫度Tf Tsat
P ，此假 設為一般液氣相變化之數值模擬所使用，適用於 macro scale 之沸騰問題 [20]。若是對於更小尺寸(如 micro scale)的沸騰問題，因為發現介面處存在一壓力的不連續性，也就是Pv Pl，故兩相之飽和溫度也不相等
sat v sat l T P T P 。此時需要對於介面溫度作進一步修正。但本研究之模 擬尺度屬於前者，因此採用介面溫度為系統壓力之飽和溫度之假設。 考慮質量的跳躍條件即式(2.23)，將式(2.25)代入式(2.23) 可將兩相間 的速度不連續性表示成式(2.26)。
f 1 1 v l fg v l q u u n h (2.26) 對一微小控制體積內的介面 F 上之微小截面面積進行積分，則介面處因相 變化產生的速度的不連續性可改寫成式(2.27)
1 1 1 f v l f F f f f fg v l F u x x u u n dA x x q dA h
(2.27) 其中 為一 delta 函數，xf 代表位於介面上的一點，而dAf 為介面 F 上的 微小面積。之後可將式(2.27) 加入原質量守恆方程式式(2.4) 作為修正， 如式(2.28)
1 1 1 f f f fg v l F V x x q dA t h
(2.28)16 2.7 表面張力
本研究中使用 CSF(continuum surface force)的表面張力模型。在介面 處因為表面張力的存而產生一壓力差可表示成式(2.29) i o P P P (2.29) 其中Pi表凹面側的壓力，Po則為击面側的壓力，
為表面張力係數，
為 平均曲率半徑。當兩流體相接即存在一表面張力係數
，若此兩流體互不 相融，
為正值；若相融則為負值。若
值為正值，表示流體 1 位於該介 面區段的击面側；而當
值為負值，表示流體 2 位於該區段的击面側，見 圖 2.2。 一從流體 2 指向流體 1 的法線向量可用體積分率的梯度表示： n (2.30) 而平均曲率半徑
可以改寫成單位法向量之散度： n n (2.31) 我們可將式(2.29)及(2.31)代入動量守恆方程式式(2.5)中的表面張力，可 得： n f n (2.32)2.8 邊界條件 (1) 進出口邊界條件 本研究中進出口採用壓力邊界，利用給定進出口不同壓力值形成壓力 差，此壓力差將會決定流量大小。而進出口速度計算採對流邊界條件，
0
cV
t
，其中
為傳輸性質而Vc 為對流速度。 (2) 牆邊界條件 本研究中速度場的牆邊界條件採用無滑移(no-slip)邊界條件，即是(u=0, v=o)。而溫度場部分：定溫條件T=constant或絕熱條件 T 0 n 18
第三章 數值方法
3.1 前言 本章將對於第二章所敘述的數學模式進行離散，所使用的離散方法為 有限體積法。首先介紹處理介面運動的 VOF 方程式，接著為動量方程式 之離散、連續方程式與動量方程式之耦合處理。再來介紹能量方程式之離 散以及介面溫度之處理，最後則為介面處熱傳量之計算以及解題流程。 3.2 VOF 方程式 本研究求解 VOF 方程式來預測介面的運動情況與位置分佈，使用由 崔燕勇及林仕文所提出的 CISIT 介面追蹤法，此方法的解題程序主要可 分為三步驟依序進行。首先為介面重建步驟，利用上一時階的 VOF 分佈 計算出介面位置；接著求解離散後的 VOF 方程式得到新時階的 VOF 分 佈；最後對於整個計算範圍的網格之 VOF 進行修正，避免因數值計算而 造成的不符合物理現象情況。本節將依序大略介求解 VOF 步驟。 3.2.1 介面重建 在本研究中使用
0.5作為為介面所在處，首先將原本存在計算網格 中心的 VOF 值，藉由式(3.1) 之線性權重關係式差分至網格節點上。 i i i v i i w w
(3.1)其中
v 表示網格節點的 VOF 值，下標 i 表示網格節點周圍的網格，而wi 則表示節點周圍的網格所用的權重因子，此權重因子為計算網格體積之倒 數。 藉由式(3.1) 可得到各網格節點的 VOF 值，若網格的某一邊上的兩節 點之 VOF，一者大於 0.5 又另一者小於 0.5，則表示有介面通過此邊，便 可利用線性內插差分出介面(
0.5)在邊上的位置，完成此網格的所有邊 的計算後，以相同程序完成所有計算格點的計算，可以得到一連續、片段 線性的介面分佈。 3.2.2 VOF 方程式之離散 將式(2.22) 對一控制體積進行體積分，再利用高斯散度將體積分化為 面積分如式(3.2)。 0 V S d V d S t
(3.2) 非穩態項(unsteady term)之離散：
n o
D D V d t t
(3.3) 其中
表示網格的體積，t表示時間步階大小，上標 n 表示新時階的值， 上標 o 則表示為前一時階的值。 對流項(convection term)之離散：20 (3.4) 所示： w f f f S V d S V S
(3.4) 其中下標 f 表示控制體積的表面處的性質，Sfw則是表示面上被浸潤的區 域的面向量，見圖 3.2。 經由式(3.3) 及式(3.4) 可將式(3.2) 改寫成下式(3.5)。
Dn Do
f f w 0 f V S t
(3.5) 3.2.3 VOF 修正 在理想情況中，當介面通過某網格後(介面不在此網格內)，此網格之 VOF 應該為 1 或 0，但實際上，某些情況下有些網格之 VOF 分佈並未符 合以上的假設，因此必須針對這些網格之 VOF 進行修正。一般而言，可 將這些情況整理為四種類型： 3.2.3.1 過度填充(Over filling, P 1) 參考圖 3.3，當一介面前端從主格點跨入下一個網格時，因為在上一 時階(介面以實線表示之)在主格點右側的面(編號 2)上之Sfw為 0，故此面 將不會有體積通量流出，如此可能造成主格點之 VOF 超過 1(P 1)，此 時稱這種情形為過度填充。故需將其修正至 1，修正的方法為依照體積通 量的權重函數如式(3.6)將多出來的 VOF 值分配到其下游的格點，其下游格點的 VOF 即是根據式(3.7) 進行修正，修正後使主格點的 VOF 為１。
max , 0 max , 0 fi Ni fj j C w C
(3.6)
1
P Ni Ni Ni P Ni v w v (3.7) 其中w 為下游格點之權重函數，Ni Cfi VfiSfi為通過 fi面之體積通量，Ni為 下游格點之 VOF，vP、vNi分別表主格點與下游網格體積。 3.2.3.2 過度枯竭(Over depleting, P0) 參考圖 3.4， 當介面尾端離開主格點 P 而進入下一格點時，若流出體 積通量過大造成主格點之 VOF 小於 0(P 0)時，稱此種情況為過度枯竭， 必須將主格點之 VOF 修正為 0。使用類似處理過度填充的方法，依照體 積通量的權重函數將其下游格點之 VOF 往主格點補，直到主格點之 VOF 為 0。式(3.6) 為所用之權重函數，而式(3.8) 則為下游格點之 VOF 修正 式。 P Ni Ni Ni P Ni v w v (3.8) 3.2.3.3 填充不足( Under filling, P 1) 在模擬 shear flow 時，因為在旋轉方向有極大的速度梯度，可能發生22 稱之為填充不足。此時主格點四個節點之 VOF 都大於 0.5，代表介面已 經通過主格點，但實際主格點的 VOF 卻小於 1。修正方法也是藉由權重 函數將下游網格的 VOF 往主格點補直到主格點之 VOF 為 1，式(3.6)為所 用的權重函數，並依照式(3.7)修正下游網格之 VOF。 3.2.3.4 枯竭不足(Under depleting, P0) 同樣在模擬 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 動量方程式之離散化 動量方程式可表示如式(3.9)
V
Q t (3.9) 其中 表示速度，Q為動量方程式中的源項，，若將此方程式進行體積分 並且配合高斯散度定理可得到式(3.10)。
S S d V d S d S Q d t
(3.10) 3.3.1 非穩態項(unsteady term) 將式(3.10)中的非穩態項進行離散：
o n o P P P d t t
(3.11) 其中下標 P 代表計算主格點，上標 n、o 各代表新時階及上一個時階的值。 3.3.2 對流項(convection term) 對於對流項的面積分可以表示如式(3.13)：
C f f f f f f f f S V d S V S F m
(3.12) 其中mf 為通過計算面的質量通量。而對流通量
f 則是以一通量限制方程 式所決定如式(3.13)，其細節請參考 3.2 節。 ( )r A D 24
其中下標 D 表 donor 格點，A 則為 acceptor 格點。而對於對流通量的 donor, acceptor 計算如式(3.15) 所示。 , 0 , 0 D P A nb f D nb A D f for m for m (3.14) 本研究之動量方程式採用的通量限制方程式為 Van Leer ，表示如式(3.15) ( ) 1 r r r r (3.15) 3.3.3 擴散項(diffusion term) 擴散項的面積分可化簡如式(3.16)
o
D f f f f f S d S S F
(3.16) 令：Sf d
Sf d
其中d 表一從主控容中心指向鄰近控容中心之向量，可用式(3.17)表示。見 圖 3.7 。 2 f d d f S d S (3.17) 而擴散通量即可改寫成：
2 o f f D o o o o f nb P f f f d f S F S d S (3.18) 3.3.4 源項(source term) 源項的體積分可表示為：
P Q d q
(3.19) 其中動量方程式中的源項包括壓力、重力以及表面張力，敘述如下。 壓力項 壓力項的面積分之化簡近似： o o f f f S Pd S
P S P
(3.20) 重力項 重力項的體積分之化簡近似： o gd g
(3.21) 表面張力項 在處理表面張力項時，以式(3.22)表示主格點上體積分率之梯度，而 計算面上的體積分率
f 則是由相鄰兩個格點的體積分率進行內插所得， 最後表面張力項以式(3.23) 表示。
1 f f p f S
(3.22)
f f f f f f f d S S
(3.23) 3.3.5 線性代數式之整理 將上述動量方程式之各項離散式整理為式(3.24) 的型式。 n n o P P nb nb A
A
Q P (3.24)26 P nb P nb A A t
(3.25)
2 max , 0 o f f nb f d f S A m S
(3.26)
2 o o o o o o o f A D f f f P P P P f r Q m S d g t
(3.27) 3.4 連續及動量方程式耦合處理 若單純求解動量方程式，並無法確保求得之結果一定滿足連續方程式。 所以本研究將連續與動量方程式耦合處理，其主要包括一次壓力預測步驟 以及兩次壓力修正步驟，最主要藉由兩次修正步驟可以使計算結果滿足連 續方程式，而將連續方程式與速度場耦合在一起。此演算法計算程序介紹 如下。 預測步驟 首先將式(3.24) 改寫為：
o
P P nb nb P nb A V
A V S P (3.28) 其中VP表未經修正的速度項並不滿足連續方程式，S 表源項。將上式同除 以AP，並且 令 nb nb nb P p A V S H A
(3.29)可得： nb nb o o nb P P P P p P P P P A V S V P H P A A A
(3.30) 其中式(3.30)中的上標“ ”代表此值是根據一權重因子
f 將主格點 P 與相 鄰格點 C 進行內插所得。 第一次修正步驟 在進行第一次修正時，會考慮到質量守恆而得到一個更新的壓力，而 在第一次修正的所得的新速度場以上標**表示，而新壓力場以上標*表 之。
** * P P nb nb P nb A V
A V S P (3.31) 將式(3.31)同除以AP並且使用式(3.29)改寫成： ** * * nb nb nb P P P P p P P P P A V S V P H P A A A
(3.32) 而定義 ** * * , o VV V PP P 可以得到主格點的第一次速度修正方程式： P P P P V P A (3.33) 而計算面上的速度修正方程式也以同樣方式可得： f f P f V P A (3.34) 所以體積通量的修正方程式則表示成： 28 令Sf d
Sf d
將式(3.35) 改寫成式(3.36)：
2 f f f Pnb f f P f Pnb f P f S P P S d A S A (3.36) 其中 2 f nb nb f S d S 而式(3.36)中 Pf Pnb 改寫成Pnb PP，可得：
2 f f nb P f f P f Pnb f P f S P P P S d A S A (3.37) 所以第一次體積通量修正方程式可表示如：
2 ** * f f f nb P f f P f Pnb f P f S P P P S d A S A (3.38) 其中的 * f 項可以從以下步驟得到，首先將 * f 以式(3.32)的型式表示： * * o f f f P f V H P A (3.39) 其中 * * o f f f P f H V P A (3.40) 而 * f V 與Pfo則是以一權重因子f對主格點及相鄰點差分所得，如：
1
o o o f f nb f P P P P (3.41)
* * * 1 f f nb f P V V V (3.42) 將式(3.40)代回式(3.39)可得： * * o f f f f P f P f V V P P A A (3.43)其中 P f A 是採主格點與相鄰格點的幾何平均： 1 2 P f P P P nb A A A (3.44) 所以最後將計算面上的體積通量表示如下：
* * * * * * o f f f f f f f f P f o f f f f P f V S V S P P S A V S P P d A (3.45) 而根據離散後的連續方程式可知： ** * 0 f f f f f f
(3.46) 若如 2.6 節考慮相變化現象時，連續方程式將修正為式(2.40)，此時離散 後的連續方程式為式(3.47) ** * f f 1 1 f f f f f f fg v l q A h
(3.47) 其中qf 為介面處淨熱傳量，Af 為介面截面積 若將式(3.38)代入式(3.46)或(3.47)，再整理成如式(3.48)的線性代數式 * 1 1 1 2 P P P P nb nb f P P f f A P
A P
S S (3.48) 其中 1 1 P S , 1 2 P S 表示第一次修正之源項，而各項為 P P P nb f A
A (3.49) 2 f P nb P nb f S A A S (3.50)30
1 *
1
for phase change
1 1 f f P f f fg v l q A S h
(3.51)
1 2 P f f f P f S P S d A
(3.52) 第二次修正步驟 接下來進行第二次修正步驟，本次修正所得到的新速度場以上標“***” 表之，新壓力場則以上標“**”表之，而動量方程式可以表示如式(3.53)
*** ** ** P P nb nb P nb A V
A V S P (3.53) 以類似前一步驟的處理方式可將第二次速度修正方程式表示如： nb nb f P f P P f f A V V P A A
(3.54) 其中定義 *** ** ** * ** * , , VV V VV V PP P 之後得到新的體積通量修正方程式： nb nb f f f f P P f f A V P S A A
(3.55) 根據連續方程式，可將式(3.53)整理為式(3.56)的線性代數式型式： * 2 2 1 2 P P P P nb nb f P P f f A P
A P
S S (3.56) 其中 2 2 1, 2 P P S S 表示第二次修正的源項 P P P nb f A
A (3.57)2 f P nb P f Pnb f S A A S (3.58) 2 1 nb nb f P f P f A V S S A
(3.59)
2 2 P f f f P f S P S d A
(3.60) 3.5 求解動量方程式所用之流體性質 CISIT 的介面追蹤法可以將介面位置控制在一個計算網格內，但若直 接使用 VOF 作為單一流體模型之混和係數計算時，在介面處會產生流體 性質極大的不連續性(如質量)，如此易導致數值計算之不穩定。為了改善 這樣的情況，本研究利用式(3.1) 先將網格中心的 VOF 差分至網格節點上， 再將各節點的 VOF 作平均後差分回網格中心。以上稱為一次平滑過程， 經過平滑過程後介面將會散開而較為平滑，使得 VOF 不像原先類似步階 函數分佈。 將平滑處理後的 VOF 代入式(2.6) 及式(2.7)進行質量與黏滯係數的 計算，動量方程式將使用平滑後的質量與黏滯性質以改善介面處的性質不 連續性及計算之不穩定性。32 3.6 能量方程式之離散化 3.6.1 前言 關於雙流體流，原本不同流體各自遵守各自的能量方程式，而各自解 出的溫度分佈在兩流體的介面處，必須滿足熱傳量的跳越條件(jump condition)。但本研究將兩流體混合後視為一均質流體，此均質流體的速度 場及溫度場皆為連續分佈。整個雙流體流系統使用同一組動量及能量方程 式求解，所用的流體性質係依照流體之體積分率混合兩相異流體的流體性 質。 3.6.2 單一均質模型 如前言所提，本研究採用為單一均質流體進行計算，而在 2.4 節能量 守恆方程式中，能量方程式可表示為內能型式：
e v k e V e e Q t c (3.61) 其中e表示質量平均內能，而各流體性質均在 2.4 節中已定義。，Qe為源 項，若將式(3.61)進行體積分並且配合高斯散度定理可得式(3.62)
e v S S e k d Ve d S e d S Q d t c
(3.62) 之後對式(3.62)進行離散步驟，此離散過程可參考 3.4 節的動量方程式之 離散化，對流項之處理同樣採用一通量限制方程式搭配 donor-acceptor 方法，最後整理為一線性代數式。 3.6.3 線性代數式之整理 將上述能量方程式之各項離散式整理為式(3.63)的線性代數型式。 n n P P nb nb e nb A e
A e Q (3.63) P nb P nb A A t
(3.64)
2 max , 0 o f nb f v f d f s k A m c s (3.65)
2 o o o o o o e f A D f f P P f v f r k Q m e e e S d e c t
(3.66) 3.6.4 介面處之網格溫度處理 如 2.6 節中所提本研究中假設介面溫度為系統壓力之飽和溫度，為了 在數值計算時維持此假設，進行以下處理。當蒐集完全部計算網格之能量 方程式各項係數後，若某一網格為介面所在之網格將會設定其 30 1. 10 P A ， 各Anb 0及Qe 0，藉由以上處理將介面所在網格之溫度維持於飽和溫 度。 3.7 介面處熱傳量之計算 在 2.6 節中考慮相變化時，需要求得介面處的熱傳量34 為介面溫度、Tn1、Tn2 各為位於兩相之溫度參考點，n表某一定長。在 本研究中n通常設定為一至二倍網格長度。 計算 VOF 方程式後可得到介面位置，在主計算網格內介面截線的中 點上，往其法線方向各取一定長n得到位於介面兩側之兩相異點。接著 各以此兩點為圓心，n為半徑作圓，見圖 3.6。若相鄰網格中心位於此圓 內，則選擇此網格作為計算此兩點溫度之參考網格。最後以相鄰網格中心 至此兩點的距離之倒數作為權重函數，計算出此兩點的溫度Tn1及Tn2，以 1 n T 為例，如式(3.67)及圖 3.7 所示。 1 1 i i i n i i T l T l
(3.67) 其中下標 i 表示計算溫度參考點時所用的相鄰網格，l 表相鄰網格中心至 參考點之距離。 在計算介面處的熱傳量後，可進一步求出介面產生相變化的體積流率 1 f f v fg v q A h ，其為因相變化現象造成的氣體體積增加量，因此需要依 此變化量對於此時階的 VOF 進行修正。3.8 解題流程 在本章已經介紹 VOF 方程式、動量方程式與能量方程式之數值方法， 而雙流體流系統的解題流程如下： 1. 給定所有變數的初始值 2. 利用上個時階的體積通量求解 VOF 方程式得到新的體積分率
與重建 界面位置 3. 使用新的介面位置及上個時階所得溫度分佈求得界面處之熱傳量供相 變化計算使用，並且根據相變化造成的氣體體積變化修正此時階的 VOF 值 4. 求解動量方程式並且進行壓力修正以滿足連續方程式 5. 利用步驟 4 所得新的速度分佈求解能量方程式 6. 利用新的能量分佈求得溫度場分佈 7. 重複步驟 2 至步驟 6 至設定總計算時階數36
第四章 結果與討論
4.1 相變化測試(一)，Stefan’s 1st problem 為了測試本研究所使用的相變化模型的準確性，將以兩個具有理論解 的相變化作為測試。首先為由物理學家Jožef Stefan 在 1890 年所提出的 Stefan 問題。圖 4.1 為其示意圖，考慮一維分佈情況的氣態及液態流體分 別位於氣液介面左側與右側，左側氣相流體的左端為一定溫(T Twall)的牆 邊界條件；而其右端則為定溫(Tf Tsat)的氣液介面。由於氣相流體溫度高 過此系統的飽和溫度，所以介面處會產生相變化，使得氣相流體增加，而 介面開始向液相側移動。 在氣相流體區域，假設流體速度為零，並且溫度分佈須滿足一維熱擴 散方程式如式(4.1)；而在液相區域，假設溫度保持飽和溫度。而式(4.2)與 式(4.3)為此方程式的邊界條件，式(4.4)則為介面處的能量跳躍條件。 2 2 T T t x , 0 x
t (4.1)
,
sat T x t t T (4.2)
0,
wall T x t T (4.3) f v f fg x t T q u h k x (4.4) 其中uf 為介面移動速度 根據[21]，[22]此問題之 Neumann 解如式(4.5)及式(4.6)所示
t 2 vt (4.5)
,
2 sat wall wall v T T x T x t T erf erf t (4.6)
2
exp v wall sat
fg c T T erf h (4.7) 其中 t 為介面位置，為一滿足式(4.7)的根，erf 表 error function
本例子所用的流體為一大氣壓之飽和狀態純水，其流體性質列於表一。飽 和溫度為 100o C而牆溫度為 125oC ，由式(4.7)可解得=0.105612。 本例子採二維之 10mmx1mm 的均勻計算網格，共有 50x5、100x10 及 200x20 三組網格大小，所用的計算時間步階為 3 2 10 s。邊界條件設定 在左側 x=0mm 處為保持 125o C的定溫條件；右側 x=10mm 處則為 100oC 的定溫條件；至於 y=0mm 及 y=1mm 處則為絕熱邊界條件。初始的介面 位置定在 x=0.3mm 處，並且以此處作為介面移動原點。在計算過程中， 介面所在格點以及液相格點之溫度均為飽和溫度。 本研究求解能量場使用的是內能型式能量方程式，此型式能量方程式 的熱擴散係數為 v k c 。在本例子使用了三種不同型式熱擴散混合係數進行計 算，分別為調和平均(harmonic mean)、混合型式(mixed)及氣相熱擴散係數， 依序以式(4.8)至式(4.10)表示，而介面熱傳計算所用 n x。 v v v l k k c c k
