基於流體體積法之兩相流(包含熱質傳)數值方法及其應用

國 立 交 通 大 學

機 械 工 程 學 系

博士論文

基於流體體積法之兩相流(包含熱質傳)

數值方法及其應用

VOF-Based Algorithms and Their Applications for

Two-Phase Flows (Including Heat and Mass Transfers)

研 究 生：林 仕 文

指 導 教 授：崔 燕 勇 教授

基於流體體積法之兩相流(包含熱質傳)數值方法及其應用

VOF-Based Algorithms and Their Applications for

Two-Phase Flows (Including Heat and Mass Transfers)

研 究 生：林仕文 Student:Shi-Wen Lin

指導教授：崔燕勇 Advisor:Yeng-Yung Tsui

國立交通大學

機械工程學系

博士論文

A Thesis

Submitted to Department of Mechanical Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Doctor of Philosophy

In

Mechanical Engineering

January 2014

Hsinchu, Taiwan, Republic of China

基於流體體積法之兩相流(包含熱質傳)數值方法及其應用

學生：林 仕 文 指導教授：崔 燕 勇 博士

國立交通大學機械工程學系

摘 要

本研究最主要的目標是發展一套可用於非結構性網格中的含相變化的 兩相流數值方法。基於流體體積法(VOF)，本文提出兩套相關的計算方法。 第一種方法稱為通量混合介面捕捉法(FBICS)，其是透過直接求解流體體積 分率的傳輸方程式來捕捉介面的運動。在計算的過程中，為了同時維持介 面的鮮明度及體積分率的界限性，本研究是使用通量混合的方式來處理網 格控制面上所需要的對流通量，然而本方法的缺點在於其所獲得的介面通 常佔據了多個網格的寬度。另外一套方法則稱為守衡內插介面追蹤法 (CISIT)，本方法是先以內差的方式重建介面後，再以先預測後修正的方式 處理介面的移動。此方法計算後所得的流體體積分率分布，除了介面所在 的網格外，其餘皆是 0 或 1 的均勻分布，而且介面只有佔據一個網格的寬 度。不同於 PLIC 法，由於 CISIT 法的方法相當簡單，因此可輕易地在幾何 外觀較為複雜的非結構性網格中使用，而且不需任何複雜的處理便可推展 至三維的問題之中。在 PLIC 法中，其介面重建的方式並非簡單容易處理的， 而且在計算面上體積通量時必須考慮許多類型的介面形狀(二維流場中有 16 種情形而三維則有 64 種)來進行，因此也造成其介面推移的計算上相當 複雜。經由許多問題的測試，可知本研究的兩套數值方法於自由表面流中 模擬後所得的結果與理論解或實驗數據都相當吻合並且可在非結構性的網 格中所使用。

接著為了能夠模擬含相變化的兩相流問題，本研究特別將相變化所產生 的質熱傳效應導入守衡內插介面追蹤法之中。經由介面處質量及能量的跳 躍邊界條件可獲得通過介面的質傳量並將此質傳量以源項的方式導入連續 方程式之中。然後，對於溫度場而言本研究是將介面視為一組能量內邊界 的方式來處理，並且最後將以全隱性的求解方式來計算能量方程式。此外， 基於假設介面處溫度及熱傳量連續性的條件來取代跳躍邊界，本方法亦可 推展至模擬無相變化的雙流體熱傳問題。將本方法應用於鄰近臨界壓力的 水平平板薄膜沸騰問題中，其結果顯示在不同的壁過熱溫度下將有不同的 沸騰模式出現。根據不同的壁面過熱度，沸騰的模式主要可分為五種：單 氣 泡 模 式 ( T 15C ) 、 單 / 多 氣 泡 模 式 ( T 15 ~ 19C ) 、 單 噴 流 模 式 ( T 19 ~ 22C)、雙氣泡模式( T 22 ~ 27C)及雙噴流模式( T 27C)。在 單氣泡模式中，模擬所得的時均化 Nusselt 數與半經驗式之結果相當一致。 再者，透過模擬水平圓管的薄膜沸騰問題可以證明本文所提出的含相變化 之兩相流數值方法可使用於複雜幾何外觀的沸騰流動中。 本文最後將含相變化的守衡內插介面追蹤法修改為可用於計算三維兩 相流的流場之中。不同於二維的流場，在網格中本研究是使用多個不共面 的三角形介面來進行三維介面的重建。本方法首先透過許多不含相變化的 雙流體問題來驗證其在三維流場中對於介面追蹤及預測的能力。另外，也 將本方法應用至三維水平平板的薄膜沸騰流動之中。結果顯示利用本方法 模擬所獲得的時均化 Nusselt 數與 Klimenko 所提出的半經驗公式之結果相 當吻合，尤其是在壁過熱度為 10℃的情況下。最後，無相變化的熱傳模型 也用於模擬油槽內的熔解錫液滴及水中的高溫辛烷噴流之問題中。

VOF-Based Algorithms and Their Applications for

Two-Phase Flows (Including Heat and Mass Transfers)

Student: Shi-Wen Lin Advisor: Dr. Yeng-Yung Tsui

Department of Mechanical Engineering

National Chiao Tung University

ABSTRACT

This paper is aimed at developing a numerical method for two-phase flows with phase change on unstructured grids. In this article, two schemes are presented based on VOF (volume-of-fluid) method. The first scheme is to capture the interface by solving the advection equation of the volume fraction directly, termed as FBICS. In order to maintain the sharpness and boundedness of the interface, the convective flux through each cell face is determined by means of flux blending. The weakness of this method is that the interface region will occupy several grid spaces. In the other scheme (termed as CISIT), the interface is reconstructed first using interpolation practice, following by a predicted-correction procedure to handle the movement of the interface. Except for the interface cells, the VOF distribution is uniform, either in 1 or 0, and the interface occupies only one cell in its width. Unlike the PLIC method, the CISIT can be easily extended to unstructured grids with arbitrary geometry and 3-D problems without causing any further complication because its formulation is very simple. In PLIC, the reconstruction of the interface is not straightforward and the procedure to advance the interface is complicated because a large number of interface configurations (16 configurations for 2-D flows and 64 for 3-D flows) must be considered for determining of the flux across cell faces. Tests on a number of cases reveal that results via these two schemes in this study, which can be used on the unstructured grids, give good agreement with exact solutions or experimental data of free surface flows.

In order to simulate the two-phase flow with phase change, the CISIT method is extended to include heat and mass transfer due to phase change. The mass transfer across the phase boundary is determined by taking into account the mass and energy jump conditions at the interface and added as a source term in the continuity equation. Then, the interface is treated as an internal boundary condition in the temperature flied. Finally, the energy equation is solved in an implicit way. Besides, this method is also extended to simulate the heat transfer of two-fluid flows without phase change based on the assumption of the continuity conditions of the temperature and heat flux instead of jump conditions at the interface. Application to film boiling flow on a horizontal plate at a state near the critical pressure shows that the boiling mode will be different at different superheat temperatures. According to different superheat tempera- tures, the boiling flows can be divided into five modes: single-bubble mode ( T 15C ), single/multiple bubble mode (  T 15 ~ 19C ), single-jet mode

( T 19 ~ 22C), double-bubble mode ( T 22 ~ 27C), and double-jet mode

( T 27C). In the single-bubble mode, good agreement with semi-empirical

correlations was obtained in terms of averaged Nusselt number. Furthermore, simulation of film boiling flow on a cylinder demonstrates that this method is applicable to boiling flow with complex geometry.

Finally, the CISIT method with phase change is modified to calculate three-dimensional two-phase flows. Unlike two-dimensional flow, the interface is reconstructed with several non-coplanar triangular interfaces within the grid. First, this method was tested through computations of a number of two-fluid flows without phase change to validate the capability of tracking the interface in three- dimensional flows. In addition, this method was also applied to simulated film boiling flow on a horizontal plate. It can be shown that the space and time averaged Nusselt numbers obtained from the current simulations have good agreement with the semi-empirical correlations of Klimenko, especially for ₁₀0

T C

  .Finally, the heat transfer model without phase change was used to simulate the molten tin droplet in oil and the octane inlet in water.

誌 謝

首先最需要感謝的當然是指導我七年半的崔燕勇教授，自從碩士班來到 交通大學後，很慶幸地遇到了崔老師，在他耐心且用心的指導下才有今日 如此紮實的研究內容以及成果。崔老師是我見過做研究最積極且指導學生 不遺餘力的教授，而且其對於研究的態度精確且踏實是老師給我最大的影 響，將來不管是在業界或是學界相信這種研究態度會讓我受益良多。老師 對學生的好實在不是三言兩語在此可以言明，我心中除了感謝還是感謝， 非常感謝老師您這一路來的指導與付出。 接著當然是要感謝我的父母，謝謝你們在這一條漫長的路上給我最大的 支持及陪伴，如果沒有你們在我背後默默的付出，今天肯定無法到達這條 路的終點，再未來的日子裡我一定會努力的回報妳們這幾十年來的養育之 恩。接著還要感謝的就是這五年來陪我ㄧ起走過的女友，在我為了研究繁 忙的日子裡，你是我支撐下去的力量，謝謝你的包容及等待。 最後要感謝所以交大機械熱流組的教授們及學長學弟們，這一路下來感 謝大家的幫忙及照顧，沒有你們相信我在交大的日子不會這麼精彩快樂， 還有要感謝交大資管所及光電系的各位棒球隊的戰友們，讓我在閒暇的時 間能夠有棒球可以打，非常謝謝你們，希望以後黑龍旗有機會再跟你們一 同出戰。

目 錄

摘 要...i ABSTRACT ...iii 誌 謝...v 目 錄...vi 表目錄...ix 圖目錄...x 符號說明...xv 第一章 緒論...1 1.1 簡介...1 1.1.1 兩相流簡介...1 1.1.2 薄膜沸騰簡介...5 1.2 文獻回顧... 11 1.2.1 兩相流之數值方法... 11

1.2.1.1 移動網格法 (Moving Grid Method)... 11

1.2.1.2 前端追蹤法 (Front-Tracking Method)... 11

1.2.1.3 等位函數法 (Level-set method)...12

1.2.1.4 體積追蹤法 (Volume Tracking Method)...13

1.2.1.5 流體體積法 (Volume-of-Fluid Method, VOF)...14

1.2.1.5.1 介面捕捉法 (Interface Capturing Method)...15

1.2.1.5.2 介面重建法 (Interface Reconstruction Method)...17

1.2.1.6 流體體積法與等位函數法之混合...20 1.2.2 薄膜沸騰之數值模擬...21 1.2.2.1 前端追蹤法...21 1.2.2.2 等位函數法...22 1.2.2.3 流體體積法...23 1.3 研究目的...24 第二章 數學模型...34 2.1 簡介...34 2.2 流體體積(VOF)方程式...35 2.3 統御方程式...36 2.4 表面張力模型...38 2.5 相變化之熱質傳模型...39 2.5 無相變化之熱傳模型...41 2.6 邊界條件...42 第三章 通量混合介面捕捉法(FBICS)...45 3.1 簡介...45

3.2 流體體積方程式之離散...46 3.3 通量限制函數...47 3.4 通量混合介面捕捉法...48 第四章 守恆內差介面追蹤法(CISIT)...59 4.1 簡介...59 4.2 介面重建...61 4.3 流體體積分率預測步驟...61 4.4 流體體積分率修正步驟...63 4.4.1 填充過度(over-filling, _P  ) ...64 1 4.4.2 耗竭過度(over-depleting,_P  ) ...65 0 4.4.3 填充不足(under-filling,_P  ) ...65 1 4.4.4 耗竭不足(under-depleting,_P  ) ...66 0 4.5 CITSIT 計算流程...66 4.6 流體體積分率平滑化...67 第五章 速度場之數值方法...73 5.1 簡介...73 5.2 動量方程式的離散...73 5.2.1 非穩態項...74 5.2.2 對流項...74 5.2.3 擴散項...75 5.2.4 源項...76 5.2.5 動量方程式之代數方程式...78 5.3 速度與壓力之耦合...79 5.3.1 預測步驟 (predictor step)...79

5.3.2 第一次修正步驟 (first corrector step)...80

5.3.3 第二次修正步驟 (second corrector step)...83

5.3.4 PISO 演算法的計算流程...85 5.4 壓力出口邊界之處理...85 第六章 熱傳及質傳之數值方法...88 6.1 簡介...88 6.2 相變化模型之熱傳數值方法...89 6.2.1 介面熱通量之計算...89 6.2.2 能量內邊界之處理...89 6.3 無相變化之熱傳數值方法...90 6.4 能量方程式之離散...91 6.5 兩相流模型之計算流程...93 第七章 FBICS 法之驗證與自由表面流之應用...97 7.1 簡介...97

7.2 介面於均勻速度場之傳輸...98 7.3 介面於剪切流(Shear Flow)中之拉伸...102 7.4 二維壩體潰堤...104 7.5 壩體潰堤流經阻塊...107 7.6 水力湧潮 (Hydraulic Bore)...109 第八章 CISIT 法之驗證與雙流體流場之應用...141 8.1 簡介...141 8.2 介面於均勻速度場中之傳輸...142 8.3 介面於剪切流中之拉伸...146 8.4 二維壩體潰堤...148

8.5 單一上升氣泡 (Single Rising Bubble)...150

8.6 雷利-泰勒不穩定性之問題...153 第九章 CISIT 法於二維薄膜沸騰分析...188 9.1 簡介...188 9.2 一維汽化問題...189 9.3 水平平板薄膜沸騰...192 9.4 水平圓管薄膜沸騰...201 9.5 多模式(multi-mode)長平板薄膜沸騰...206 第十章 CISIT 法之三維應用...236 10.1 簡介...236 10.2 CISIT 法之三維化處理...238 10.2.1 三維浸潤面積之計算...238 10.2.2 三維介面類型及重建...239 10.2.3 三維介面處熱通量計算及能量內邊界...240 10.2.4 平行化運算...242 10.3 三維壩體潰堤...244 10.4 氣泡運動分析...246 10.4.1 單一上升氣泡...246 10.4.2 氣泡融合...251 10.5 液滴撞擊薄層液面...254 10.6 三維水平平板薄膜沸騰...256 10.7 無相變化之熱傳問題...260 第十一章 結論...294 第十二章 參考文獻...298 簡 歷...308 論文發表...309

表目錄

表一 介面傳輸之L 誤差與網格收斂速度表(四邊形網格1 u1 v0) ... 155 表二 介面傳輸之L 誤差與網格收斂速度表(四邊形網格1 u1 v1) ... 155 表三 介面傳輸之L 誤差與網格收斂速度表(三角形網格1 u1 v1) ... 156 表四 介面於剪切流中之L 誤差與網格收斂速度表 ... 156 1 表五 雷利-泰勒不穩定性之之L 誤差與網格收斂速度表 ... 156 1 表六 二維平板薄膜沸騰於不同網格之氣泡生成周期及時均化 Nusselt 數 ... 210

表七 過熱氣體性質表(P_sat 218bar，Tsat 3730C) ... 210

表八 二維平板薄膜沸騰之時均化 Nusselt 數與半經驗式對照表 ... 210 表九 二維水平圓管薄膜沸騰之時均化 Nusselt 數與半經驗式對照表 ... 210 表十 單執行緒與多執行緒平行計算之 CPU 運算時間 ... 264 表十一 單一上升氣泡之終端雷諾數表(*：2D 模擬結果) ... 264 表十二 單一上升氣泡底部與頂端沿 z 方向上的壓力變化率 ... 264 表十三 三維平板薄膜沸騰之時均化 Nusselt 數與半經驗式對照表 ... 264

圖目錄

圖 1.1 兩相流的基本流動形態 ... 26 圖 1.2 兩相流模擬所使用的網格種類 ... 27 圖 1.3 水沸騰曲線 ... 28 圖 1.4 Berenson 假設之薄膜模型 ... 28 圖 1.5 前端追蹤法示意圖 ... 29 圖 1.6 等位函數法示意圖 ... 29 圖 1.7 MAC 法示意圖 ... 30 圖 1.8 流體體積(VOF)法示意圖 ... 30 圖 1.9 SLIC 法沿 X 方向掃視之介面重建示意圖 ... 31 圖 1.10 PLIC 法簡化後介面於網格內分布情形(介面方向與 x 軸之夾角為0 ~ / 2 ). 32 圖 1.11 不同 VOF 之介面重建法比較示意圖 ... 33 圖 2.1 介面曲率大小之示意圖 ... 44 圖 3.1 任意形狀之控制體積與鄰近網格示意圖 ... 54 圖 3.2 網格控制面之上、下游和上上游之示意圖 ... 54 圖 3.3 控制面的上上游格點於非結構性網格中之示意圖 ... 55 圖 3.4 正規化視圖之 CBC 和 TVD 界限準則 ... 56 圖 3.5 FBICS-A 法所使用的非線性通量限制函數 ... 57 圖 3.6 FBICS-B 法所使用的非線性通量限制函數 ... 58 圖 4.1 CISIT 法之介面重建示意圖 ... 69 圖 4.2 介面橫跨不規則網格及所構成的浸潤區域之示意圖(紅色為流體 1) ... 69 圖 4.3 流動前端介面跨越網格面往鄰近網格推進之示意圖 ... 70 圖 4.4 流動尾端介面跨越網格面離開原網格之示意圖 ... 71 圖 4.5 使用兩次平均平滑化處理後的體積分率分布示意圖 ... 72 圖 5.1 非結構性網格中主格點與鄰近格點及面向量(S_f)之示意圖 ... 87 圖 5.2 邊界壓力計算及出口壓力邊界處理方式之示意圖 ... 87 圖 6.1 介面處熱通量計算之示意圖 ... 95 圖 6.2 二維介面處能量內邊界處理之示意圖 ... 95 圖 6.3 兩相流計算流程圖 ... 96 圖 7.1 介面於均勻速度場中所使用的三角網格 ... 111 圖 7.2 介面於均勻速度場中傳輸之理論解 ... 112 圖 7.3 介面於均勻速度場中傳輸之結果 (C_N 0.75，四邊形網格) ... 113 圖 7.4 介面於均勻速度場中傳輸之結果 (C_N 0.5，四邊形網格) ... 114 圖 7.5 介面於均勻速度場中傳輸之結果 (C_N 0.1，四邊形網格) ... 115 圖 7.6 介面於均勻速度場中傳輸之結果 (C_N 0.75，三角形網格) ... 116 圖 7.7 介面傳輸之數值誤差與庫倫數關係圖 ... 117

圖 7.8 剪切流中初始介面與速度場分布之示意圖 ... 118 圖 7.9 剪切流中介面正轉與反轉之結果(C_N 0.75，N=2) ... 119 圖 7.10 剪切流中介面正轉與反轉之結果(C_N 0.75，N=4) ... 120 圖 7.11 剪切流中介面正轉與反轉之結果(C_N 0.75，N=8) ... 121 圖 7.12 剪切流中介面正轉與反轉之結果(C_N 0.75，N=16) ... 122 圖 7.13 剪切流中介面正轉與反轉之結果(C_N 0.25，N=2) ... 123 圖 7.14 剪切流中介面正轉與反轉之結果(C_N 0.25，N=4) ... 124 圖 7.15 剪切流中介面正轉與反轉之結果(C_N 0.25，N=8) ... 125 圖 7.16 剪切流中介面正轉與反轉之結果(C_N 0.25，N=16) ... 126 圖 7.17 介面於剪切流中之數值誤差與庫倫數關係圖 ... 127 圖 7.18 二維壩體潰堤之幾何尺寸示意圖 ... 128 圖 7.19 二維壩體潰堤之實驗結果[74] ... 129 圖 7.20 二維壩體潰堤模擬所使用的三角網格 ... 130 圖 7.21 二維壩體潰堤隨時間之演進圖 (網格 80x50) ... 131 圖 7.22 二維壩體潰堤隨時間之演進圖 (網格 120x70) ... 132 圖 7.23 二維壩體潰堤隨時間之演進圖 (網格 6218) ... 133 圖 7.24 二維壩體潰堤隨時間之演進圖 (網格 9698) ... 134 圖 7.25 底部水流前端與時間無因次化關係圖 ... 135 圖 7.26 左側水柱高度與時間無因次化關係圖 ... 135 圖 7.27 壩體潰堤流經阻塊之實驗結果[74] ... 136 圖 7.28 壩體潰堤流經阻塊隨時間之演進圖 ... 137 圖 7.29 水力湧潮物理模型之示意圖 ... 138 圖 7.30 水力湧潮隨時間之演進圖 ... 139 圖 7.31 水力湧潮右側水位高度與理論解 ... 140 圖 8.1 介面於均勻速度場中所使用的三角網格 ... 157 圖 8.2 介面於均勻速度場中傳輸之結果 (u1,v0,C_N 0.75，四邊形網格) .... 158 圖 8.3 介面於均勻速度場中傳輸之結果 (u1,v0,C_N 0.5，四邊形網格) ... 159 圖 8.4 介面於均勻速度場中傳輸之結果 (u1,v0,C_N 0.25，四邊形網格) .... 160 圖 8.5 介面於均勻速度場中傳輸之結果 (u1,v0,C_N 0.1，四邊形網格) ... 161 圖 8.6 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.75，四邊形網格) .... 162 圖 8.7 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.5，四邊形網格) ... 163 圖 8.8 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.25，四邊形網格) .... 164 圖 8.9 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.1，四邊形網格) ... 165 圖 8.10 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.75，三角形網格) ... 166 圖 8.11 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.5，三角形網格) .... 167 圖 8.12 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.25，三角形網格) ... 168 圖 8.13 介面於均勻速度場中傳輸之結果 (u1,v1,C_N 0.1，三角形網格) .... 169 圖 8.14 介面傳輸之數值誤差與網格數目之關係圖(u1,v ，四邊形網格) .... 170 0

圖 8.15 介面傳輸之數值誤差與網格數目之關係圖 (u1,v ，四邊形網格) .... 171 1 圖 8.16 介面傳輸之數值誤差與網格數目之關係圖 (u1,v ，三角形網格) .... 172 1 圖 8.17 介面於剪切流中拉伸之結果(四邊形網格) ... 173 圖 8.18 介面於剪切流中拉伸之結果(三角形網格) ... 174 圖 8.19 介面於剪切流中誤差與網格數目之關係圖 (四邊形網格) ... 175 圖 8.20 介面於剪切流中誤差與網格數目之關係圖 (三角形網格) ... 175 圖 8.21 二維壩體潰堤隨時間之演進圖 (網格 60x35，平滑化 1 次) ... 176 圖 8.22 二維壩體潰堤隨時間之演進圖 (網格 60x35，平滑化 2 次) ... 177 圖 8.23 二維壩體潰堤隨時間之演進圖 (網格 120x70，平滑化 1 次) ... 178 圖 8.24 二維壩體潰堤隨時間之演進圖 (網格 120x70，平滑化 2 次) ... 179 圖 8.25 底部水流前端與時間無因次化關係圖 ... 180 圖 8.26 左側水柱高度與時間無因次化關係圖 ... 180 圖 8.27 氣泡形狀於靜止流體中上升之無因次參數關係圖[77] ... 181 圖 8.28 上升氣泡隨時間之介面與速度場演進圖(Eo1 , Mo0.001) ... 182 圖 8.29 上升氣泡隨時間之介面與速度場演進圖(Eo10 , Mo ) ... 182 1 圖 8.30 上升氣泡隨時間之介面與速度場演進圖(Eo100 , Mo1000)... 183 圖 8.31 氣泡形狀於 t=0.5 秒之局部放大圖 ... 183 圖 8.32 氣泡中心高度隨時間變化之關係圖 ... 184 圖 8.33 雷利-泰勒不穩定性之介面演進圖 (網格 48x144) ... 185 圖 8.34 雷利-泰勒不穩定性之介面演進圖 (網格 96x288) ... 186 圖 8.35 雷利-泰勒不穩定性之介面演進圖 (網格 192x576) ... 187 圖 9.1 一維飽和液體汽化問題之示意圖 ... 211 圖 9.2 一維飽和液體汽化問題之介面位置隨時間變化圖 ... 211 圖 9.3 一維飽和液體汽化問題之溫度分布圖 (t=1s) ... 212 圖 9.4 一維過熱液體汽化問題之示意圖 ... 212 圖 9.5 一維過熱液體汽化問題之介面位置隨時間變化圖 ... 213 圖 9.6 一維過熱液體汽化問題之溫度分布圖 (t=10s) ... 213 圖 9.7 水平平板薄膜沸騰之示意圖 ... 214 圖 9.8 不同網格解析度之氣泡外型比較圖 ... 214 圖 9.9 不同網格解析度之 Nusselt 數隨時間變化圖 ... 215 圖 9.10 在近臨界壓力下不同熱通量的薄膜沸騰之實驗結果[8,78] ... 215 圖 9.11 水平平板薄膜沸騰 Nusselt 數隨時間之變化圖 ... 216 圖 9.12 水平平板薄膜沸騰之動態圖( T 10C) ... 217 圖 9.13 介面和速度向量之分布圖( T 10C) ... 218 圖 9.14 水平平板薄膜沸騰之動態圖( T 17C) ... 219 圖 9.15 水平平板薄膜沸騰之動態圖( T 18C) ... 220 圖 9.16 水平平板薄膜沸騰之動態圖( T 20C) ... 221 圖 9.17 水平平板薄膜沸騰之動態圖( T 25C) ... 222

圖 9.18 水平平板薄膜沸騰之動態圖( T 30C) ... 223 圖 9.19 介面及溫度分布圖 ... 224 圖 9.20 壁面處 Nusselt 數沿 x 軸分布圖 ... 224 圖 9.21 壁面過熱度與時均化 Nusselt 數之關係圖 ... 225 圖 9.22 水平圓管薄膜沸騰之示意圖 ... 226 圖 9.23 水平圓管薄膜沸騰所使用之網格(總數 12200) ... 226 圖 9.24 水平圓管薄膜沸騰 Nusselt 數隨時間之變化圖(D0.211mm) ... 227 圖 9.25 水平圓管薄膜沸騰之動態圖(D0.211mm， T 10C)... 228 圖 9.26 水平圓管薄膜沸騰之動態圖(D0.211mm， T 20C) ... 229 圖 9.27 水平圓管薄膜沸騰之動態圖(D0.211mm， T 30C) ... 230 圖 9.28 水平圓管薄膜沸騰之動態圖(D21.6mm， T 5C) ... 231 圖 9.29 水平圓管薄膜沸騰 Nusselt 數隨時間之變化圖(D21.6mm)... 232 圖 9.30 長平板薄膜沸騰之過熱氣體初始介面分布函數 ... 233 圖 9.31 長平板薄膜沸騰 Nusselt 數隨時間變化圖 ... 233 圖 9.32 長平板薄膜沸騰之動態圖( T 30C) ... 234 圖 9.33 長平板薄膜沸騰之溫度分布( T 30C) ... 235 圖 10.1 三維網格內浸潤面向量(Swf  )之示意圖... 265 圖 10.2 三維網格內介面所跨越控制面之類型 ... 265 圖 10.3 介面通過六面體網格之類型 ... 266 圖 10.4 空間中重建後所得的三角形介面之示意圖 ... 267 圖 10.5 三維介面處總熱通量計算之示意圖 ... 267 圖 10.6 三維介面的平均面向量及平均溫度變化量計算之示意圖 ... 268 圖 10.7 三維介面處能量內邊界處理之示意圖 ... 269 圖 10.8 三維壩體潰堤介面之動態圖(300 90 90  ) ... 270 圖 10.9 三維壩體潰堤壓力及速度場之動態圖(300 90 90  ) ... 271 圖 10.10 底部水流前端與時間無因次化關係圖 ... 272 圖 10.11 水柱高度與時間無因次化關係圖 ... 272 圖 10.12 單一上升氣泡之氣泡外形 ... 273 圖 10.13 單一上升氣泡之動態圖(50 50 300  ，0.1~0.5 秒) ... 274 圖 10.14 氣泡中心高度隨間變化之關係圖(50 50 300  ) ... 274 圖 10.15 氣泡在 x=0 平面上的壓力及流線之分布圖(50 50 300  ) ... 275 圖 10.16 氣泡在中心軸上沿 z 方向的壓力變化 ... 276 圖 10.17 在 0.5 秒時氣泡中心位置上沿水平方向的壓力變化 ... 276 圖 10.18 同軸雙氣泡上升融合之動態圖 ... 277 圖 10.19 同軸雙氣泡上升融合之壓力及流線分布圖 ... 278 圖 10.20 不同軸雙氣泡上升融合之動態圖 ... 279 圖 10.21 不同軸雙氣泡上升融合之壓力及流線分布圖 ... 280

圖 10.22 雙氣泡上升融合之實驗結果[88] ... 281 圖 10.23 液滴撞擊薄層液面之動態圖 ... 282 圖 10.24 環狀水波在 y=0 截面上介面之動態圖(t=0~7ms,間隔為 1ms) ... 283 圖 10.25 環狀水波半徑隨時間之變化圖 ... 283 圖 10.26 三維水平平板薄膜沸騰平均 Nusselts 數隨時間之變化圖 ... 284 圖 10.27 三維水平平板薄膜沸騰之動態圖( T 5C) ... 285 圖 10.28 三維平板薄膜沸騰於 x=0.04m 平面上之溫度分布圖( T 5C) ... 286 圖 10.29 三維水平平板薄膜沸騰之動態圖( T 10C) ... 287 圖 10.30 三維平板薄膜沸騰於 x=0.04m 平面上之溫度分布圖( T 10C) ... 288 圖 10.31 壁面過熱度與時均化 Nusselt 數之關係圖 ... 289 圖 10.32 高溫錫液滴於油槽內介面及溫度分布之動態圖 ... 290 圖 10.33 高溫錫液滴於 x=0 平面上之溫度及流線分布 ... 291 圖 10.34 高溫正辛烷注入水槽之動態圖 ... 292 圖 10.35 高溫正辛烷注入水槽於 x=0 平面上之溫度及介面分布圖 ... 293

符號說明

符號 定義 P A ,A_C 主格點係數矩陣,相鄰格點係數矩陣 C 比熱 N C 庫倫數 d  定義在PC  方向上之向量 e d 氣泡或液滴初使直徑 D 圓管直徑 Eo Eovots 數( 2_{( ) /} e l g Eogd   ) f_ 表面張力 F 體積通量 C F 對流通量 D F 擴散通量 , g g  _{重力加速度} Gr Grashof 數( _{( )} 3_/ 2 g l g g Gr   g  ) lg h 汽化潛熱 Ja Jacob 數(JaC T_g( _wT_sat) /h_lg) * Ja 修正後之Jacob 數 (Ja*Ja/(1 0.34 ) Ja ) k 熱傳導係數 m 質傳量 Mo Morton 數( 4_/ 3 l l Mog   ) n  _{單位法向量} Nu Nusselt 數 Nu 時均化Nusselt 數

P 壓力 Pr _{Prandtl 數( /} g g g PrC  k ) q 熱傳量 Re 雷諾數 S 代數方程之源項 M S 連續方程式中之源項 f S  _{控制面之面向量} w f S  _{浸潤區域之面向量} t 時間 T 溫度 , V V  _速度 i V 在卡式座標系統之速度分量 w 權重因子 , , x y z 卡式座標系統 希臘符號  流體體積分率 ( )r  通量限制函數(r為體積分率梯度比)  平均曲率  薄膜沸騰之特徵長度 2 d  二維最危險泰勒波長 3 d  三維最危險泰勒波長  _黏滯係數  網格控制面與介面之夾角  _密度

 表面張力係數 ij  黏性應力張量 ( )   與相關的權重因子 A  控制體積內之介面面積 l  網格中心或網格面至介面之垂直距離 n  沿介面垂直方向上之距離 t  時間步階 T  壁面過熱度 V  網格體積 下標 C 主格點之相鄰格點 D 網格面之下游格點 f 網格控制面 g _氣體 int 液氣介面 l 液體 P 主格點 sat 飽和狀態 U 網格面之上游格點 UU 網格面之上上游格點 v 網格節點 w 壁面 上標 n 新時階

o 舊時階

 平均值

~ 平均平滑化之值

*,** 預測值及修正後之值

第一章 緒論

1.1 簡介

1.1.1 兩相流簡介

所謂兩相流亦即在同一系統中同時包含兩種不同狀態的物質之流動現

象並有一介面存在於兩者之間。當兩種不同性質的物質相互接觸時，交接

面處的不連續性將使其產生介面自由能(interface free energy)，此能量

係因兩物質內部的分子個別所產生的內聚力相互作用之結果。而為了維持 此介面穩定的存在，其介面自由能必須能維持正值，方可避免兩物質在介 面處產生混合之行為。因此在雙流體流場中是否發生混合，都是基於兩流 體間之介面自由能而定，而在不可混合的兩相流系統中，此介面自由能即 大家所熟知的表面張力。 根據物質的基本三態：固體、液體、氣體，兩相流可依其所含物質狀態 的種類分成：固液系統、固氣系統、液氣系統、不可混合液-液系統。其中 不可混合液-液系統，雖然並不是標準的兩相流，但在處理上由於與液氣系 統相近，因此也將其歸類為兩相流的系統之ㄧ。此外，根據流動的形態以 及介面的結構，兩相流又可分為三種基本的類型(如圖 1.1)：分層流

(separated flow)、混合流(mixed flow)、分散流(dispersed flow)[1]。

若以一個局部被水填滿的管道流為例，當管道內的水和空氣是分層流動，

流則是管內的空氣都變成小氣泡而充滿在水中(圖 1.1(c))。至於上述兩者 之間的情況則稱為混合流(圖 1.1(b))，其流動過程中包含了空氣和水的分 層流動以及小氣泡的產生。不同類型的兩相流，不僅在於介面形狀分布有 所不同，而且對於流場整體的質量、動量以及能量的傳遞也有相當大的影 響。雖然兩相流系統可分為四大類，但就屬液氣和不可混合液-液系統在工 程或是自然界中存在最為廣泛，因此本文所發展的數值方法主要以此兩系 統為主，而所分析的流動型態則為分層流，然而當所使用網格之解析度足 夠時分散流亦可視為分層流，因此本文所提出的方法也可適用於分散流之 中。 兩相流之應用常見於各式的工業系統之中，像是熱傳系統、動力系統、 生化系統、微機電系統等等。隨著科技的進步，工業產品的尺寸也逐漸縮 小，像是在燃料電池的微流道之中，流體流動時表面張力產生的效應，因 為小尺度而顯得更為重要，所以傳統上對於單一流體進行分析的方式已經 不再適宜，必須針對兩相流的系統來加以分析。 早期對於兩相流的研究主要是採實驗的方式來得到兩相流中相關的流 場流動現象和性質，然而隨著電腦科技的快速發展以及計算流體力學的邁 進，利用數值模擬來分析兩相流已經逐漸取代實驗，而一套有效而且準確 的兩相流數值方法不僅可節省實驗所需要龐大經費，且對於工業上要求的 最佳化設計以及系統安全性的預估都有相當大的幫助。然而在兩相流的數

值計算中，其困難在於它本身有著隨時間移動且多變的介面，且在介面處 流體性質存在極大的不連續性，如密度或是黏滯性，理論上此不連續性是 以階梯函數的形式來分布。在高密度比或高黏滯係數比的兩相流系統中， 此流體性質的不連續性將造成介面處的速度及壓力相當不穩定，進而導致 數值方法的發散。故在兩相流數值模擬中，絕大部分的困難來自於介面處 不連續的流體性質所造成，因此如何有效的呈現介面處的不連續性同時兼 顧計算上的穩定，成為兩相流數值方法中最棘手的問題。 在兩相流或自由表面流的數值模擬中，最困難之處在於必須考慮介面處 的跳躍邊界條件(jump condition)。假設兩流體間存在物質的不連續性且 不考慮介面處相變化所產生的質傳和熱傳，跳躍邊界可依介面上的質量與

動量守恆，分為運動邊界條件(kinematic boundary condition)及動力邊

界條件(dynamic boundary condition)，前者代表在介面的垂直方向上，

介面兩邊的流體分子與介面沒有相對運動，也就是說流體不能跨越介面； 後者則是表示作用在介面的應力必須達成平衡。在計算過程中，若介面位 置已知，則只須將上述兩組邊界導入即可。然而在每個計算的時階下，介 面位置並不是已知的變量，所以無法直接導入邊界條件，因此通常處理兩 相流的數值模擬，必須藉由運動邊界條件來求得介面的位置，再利用動力 邊界來滿足介面處的邊界條件。綜合上面所述，在兩相流模擬中，最主要 的課題有下列幾點：1.介面隨時間移動的問題 2.精確介面位置的預測 3.處

理介面上流體性質的不連續性 4.跳躍邊界以及表面張力的計算。 對於兩相流介面的計算方法，主要可根據其使用的網格類型分成兩大 類：拉格郎日(Lagranian)和尤拉(Eulerian)網格法(如圖 1.2)。拉格郎日 網格法是透過一組與介面貼齊的網格系統來計算介面隨時間移動的情況， 當介面開始變化，網格也隨之產生變形。至於尤拉網格法則是在一套固定 的網格之中，透過特殊的方式來追蹤介面運動的情形。前者常用的方法為

移動網法(moving grid mehtod)，然而此類的方法不僅在於複雜的流場內

網格隨流場嚴重變形而產生極大的誤差，並且對於其移動的邊界條件設定

上也是一大困難，因此本方法通常應用於較單純且介面無太大變形的自由

表面流之中。至於後者因其方法較為實用，所以發展出許多對於介面處理

的 特 殊 方 法 ， 如 前 端 追 蹤 法 (front-tracking method) 、 等 位 函 數 法

(level-set method)、體積追蹤法(volume tracking method)、流體體積

法(volume-of-fluid method)等。上述幾種尤拉網格法中最常被採用的是 流體體積法與等位函數法，因此後來又有許多學者擷取了兩者之優點而發 展出兩套方法結合的計算方式。而對於流體體積法而言，其根據介面計算 的方式不同又可分為介面捕捉法(interface capturing)以及介面重建法 (interface reconstruction)。前者是以類似捕捉超音速流中震波之手法 來直接計算在介面處流體體積分率跳躍的特性，而後者則採用不同的幾何 形狀來重建流體之介面，再透過體積守衡的概念來達成流體介面的推進。

1.1.2 薄膜沸騰簡介

對於一個實際的兩相流系統而言，介面處因熱傳所產生的相變化過程是 相當值得分析的問題，而沸騰(boiling)即為其中一種相變化的機制。沸騰 現象除了液體轉換為氣體的質傳行為外，更重要的是在相變化過程中，能 量以潛熱(latent heat)的形式儲存，使整個系統能夠產生高密度的熱傳。 此特點讓沸騰流動在工業應用上扮演著相當重要的角色，例如鍋爐、蒸發 器、金屬淬火、核子反應爐的冷卻、電子散熱的熱管等等。 當液體覆蓋於加熱面之上而加熱面的溫度能維持高於液體的飽和溫度 時，沸騰現象便會發生。依據流動中是否有強制對流效應，可將其分為兩

大 類 ： 池 沸 騰 (pool boiling) 和 強 制 對 流 沸 騰 (forced-convection

boiling)。當加熱面沉浸在初始為靜止狀態的液體中且沒有強制對流的效 應時，稱之為池沸騰。反之，若有強制對流的效應則為後者。 依照水的沸騰曲線(圖 1.3)，主要可以分為四大區塊，包含自然對流區、 核沸騰區、過度沸騰區及薄膜沸騰區。在自然對流區內，因壁面過熱度(加 熱面溫度與液體平均溫度之差值)較低且液體的平均溫度低於飽和溫度，雖 然可能有極少數的氣泡會在壁面上產生，但此刻整體的熱傳機制仍以液體 本身的自然對流為主。當過熱度逐漸地提高，壁面上的小穴孔處將開始產 生小氣泡，這些小氣泡稱為汽化核心。此時因相變化所產生的質傳將使小 氣泡開始成長而後脫離加熱面並在液體中上升，在氣泡成長的過程中，壁

面上的熱能將以潛熱的形式逐漸被吸收並帶離加熱面，這是沸騰熱傳能有 高效率除熱能力的主因，而此階段稱為核沸騰區。隨著過熱度再提高，氣 泡的成長速度太快以致於無法順利脫離壁面而開始在壁面處形成氣膜，此 時的氣膜不甚穩定很容易形成較大的氣泡而脫離表面，稱為過渡沸騰區。 當壁面過熱度持續提高到最低薄膜沸騰溫度(Leidenfrost point, B 點) 時，氣膜將穩定的存在並且包覆整個加熱的表面，此刻則由過度沸騰區轉 換為薄膜沸騰。 沸騰狀態屬於低溫核沸騰時，液體不斷與加熱面接觸而產生相變化，此 時熱量將透過自然對流和潛熱的形式被移除，進而增加了沸騰熱傳的強

度。然而，當熱通量高於臨界熱通量(critical heat flux, A 點)時，將由

核沸騰模式轉為過渡沸騰，此時液體因被氣膜所隔絕而鮮少與加熱表面接 觸，熱量只能依靠氣體的傳導以及輻射的方式傳遞，然而氣體本身的傳導 係數較低，因此熱通量將大大的下降。而在高溫薄膜沸騰中，壁面過熱度 的提高將使輻射效應越趨顯著，因此整體的熱通量又逐漸的提升。總而言 之，一般工業或是電子元件的操作條件都是期望在較低溫且有高熱通量環 境下，因此通常都是將整個熱通量控制在臨界熱通量以下，以維持核沸騰 的形式。 在過去的文獻中，已經有相當多關於薄膜沸騰的研究和理論被提出，然 而由於薄膜沸騰的許多現象在空間和時間尺度極小，因此大大侷限了實驗

上量測的精確度。雖然實驗上有所限制，不過也有相當多的學者已經提出 針對在某些特定條件下薄膜沸騰的半經驗式。而這些半經驗式中，較常被 後 人 所引 用 的是 Berenson[2] 所提 出 在水 平 平板 上 薄膜 沸 騰的 模 型和 Bromley[3]所推導出對於水平圓管上的薄膜沸騰之理論模型。 在水平平板的薄膜沸騰中，流體是因重力所產生的介面不穩定性而導致 氣體開始以氣泡的型態，週期性地離開平板。因此 Berenson 的理論模型是

在水平平板上一個長度為二維最危險泰勒波長(most dangerous Taylor

wavelength)的正方形區域內，存在一個半徑為R的氣泡並連接均勻厚度 的薄膜(圖 1.4)。而二維最危險泰勒波長定義如下： 2 2 3 d     (1.1) 其中為特徵長度，其值與重力(g)、兩流體密度差( _l _g)及表面張力係 數( )有關，如下式： ( l g)g       (1.2) 位在液體下方的氣體因為薄膜和氣泡兩側的壓力差而開始產生朝向氣泡之 流動，Berenson 假設此氣體的流動為層流並忽略液體側的流場影響，同時 氣泡間的距離和氣泡半徑與d₂成正比，而導出 Nusselt 數之半經驗式，如 下： 1 4 0.425( ) B GrPr Nu Ja  (1.3)

3 2 ( ) g l g g g Gr        (1.4) g g g C Pr k   (1.5) ( ) g w sat lg C T T Ja h   (1.6) 其中為密度，為黏滯係數，g為重力加速度，C為比熱，k為熱傳導係 數，T_w為壁面溫度，T_sat流體飽和溫度，h_lg為汽化潛熱。下標l和g分別代表 液體和氣體。 在 與 Berenson 相 同 的 薄 膜 模 型 下 ， Klimenko[4] 利 用 雷 諾 類 比 (Reynolds analogy)將模型從層流推展至紊流並，其半經驗式表示如下： 1 1 5 3 3 1 1 1 5 3 2 2 0.19 4.03 10 0.0216 4.03 10 K Gr Pr f for Gr Nu Gr Pr f for Gr        _{ }  (1.7) 其中 1 1/ 3 2 1/ 2 1 0.71 0.89 0.71 1 0.5 0.71 0.5 for Ja f Ja for Ja for Ja f Ja for Ja               (1.8) 上述的半經驗式表示當Ja較大的情況下，不需考慮壁面過熱度的影響。然 而，當Ja較小時則須考慮過熱度來做修正。 至於水平圓管的薄膜沸騰而言，Bromley 是最早提出理論模型的學者。

1 4 0.62( ) C GrPr Nu Ja  (1.9) 此處 Grashof 數為 3 2 ( ) g l g g gD Gr       (1.10) 其中D是圓管直徑。上式與平板薄膜沸騰的半經驗式(式 1.3)做比較，可知 Bromley 學者是透過實驗數據將前方係數由 0.425 調整為 0.62 並將特徵長 度換為圓管直徑D。 兩相流由於介面隨時間的形變及熱傳相變化機制相當複雜，因此早期實 驗上，常採用簡化後或特定狀態的模型來分析，對於數值模擬上也如同實 驗一般。一個有效而準確的兩相流數值方法，不僅要能夠精確的預測介面 的位置和形狀，更要能夠處理介面處動量和能量的跳躍邊界。早期針對兩 相流相變化的數值研究，透過假設氣液介面在變化過程中維持球形來模 擬。直到 1995 年，Welch[5]採用移動的三角網格，方能模擬二維氣液介面 可隨時間任意變形的情況。而 1997 年，Son 和 Dhir[6]則利用座標轉換的 方法，模擬一組二維軸對稱的平板薄膜沸騰現象。然而上述此兩種方法雖 然能夠維持介面鮮明的程度，然而當計算時間拉長介面開始產生大量變 形，此兩方法便無法有效的呈現薄膜沸騰的現象。隨著介面運動計算方法 的精進，許多有關薄膜沸騰的數值方法也相繼被提出。常見的方法像是

間發展了相當多有關薄膜沸騰的計算方法，包含二維水平平板及圓管的薄

膜沸騰計算乃至於三維的數值模擬。本文將再下一節中，針對兩相流介面

1.2 文獻回顧

1.2.1 兩相流之數值方法

1.2.1.1 移動網格法 (Moving Grid Method)

此方法使用一組隨流體運動的網格來計算介面的變化，亦即在所考慮的 流體中置入隨時間變形的網格，而網格邊界也就成為了所考慮流體的邊界 且各網格內的流體永遠都保持在該網格中[10,11]。本方法的優點在於其網 格邊界就是流體的介面，因此其計算後所獲得之介面可以維持相當好的鮮 明度(sharpness)和準確性。然上述的情況只發生在較簡單的流場中，當整 個的流體流動造成大量的介面變形(deforming)、破壞(breaking)或是融合 (merging)時，將使網格產生極度的扭曲，此舉將帶來計算上極大的數值誤 差且大大降低了整個運算的穩定性。在[12]中，Muzaferija 等人透過有限 體積法和移動網格配合空間守衡的概念預測自由表面流之介面形狀，雖然 在過程中流場內的網格是改以拉格郎日-尤拉法來移動，然而其計算還是侷 限在一些簡單且介面無太大變形的流場之中，並無法在一些有介面破碎或 翻轉的問題中使用。

1.2.1.2 前端追蹤法 (Front-Tracking Method)

Tryggvason 等人[13,14]提出了在尤拉網格系統中使用另一組移動的網 格來追蹤介面移動的技術，稱之為前端追蹤法。其方法主要是在介面上安 排一組相連的網格點(如圖 1.5)，利用這些網格點來追蹤介面運動的情形，

但其速度場則採用尤拉網格系統來求解動量方程式。上述的網格點構成了 固定網格系統中的一組移動邊界，此邊界依照所計算物理模型的維度可以 是二維曲線或是三維中之曲面。透過在尤拉網格上計算所得到的速度分佈 內差之後，可獲得流體前端在拉格朗日系統上所需要的傳播速度，經由此 速度可計算出流體前端網格在下一個時間的位置。本方法雖可透過前端網 格的分布獲得鮮明之介面，然而在介面傳遞的過程中其兩流體的質量並非 守恆，這並不符合流體本身的物理意義。另外，前端上的網格經過計算後 可能會被扭曲或拉伸，因此必須隨時執行重置網格(re-meshing)的動作來 調整介面上之網格，這將使整個方法便的相當繁複。再者，由於前端與速 度場是處在兩套不同的系統網格內，因此計算過程中必須隨時在尤拉和拉 格朗日系統之間作轉換，這動作也使程式化上變得相當複雜。

1.2.1.3 等位函數法 (Level-set method)

等位函數法是由 Osher 和 Sethian[15]所提出，其透過定義一組與網格 點到介面最短距離有關的連續函數()來計算介面的運動，如下式： - ( , ) 0 ( ) d x gas x t x t d x liquid     _    _  (1.11) 其中d為位置函數、表示介面所在位置。此函數根據其所在位置的流體性 質不同而具有正負號的形式，其在氣體側通常定義為負而液體側則定義為 正，因此介面將落在等位函數為零的位置上(圖 1.6)。由於流體介面是跟隨

著流體流動，因此等位函數法在介面傳輸上是直接計算等位函數的傳輸方 程式。不同於上一類的方法，等位函數只須在尤拉網格系統中直接求解此 組雙曲型式的傳輸方程式，因此在處理上方便了許多。然而，在直接求解 的過程中將會有兩個問題產生，一為數值誤差的問題，另一則是過程中會 導致質量不守恆的問題。第一個問題是來自於計算傳輸方程式時對流項之 離散所使用的數值方法造成，此問題在介面劇烈拉伸或破裂的情況甚為嚴 重。至於第二個問題則是當所定義的距離函數經由傳輸方程式計算一段時 間之後，雖然介面的等位函數還是可以維持在零階上，但是其等位函數已 無法保持距離函數的特性[16]，尤其是在介面受到扭曲變形或介面融合的 情況下。為了改善上述質量無法守恆的問題，因此[17,18]提出了在每個時 間步階下必須針對計算完成後的等位函數進行重距離化(re-distancing) 及重初始化(reinitialize)的動作來降低質量不守恆所產生的誤差。此 外，由於使用等位函數法計算時，流體的性質受到其分布所影響，在兩流 體交界處常會因其密度或是黏滯係數差異過大而導致流場計算的不穩定， 因此通常會導入 Heaviside 函數將介面拓展到 1 到 3 網格厚度，來對流體 性質作平滑化處理。

1.2.1.4 體積追蹤法 (Volume Tracking Method)

本類的方法主要是在流體中置入標記用的質點，透過計算質點的運動及

所提出的標記質點網格法(Marker-and-cell, MAC)以及 Monaghan[20]所法

展的光滑流體粒子流體力學法(Smoothed Particle Hydrodynamics, SPH)。

前者是在一套固定網格系統之中，將無質量的粒子安排在某一特定流體及

流體介面上(圖 1.7)，而這些粒子透過其鄰近固定網格上求解動量方程式所

獲得之速度來移動。至於後者則是在不需使用網格的情況下，以拉格郎日

的計算方式直接求解粒子本身的運動，與後者類似方法還有 Koshizuka[21]

等人所提出的移動粒子半隱式法(Moving Particle Semi-implicit, MPH)。

這類方法主要的特色在於可直接透過粒子分布的情形獲得介面位置，然而

當流體處在高度剪切的流場中，介面因受到拉伸而變形，此時必須使用大

量的粒子來處理介面定位之問題，這導致整個計算將相當耗費計算資源與

時間，這情況也發生在三維流場的應用上。

1.2.1.5 流體體積法 (Volume-of-Fluid Method, VOF)

流體體積法法主要是在一固定網格系統中使用體積分率來標示介面的 位置，而體積分率則是由兩流體其中之ㄧ的體積佔其網格總體積的比例來 定義。根據上述定義，可將整個計算域分為三種型態，當體積分率為 1 表 示第一種流體而 0 則表示另一種，至於體積分率介在 0 至 1 之間的網格則 表示介面存在的網格(圖 1.8 所示)。經由計算流體體積分率之傳輸方程式 (即 VOF 方程式)可以獲得介面隨時間傳遞的情形，雖然此傳輸方程式與等 位函數所使用的形式相同，然而處理上卻是複雜很多，這起因於體積分率

分布的情況在介面處會呈現一個階梯函數的形式，相對於等位函數則是呈 現連續分布的情況。自從七零年代開始已經有許多關於流體體積法的研究 被提出來，而大致上可以分為兩大類：介面捕捉法與介面重建法。前者是 以數值離散的觀點來直接求解流體體積方程式，後者則是採幾何的方式來 重建介面並配合體積守恆的概念來計算流體介面之推進。以下將針對兩種 方式做深入的介紹。

1.2.1.5.1 介面捕捉法 (Interface Capturing Method)

所謂的介面捕捉法亦即利用類似於在高速流動中捕捉可壓縮流震波之

數值差分方式來離散 VOF 方程式中的對流項，以此方式直接來獲取流體介

面 處 體 積 分 率 跳 躍 的 特 性 。 此 類 的 方 法 所 遭 遇 到 的 問 題 有 數 值 擴 散

(numerical diffusion)以及數值震盪(numerical dispersion)。前者造成

介面跳躍的特性被抹滅(smearing)，後者則常使介面產生產生不穩定性且

使體積分率產生不合其物理特性的無界限性(unboundness)。這類方法較低

階精度的有一階上風法(upwind differencing scheme)和下風法(downwind

differencing scheme)，上風法是擁有強烈數值擴散性的差分法，其可使

計算過程相當穩定而不發散，然而此特性將使介面無法達到所需要的鮮明

度。而下風法本身具有離散時所給予的反擴散性，因此對於介面有高度的

壓縮性，此性質可使介面受到強烈的壓縮而成一階梯(step)函數的形式，

不穩定並使體積分率無法維持界限性。單獨使用上述的兩種差分方式無法 獲得適用於兩相流或雙流體流場之介面，因此 Hirt 和 Nichols [22]在其研 究中提出以主體受體法(donor-acceptor scheme)的方式來混合上風法及下 風法，來結合上述兩種差分法之優點。另外，Rudman[23,24]提出以通量修 正傳輸法(Flux-Corrected Transport, FCT)[25]來計算介面的傳輸，此方 法分為兩個步驟，首先利用上風法來求取體積分率之傳輸，而修正步驟則 採下風法的方式來減少上風法所造成數值擴散的問題，然而此方法以一維 的定理所推出，較不適用於非結構性網格之中。 上述皆使用較為低階的離散方式來處理 VOF 方程式，因此同樣也可採高

階的差分方式來計算 VOF 方程式，例如：正規化變數(Normalizes Variables

Diagram, NVD)或全變量縮減(Total Variation Diminishing, TVD)之差

分法。然而直接採用上述這兩類的差分法並無法完全的避免數值擴散發 生，因此衍生了另一套混合(blending)的方法來減少數值擴散以維持介面 的 鮮 明 度 。 此 類 型 的 混 合 法 其 主 要 的 手 法 是 選 用 一 組 高 解 析 度 (high-resolution)差分法和另一組壓縮性(compressive)差分法來加以混 合，而其混合的方式為使用一個與網格控制面和流動方向之夾角有關的轉 換函數來做兩者之間的切換。近年來較著名的有 Ubbink 和 Issa[26]所提出

的 壓 縮 性 介 面 捕 捉 法 (compressive interface capturing scheme for

法(high- resolution interface capturing scheme, HRIC)及 Darwish 和

Moukalled [28] 對 流 切 換 介 面 捕 捉 法 的 (switching technique for

advection and capturing surface, STACS)等。上述的混合法都是符合所

謂的對流界限準則(Convection Boundness Criterion, CBC)[29]，屬正規

化變數差分法之類型，雖然透過混合的方式可減少數值擴散，然而其所計

算出的介面分布通常仍分散於數格網格的厚度之間。

1.2.1.5.2 介面重建法 (Interface Reconstruction Method)

介面重建法主要是經由定義流體介面幾何形狀與介面的法向量，在透過 局部體積守恆的概念加以重建介面。在介面重建完成之後，根據局部的速 度分布與介面形狀來計算介面推進時所需的體積通量，最後可依體積通量 來估算介面移動時所需的傳輸量，進而可獲得新的體積分率。因此本類型 之方法主要分為三個步驟：介面幾何重建、網格面上體積通量計算、計算 介面傳輸方程式。大部份介面重建法可依所使用的介面幾何形狀和介面法

向量定義的方式，主要可分為下列三種：SLIC 法 [30]、Hirt-Nichols＇ VOF

法 [22]及 PLIC 法[31]。

SLIC(simple line interface calculation)法[30]是由 Noh 和 Woodward

所提出。在每個網格中，介面是假設為片段常數(piece-constant)分布之

情形，亦即介面幾何是被假設為與網格線平行的直線。本方法是屬於方向

向有關，掃視的方向可分為 x 方向掃視(x-sweep)和 y 方向掃視(y-sweep)。 對於 x 方向掃視而言，介面法向量是經由介面所在網格的右方和左方的網 格之體積分率來決定；同理，y 方向則由上下方鄰近格點來決定。在雙流體 流場中，使用 SLIC 法沿單一方向掃視下，其介面分布的情形可分為下列四 種情況(圖 1.9)：(1)兩側為不同流體。(2)兩側為同一流體。(3)兩側皆含 介面。(4)單側含介面。上述為 SLIC 法中介面重建的方式，當此步驟完成 後，接下來便是求取介面移動時所需網格面上的體積通量同樣採用方向分 離的方式來進行運算，亦即在當下所掃視的方向來進行面上體積通量的計 算。最後必須計算傳輸方程式以得到下一時階的體積分率，在 SLIC 法中通 常使用運算子分離(operator split)的方式來處理。

Hirt 和 Nichols＇VOF 法[22]中介面重建的方式與 SLIC 同樣是採用水

平線或垂直線並將體積分率視為單值函數(single-value function)，然後

透過計算主格點周圍九個網格的斜率，由斜率大小來判斷介面的方向並對

主格點內的介面做定位。至於面上體積流率的方式則如同前小節所介紹是

採用主體受體法來計算，根據重建後介面的方向與流動方向相對關係來決

定使用上風法或下風和上風之混合。

PLIC(piecewise linear interface calculation)法[31]是由 Youngs

所提出的一種高階精度的介面重建法。其所假設的介面幾何為兩端點在網

向量可推導出介面在網格內所截斷的體積，然後透過截斷體積與網格內局 部體積分率來進行介面重建。上述的重建過程中，介面所截斷體積必須與 網格內所考慮的流體體積相等，因此有直接或疊代[32,33]求解的方式來重 建介面。直接求解以 case-by-case 的方式來重建介面，在一個二維系統中 根據介面法向量與 x 座標軸的夾角和網格內所考慮流體之體積可將介面在 網格內的分布情形分為 16 種類型(三維可達 64 種情形)，經由翻轉與對稱 可分為 4 種基本類型(圖 1.10，夾角為0 ~ / 2 )，因此當夾角決定後可知介 面所屬類型，再進一步由網格所含的流體體積來重建介面。直接求解必須 逐一針對包含介面的網格進行判斷，這方式相當的繁瑣，因此有所謂疊代 的求解方式來取代之。疊代法主要是透過介面法向量來求取介面的方程 式，再由介面方程式推出介面在網格內所截斷的體積，最後獲得截斷體積 與網格內流體體積的關係式，經由疊代法計算此關係式以重建介面。當上 述介面重建完成之後，同樣必須計算網格面上所需的體積通量，其方式可 分為尤拉(Eulerian)式和拉格郎日式(Lagrangian)。所謂尤拉式就是根據 尤拉網格系統中的相關已知速度場的資訊配合所重建出的介面幾何形狀來 進行計算，而拉格朗日式則是透過直接計算介面移動後的位置來獲取所需 的面上體積通量。當計算出介面傳遞時所需的體積通量之後，便可求解介

面傳輸方程式。Ashgriz 和 Poo[34]改進 PLIC 法而發展成為 FLAIR(flux

其方法主要是根據網格面相鄰之兩個網格的體積分率來獲得可穿越網格面

的斜直線以進行介面的重建，此介面重建完成後，網格面上所需的體積通

量則可直接沿此介面進行積分而得。圖 1.11(a)表示實際流體介面分布的情

形，而(b)至(e)則分別顯示使用 x 和 y 方向掃視的 SLIC 法、PLIC 法以及

FLAIR 法所重建後的介面形狀示意圖。 綜觀上述所介紹的介面重建法，其基礎都是來自於介面法向量的定義， 此法向量的準確與否將直接影響整個介面重建過程的精準度，間接也影響 到整個速度場及表面張力的計算。然而體積分率在計算域中是呈現一個極 度不連續的分布，此現象將造成介面法向量在計算上很難達到很高的精確 度，因此有許多研究都著重於上述這個困難的突破上。雖然精度較高的 PLIC 法以被使用在相當多的應用上，然而其方法基本上仍有兩個問題，第一是 若是在不規則的網格中使用，其介面推進中所需的網格面上體積通量在計 算上將會變得相當複雜且困難。如前所述，第二個困難將會發生在三維的 流場之中，介面在網格內分佈的情況將會有 64 種，雖然根據對稱及旋轉的 簡化可化約 5 種[35]，然而在程式化上的步驟仍然相當繁瑣。針對第一個 問題 Mosso[36]等人透過 Lagrangian-Eulerian 的方式來進行處理，並有許 多相關的研究將其應用到三角型的網格系統之中[37,38]。

1.2.1.6 流體體積法與等位函數法之混合

上述已經談過等位函數法雖然可以得到準確的介面位置並且其函數是呈

現連續性的分布可以有效地計算出介面的曲率以及方向，然而在計算過程 中會有質量損失的問題。至於在流體體積法中若使用介面追蹤的方式，將 會因為數值擴散的問題導致介面位置不明確，然採用介面追蹤的方式來計 算雖可得到確切的介面幾何和形狀，但是由於體積分率在介面處為梯度函 數的分布，這將造成在計算介面法向量所需的體積分率梯度及介面曲率上 之困難，因此有許多學者提出了結合等位函數及流體體積的方式來進行介 面傳輸的計算[39-43]，其方式就是透過求解 VOF 方程式來計算介面推進的 情形而介面曲率及法向量則由等位函數法所計算出。

1.2.2 薄膜沸騰之數值模擬

1.2.2.1 前端追蹤法

1998 年，Juric 和 Tryggvason[7]採取前端追蹤法[13]並在連續及能量 方程式中加入介面處跳躍邊界所導出之源項的概念，提出了一些有關二維 水平平板薄膜沸騰的計算結果。Shin 和 Juric[44]採前端追蹤法配合等位 面法(level-contour)重建介面的方式，使氣液介面可在三維的流動中有效 的合併或破裂，進而模擬三維池沸騰現象。而 2004 年，Esmaeeli 和 Tryggvason[45,46]以相同的方法研究有關飽和液態水與過熱蒸汽所形成 的池膜沸騰現象，其不僅將二維單顆氣泡之薄膜沸騰推展至三維，並著手 研究不同的流體性質對多模式(multi-mode)薄膜沸騰的影響。上述均為平 板膜沸騰的問題，而為了能處理幾何外觀較為複雜的圓管薄膜沸騰，

Esmaeeli 和 Tryggvason[47]導入嵌入邊界法(immersed boundary method) 來處理水平圓管上的沸騰現象。雖然有關前端追蹤法於沸騰熱傳的計算已 經有相當程度的發展，但是如同前小節所介紹，此法在兩種座標系統中的 轉換是相當困難而且複雜，尤其是在三維的問題之中。

1.2.2.2 等位函數法

有別於使用前端追蹤法，Son 和 Dhir[8]修改了等位函數法並考慮軸對 稱的假設下，來模擬近臨界壓力下薄膜沸騰在不同的壁面過熱度所產生的 影響。其發現在低過熱度的情況下，氣體以氣泡的方式離開加熱面。然而， 在高過熱度的情況下，則會開始以噴流(jet)的方式離開。Son[48]引入介 面處相變化的效應並改善整體質量守恆的概念修改此法後，模擬在靜止水 中上升氣泡及成長的問題。本方法因流體性質不連續性將造成的數值不穩 定 問 題 ， 因 此 常 在 氣 液 介 面 處 加 入 平 滑 脈 衝 函 數 (smoothed delta function)[17]以增加介面處的擴散性。然而此舉將使能量的跳躍邊界條 件，無法準確的導入整個系統之中。另外一套在等位函數中加入虛擬網格

(ghost fluid method)[49,50]的方法常被用來改善上述的問題，此方法可

有效的維持介面處性質的不連續性。因此許多學者[51-53]相繼採用這個方

式來模擬含有相變化的流場。而 2007 年，Son 和 Dhir[54]也使用相似概念

模擬了二維和三維水平圓管上的薄膜沸騰。更在[55]的文章中，詳盡地討

1.2.2.3 流體體積法

Welch 和 Wilson[9]是最早將流體體積法實行於相變化流場的學者，其 採用 PLIC 法來重建介面並在連續方程式中導入相變化所產生的質傳量來模 擬水平板上薄膜沸騰的行為。在 2002 年，此方法更進一步用於分析包含固 體平板和薄膜沸騰之間的熱傳現象[56]。Agarwal[57]等人則是利用此法模 擬近臨界壓力下的平板薄膜沸騰，在模擬過程中，其考慮液體為飽和態且 過熱氣體的熱力性質隨溫度而改變，結果顯示過熱態性質是否為常數，對 於整體熱通量的影響甚鉅。2011 年，Guo[58]等人更進一步結合上述的方 式，再配合等位函數法求取更精確的介面定位，來模擬不同性質下不同的 薄膜沸騰現象。

1.3 研究目的

本研究主要在發展可應用在非結構性網格中之兩相流的數值模擬系

統，其主要處理內容包含了介面傳輸的計算、介面運動和動量與能量方程

式之結合及相變化所產生的熱傳和質傳行為。

本文將提出兩套基於流體體積法的介面計算方式，包含通量混合介面捕

捉法(flux blending interface capturing scheme, FBICS)[59]及守恆內

差 介 面 捕 捉 法 (conservative interpolation scheme for interface

tracking, CISIT)[60]。前者是類似於捕捉可壓縮流中的震波之方式，直 接計算 VOF 方程式中所需要的體積分率的對流通量來計算介面移動的情 形，而其對流通量主要採用通量混合的方式，透過一連續性的函數在高解 析差分法及壓縮性差分法之間做切換，以同時達到抑制數值擴散和流場穩 定之效果。後者則屬於介面追蹤法的範疇，其主要是採用一種簡單且有效 的介面重建方式，並透浸潤面積的定義來獲取通過網格面之流體體積通 量，再以先預測再修正的方式來求解介面傳輸方程式。 本文在流體體積方程式與速度場結合的處理上，主要是考慮整個流場為 均質且連續之流體，亦即流體在介面處的速度及壓力是連續的。經由體積 分率來控制流場中流體性質的分布情形並在動量方程式中加入由體積分率 梯度所表示的表面張力項來模擬雙流體或兩相流流場之流動。至於動量方 程式與連續方程式的耦合，則是採用 PISO 演算法的方式來進行。

對於含相變化的兩相流之計算則是將相變化模型導入 CISIT 法中來處 理。由能量和質量的跳躍邊界條件，推導後可得介面處因相變化所產生的 質傳與熱傳關係式，將此質傳量導入連續方程式的源項中即可模擬含相變 化的兩相流，而能量方程式則是將 CISIT 法所重建的介面視為能量一組內 邊界來處理，再以全隱式法求解離散後的能量方程式。另外，藉由介面處 溫度及熱傳量連續性的條件來取代跳躍邊界條件，亦可將本方法推展至無 相變化的雙流體熱傳問題之中。下列章節將針對本文所使用的統御方程式 及詳細的數值方法做介紹，而結果討論中主要可以分為四大部分：

1. FBICS 法於自由表面流(free surface flows)之應用。

2. CISIT 法於雙流體流場(two-fluid flows)之應用。

3. CISIT 法結合相變化模型於二維薄膜沸騰之模擬分析。 4. CISIT 法於三維兩相流問題之應用。 前兩個部分是藉由一些有理論解的模型來驗證兩套方法的精確度，並模擬 流場結構較為複雜之實際模型與實驗上的結果作比較，以說明此兩種數值 方法於真實流場中的可行性。第三個部份則是將 CISIT 法結合本文所提出 的相變化模型及能量場的求解方式，來模擬二維水平平板及圓管之薄膜沸 騰現象，驗證本文所提出之相變化模型在實際沸騰問題中的效能。最後的 部份則是將整個基於 CISIT 法的模型推展至三維的數值模擬之上，以說明 本文提出的數值方法在實際問題的應用分析上之廣泛性。

(a) 分層流

(b) 混合流

(c) 分散流

(a) Lagranian

(b) Eulerian

圖 1.3 水沸騰曲線

圖 1.5 前端追蹤法示意圖

圖 1.7 MAC 法示意圖

(1)兩側為不同流體 (2)兩側為同一流體

(3)兩側皆含介面 (4)單側含介面

圖 1.10 PLIC 法簡化後介面於網格內分布情形(介面方向與 x 軸之夾角為

(a) 真實介面

(b) SLIC 法(x-sweep) (c) SLIC 法(y-sweep)

(d) PLIC 法 (e) FLAIR 法

第二章 數學模型

2.1 簡介

本章節將介紹本文中所使用含相變化之兩相流模型的統御方程式，包含 連續方程式、動量方程式、能量方程式及流體體積(Volume-of-Fluid, VOF) 方程式。本研究採流體體積分率(volume fraction)來定義兩種不可互溶之 流體，並且利用一步階函數(Step function)的流體體積分率來表示兩流體 之介面。在動量方程式計算上，將兩流體視為均質且連續來處理，並且透 過流體體積分率來使各流體有著不同的流體性質；在介面上則是採用

CSF(Continuum Surface Force)表面張力模型[61]，此模型是利用一個表

面力的方式來近似兩流體在介面處的跳躍邊界條件(jump condition)。 至於相變化的模型則是假設溫度分布為連續，亦即兩種流體在介面處的 溫度相同且與系統壓力對應之飽和溫度相等，如此可由介面兩側流體各別 的於介面的法線方向上之溫度變化量來計算介面處的熱通量，再透過潛熱 便可獲得介面處因相變化所產生的質傳量，最後將其導入至連續方程式中。 傳統上，對於能量方程式處理的方式有兩種，一者是將上述熱通量視為 能量方程式的源項，另一則是透過數值近似直接處理介面處的熱傳。由於 介面處常呈現流體及熱力性質極度不連續的情形，此現象造成前者在數值 計算上相當不穩定，因此本文採用後者的方式，將介面視為一內邊界的概 念來處理。