超臨界二氧化碳在板式熱交換器中的熱流場之數值分析

國立交通大學

機械工程學系

碩士論文

超臨界二氧化碳在板式熱交換器中的熱流場之

數值分析

Numerical Analysis on the Thermal Flow Field of Supercritical

Carbon Dioxide in Plate Heat Exchanger

研 究 生：湯宜群

指導教授：王啟川 博士

國立交通大學

機械工程學系

碩士論文

超臨界二氧化碳在板式熱交換器中的熱流場之

數值分析

Numerical Analysis on the Thermal Flow Field of Supercritical

Carbon Dioxide in Plate Heat Exchanger

研 究 生：湯宜群

指導教授：王啟川 博士

超臨界二氧化碳在板式熱交換器中的熱流場之數值分析

Numerical Analysis on the Thermal Flow Field of Supercritical

Carbon Dioxide in Plate Heat Exchanger

研 究 生：湯宜群 Student：Yi-Chun Tang 指導教授：王啟川 Advisor：Chi-Chuan Wang 國 立 交 通 大 學 機 械 工 程 學 系 碩 士 論 文 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 Master of Science

In

Mechanical Engineering July 2013

Hsinchu, Taiwan, Republic of China

超臨界二氧化碳在板式熱交換器的熱流場之數值分析 研究生:湯宜群 指導教授:王啟川 國立交通大學機械工程學系碩士班 摘要 本研究開發以二氧化碳為工作流體的板式熱交換器模擬程式，模 擬內容包括二氧化碳在板式熱交換器中的流速分佈、壓降及熱傳等。 模擬條件包括冷媒及水側的質量流率、入口壓力、板片數量以及不同 的板片等。模擬結果顯示在水對二氧化碳板式熱交換器中，由於水的 物性隨溫度變化並不明顯，因此影響其流場分佈結果主要為溫度變化 所導致。而在二氧化碳流場中，由於超臨界二氧化碳的物性變化相當 劇烈，因此密度及黏滯係數所造成的影響將不可忽略。除此之外，二 氧化碳在靠近臨界溫度時，其熱傳效果將有顯著的提升，未來在設計 熱交換器時，必須將此一特性考慮進去。透過此一程式，能夠對二氧 化碳應用於熱交換器上能夠有更進一步的瞭解，為熱交換的設計及改 進性能有相當大的幫助。 關鍵字：二氧化碳、板式熱交換器、流速分佈、熱傳

Numerical Analysis on the Thermal Flow Field of Supercritical Carbon Dioxide in Plate Heat Exchanger

Student: Yi-Chun Tang Advisor: Chi-Chuan Wang Abstract

In the present study, a numerical simulation program is developed to calculate the pressure drop, velocity distribution and temperature of both water and supercritical carbon dioxide in a plate heat exchanger. The critical state of carbon dioxide is at 7.8 MPa and 330 K, respectively. Therefore the heat exchanger is operated above critical point as a gas cooler. Since the physical properties of carbon dioxide vary drastically near the critical point, therefore the conventional numerical simulation methods for plate heat exchanger are not available. The objective of this study is to develop a carbon dioxide based program which includes the pressure drop, velocity distribution, temperature and physical properties such as density, viscosity and heat capacity. By using appropriate

correlations of the Nusselt number and the friction coefficient, the flow and the temperature distribution can be correctly calculated. The effect of the pressure drop, velocity distribution and temperature was calculated with different mass flow rate, inlet pressure, inlet velocity, inlet

temperature, number of channels and different plates. In a water to carbon-dioxide plate heat changer, the simulation results show that the flow distribution of water did not change significantly with temperature. Due to the physical properties of carbon dioxide vary drastically at supercritical state, the density and viscosity of the impact cannot be ignored. For an example, when temperature variation between the plates change rapidly, it will cause considerable flow distribution effect ,and may change the flow velocity distribution pattern change from decreasing type into increasing type. In addition, the heat transfer effect will be significantly improved when the carbon dioxide near the critical temperature, when designing of heat transfer, it could be taken into account of the characteristics.

Keywords: carbon dioxide; plate heat exchanger; flow distribution; heat transfer

誌謝 為期兩年的研究所生涯眨眼就過去了，在這兩年中，我學習了許 多事情，也擁有了許多珍貴的回憶，同時，也要感謝這一路上支持幫 助以及為我加油的師長、朋友以及家人們。 首先要感謝我的恩師王啟川 博士的辛勤指導，老師學識淵博、 待人謙和、處事圓融有條理，培養我養成獨立思考的習慣以及讓我認 識什麼才叫做研究，使我獲益良多，在此向老師致上最深摯的謝意。 感謝口試委員楊建裕 博士以及謝文德 博士在研究以及論文的 寫作上的建議以及指導，使本文能夠更加完整，也使學生受益匪淺。 感謝實驗室的博班學長懷保及昆壕在課業上以及課業下的幫助， 在我有困難時，學長們總是不遺餘力的幫助我。感謝實驗室的同學們， 建宏、宥澄、慧敏，在研究所生涯的支持與幫助，沒有你們，也許這 篇誌謝不會這麼早出來。還有感謝來自各地的同學們，白雲哥建宇、 仁甫哥茂銓、帥哥資翰、蘿蔔敦仁、小翼與榮致以及最後還有在我忙 碌時陪我聊天、在我結案時揪我出遊、在我吃飽時找我吃飯的白爛哥 建任。還有實驗室的學弟們，丁博、啓善、博元、宗瀚、能傑、裕傑、 屁王創褘及顧家冠丞，為實驗室帶來了許多笑聲，使我不空虛寂寞覺 得冷。 最後感謝遠在故鄉，台南新市的所以好朋友們，GOGOGOGO!!!! 阿，還有我家人，我永遠愛你們!! 謹以本篇論文獻給我的家人，師長以及所有的朋友們。

目錄

摘要………i Abstract………ii 誌謝………iii 目錄………iv 圖目錄……… vi 第一章 緒論………1 1-1 研究背景………1 1-1.1 二氧化碳………1 1-1.2 板式熱交換器………2 1-2 研究動機與目的………8 1-3 文獻回顧………10 1-3.1 超臨界二氧化碳熱傳特性………10 1-3.2 流場計算之文獻………17 1-3.3 熱傳計算之文獻………22 第二章 物理模式………25 2-1 板式熱交換器幾何尺寸………25 2-2 分析假設與統御方程式………27 第三章 數值計算模式與分析方法………29 3-1 模型介紹………30

3-2 方程組的建立………32 3-3 解題流程………39 3-4 牛頓法的介紹………40 3-4.1 牛頓-拉弗森法………41 3-4.2 牛頓下降法………43 3-5 程式流程………44 第四章 研究成果與討論………49 4-1 格點測試………50 4-2 模擬驗證………53 4-2.1 流場驗證………53 4-2.2 溫度場驗證………57 4-3 流場結果………60 4-4 溫度場結果……….64 4-4.1 溫度場結果分析………64 4-4.2 改變冷媒流量………70 4-4.3 改變二氧化碳壓力………75 第五章 結論………79 5-1 結論………79 5-2 未來改進方向 ………80 第六章 參考文獻………81

圖

目

錄

圖 1.1 板式熱交換器實體圖……… 4 圖 1.2 板式熱交換器內部構造示意圖……… 4 圖 1.3 板式熱交換器內部構造示意圖……… 5 圖 1.4 板式熱交換器流動示意圖……… 6 圖 1.5 (a) U 型板式熱交換器配置圖……… 7 圖 1.5 (b) 型板式熱交換器配置圖……… 7 圖 1.6 典型超臨界流體 CO2 的系統循環圖………9 圖 1.7 超臨界下二氧化碳密度變化……… 9 圖 1.8 Gnielinski Correlation 和實驗數據比較………14 圖 1.9 超臨界下二氧化碳比熱變化………15 圖 1.10 超臨界下二氧化碳導熱性變化………15 圖 1.11 超臨界二氧化碳黏滯性變化………16 圖 1.12 超臨界下二氧化碳密度變化………16 圖 1.13 流體進入流道示意圖………21 圖 1.14 入口控制體積示意圖………21 圖 1.15 出口控制體積示意圖………21 圖 1.16 熱傳控制體積示意圖………22 圖 2.1 常見板式熱交換器幾何參數其及定義………25

圖 3.1 冷側數值模型配置圖………31 圖 3.2 熱側數值模型配置圖………31 圖 3.3 分流區控制體積示意圖………34 圖 3.4 匯流區控制體積示意圖………34 圖 3.5 流道控制體積示意圖………34 圖 3.6 冷側熱場控制體積示意圖………37 圖 3.7 熱側熱場控制體積示意圖………37 圖 3.8 分流區不同流速下主流方向所對應的損耗係數………44 圖 3.9 匯流區不同流速下主流方向所對應的損耗係數………45 圖 3.10 流場程式流程圖………46 圖 3.11 熱場程式流程圖………47 圖 3.12 整合後程式流程圖………48 圖 4.1 水對水熱交換器格點測試結果………52 圖 4.2 水對 CO2 熱交換器格點測試結果………52 圖 4.3 文獻所使用板片幾何圖形……… ………54 圖 4.4 摩擦因子與雷諾數關係圖………54 圖 4.5 壓降與流道編號關係圖………55 圖 4.6 無因次化後速度與流道編號關係圖………56 圖 4.7 溫度場結果與文獻比較………59

圖 4.8 等溫時不同入口雷諾數對流速分佈的影響(n=10)…………61 圖 4.9 等溫時不同入口雷諾數對流速分佈的影響(n=40)………61 圖 4.10 質量流率相同、溫度不同流場分佈情形(n=10)………63 圖 4.11 質量流率相同、溫度不同流場分佈情形(n=40)………63 圖 4.12 水側出口溫度分佈………67 圖 4.13 二氧化碳側出口溫度分佈………67 圖 4.14 水側無因次化流速分佈………68 圖 4.15 二氧化碳側無因次化流速分佈………68 圖 4.16 水側熱傳量與流道關係圖……….69 圖 4.17 二氧化碳側熱傳量與流道關係圖………69 圖 4.18 水側速度分佈圖………72 圖 4.19 二氧化碳側速度分佈圖………72 圖 4.20 C=0.1 時，分流區及匯流區壓力變化情形………73 圖 4.21 C=1 時，分流區及匯流區壓力變化情形………73 圖 4.22 不同 C 值時，二氧化碳流道熱傳量示意圖………74 圖 4.23 C=1 時在不同壓力下各個流道的熱傳情形………77 圖 4.24 C=0.25 時在不同壓力下各個流道的熱傳情形………77 圖 4.25 C=0.25 時，壓力為 10 MPa 及 12 MPa 之比熱變化情形………78 圖 4.26 C=1 時，壓力為 10 MPa 及 12 MPa 之比熱變化情形………78

符號說明

A ：截面積，m2

Aw ：流道有效熱傳截面積，m2

C ：熱容比(heat capacity ratio)，C =(mCp)min / (mCp)max

CP ：比熱，J/kg⋅k Ct ：流道入口水頭損失候係數 Cto ：流道入口水頭損失候係數 D ：入/出口處導管直徑，m Dh ：水力直徑, m De ：等效直徑，m F ：Dracy 摩擦因子 h ：熱傳係數，W/m2_⋅ K KL ：分流區水頭損失係數 KC ：匯流區水頭損失係數 l ：格點長度 (l =L N/ )，m/個 L ：板片長度，m Lt ：[9]中入口處板片間距(即入口導管長度)，m ṁ ：質量流率，kg/s n ：流道數

N ：流道格點數，個 P ：流道截面濕周長，m 𝑃_𝑖𝑖 ：入口壓力，Pa 𝑃_𝑜𝑜𝑜 ：出口壓力，Pa 𝑃_𝑐ℎ ：流道壓力，Pa ∆P_ch ：流道壓力差 t ：板片厚度，m 𝑢_𝑐 ：[9]中流道方向的速度，m/s U ：總熱傳係數，W/m2_⋅ K 𝑉_𝑖𝑖 ：入口處(分流區)速度，m/s 𝑉_𝑜𝑜𝑜 ：出口處(匯流區)速度，m/s 𝑉_𝑐ℎ ：流道速度，m/s 𝑉_𝑐ℎ𝑚 ：流道平均速度，m/s w ：板片寬度，m W ：[9]中主流方向的速度，m/s Z ：流道入/出口處軸方向 z ：無因次化後流道位置 ，z= x L/ ，無因次 特殊符號 β ：由主流流入流道之平均速度比， Wc W β =

λ ：山型紋間距，m

θ

：無因次化後的溫度 C H C T T T T θ = − − ，無因次 w

τ

_{：管壁摩擦剪力，Pa} c

υ

：流道無因次化速度，

υ

c=Vch / Vchm，無因次 ζ_c ：流道中流體平均總水頭損失係數 上標 * ：出口狀態 下標 c ：流道 C ：二氧化碳側(CO2) in/out ：入/出口 H ：水側(H2O) i ：表示第 i 個流道處 j ：表示流道中第 j 個格點

第一章

緒論

1-1 研究背景 1-1.1 二氧化碳 在早期的冷凍空調系統中，所使用的冷媒皆以天然冷媒為主，包 括氨、丙烷及二氧化碳等。由於氨及丙烷因具有可燃性及毒性，當考 慮到使用場合以及安全性等因素時，則受到相當大的限制。至於二氧 化碳則受限於操作壓力過大以及早期的工業水準等因素，應用上十分 有限。1930 年代起，隨著 CFC（合成冷媒）的應用與發展，人們認 為此類冷媒具有無毒、無味、安定性高、不可燃等優點，故被廣泛應 用在冷凍空調系統之中，使得天然冷媒逐漸被淘汰。而直至 1980 年 代起，臭氧層的破壞與溫式效應等議題逐一浮現，使得人們對於使用 CFC、HCFC 甚至是 HFC 這類合成冷媒產生質疑，而目前較常使用 之冷媒如 R-134a，雖無破壞臭氧層等問題，但其仍具有極高之溫室 效應潛力。而其它研究中之各種鹵烷冷媒等，亦具有相同之問題。為 了解決合成冷媒對於環境的影響與衝擊，天然冷媒重新被人們所重視 及發展。其中，二氧化碳因為其取得容易，無毒、不可燃、不助燃且 無破壞臭氧層及溫室效應潛力低等優點，為一理想無汙染冷媒。

1-1.2 板式熱交換器 常見的製程用熱交換器依其結構約可分為套管式熱交換器、殼 管式熱交換器及板式熱交換器(外觀及配置如圖 1.1-1.3)等三類。由於 熱交換器高效率及小型化等需求，板式熱交換器被普遍的應用。而本 研究所探討的對象為板式熱交換器。 板式熱交換器主要由分、匯流管、板片以及墊片組成，藉由墊片 和板片不同的組合方式，冷熱流體將經由分流管分流至相鄰的板片通 道流動，而達到熱交換的效果(如圖 1.4)，再於匯流管結合將冷熱流 體導出。板式熱交換器若以流動方向區分，則可以分為順流和逆流兩 種；而若以分、匯流管出入口位置的不同則可分成 U 型和 Z 型兩種(圖 1.5)。 板式熱交換器擁有以下優點[1]： 1.易清潔、檢查及保養。 2.可隨負載而增減熱傳面積→藉由板片數、板片大小、板片型式、 流場安排等因素之變化(針對組合式而言，硬焊式無此優點)。 3.低污垢阻抗→因內部流場通常是在高度紊流情況下，故其污垢 阻抗只有殼管式之 10~25%。 4.熱傳面積大→具高熱傳係數、低污垢阻抗、純逆向流動，故在 同熱傳量下，熱傳面積約為殼管式之 1/2~1/3。

5.低成本。 6.體積小→同熱傳量下，體積約為殼管式之 1/4~1/5。 7.重量輕→在相同熱傳量下，重量約為殼管式之 1/2。 8.流體滯留時間短且混合佳→可達到均勻之熱交換。 9.容積小→含液量少、快速反應、製程易控制。 10.熱力性能高→溫度回復率可達 1°C，有效度可達 93%。 11.無殼管式中流體所引起之振動、噪音、熱應力及入口沖擊等問 題。 12.適合液對液之熱交換、需要均勻加熱、快速加熱或冷卻之場 合。

圖 1.1 板式熱交換器實體圖[11]

圖 1.3 板式熱交換器內部構造示意圖

圖 1.5 (a) U 型板式熱交換器配置圖[9]

1-2 研究動機與目的 由於二氧化碳的臨界溫度與臨界壓力分別為 31.1°C與 78.7 Bar， 在一般的操作環境下，二氧化碳將持續維持在臨界點以上(圖 1.6)， 而二氧化碳在此區域的液態、氣態之分界並不明顯，且其在熱交換器 內部的熱傳現象並不同一般冷媒經由冷凝器以冷凝的方式散熱，而是 以近乎氣體的方式散熱，此與傳統冷媒有相當大的差別。而由於在散 熱過程中，二氧化碳並無實際冷凝現象發生，一般而言，又稱氣體冷 卻器(Gas cooler)而不稱為冷凝器。另外，由於二氧化碳位於超臨界狀 態的物性隨溫度及壓力的變化相當劇烈(圖 1.7)，因此在計算上，傳 統氣對氣熱交換器理論方程式，在此並不適用，綜合上述，本研究目 的在於建立一個描述二氧化碳在臨界點附近之程式，藉由此程式模擬 二氧化碳在Ｕ型板式熱交換器中的流速、壓力與溫度分佈，並從中找 出增加熱傳性能的方法。

圖 1.6典型超臨界流體 CO2的系統循環圖

1-3 文獻回顧

1-3.1 超臨界二氧化碳熱傳特性

由於超臨界二氧化碳在穿越臨界時，許多熱力、流力特性諸如密 度、焓、黏性等皆會有極大的變化，傳統適用於單相熱傳的經驗式， 如 Dittus-Boelter[2]所提出的和 Gnielinski Correlation[3]

Dittus-Boelter： 3 . 0 8 . 0

Pr

Re

023

.

0

_{b b} b

Nu

=

(適用於流體被冷卻情況，Re=10000~120000，Pr=0.7~120) Gnielinski： ) 1 (Pr ) 8 / ( 7 . 12 1 Pr ) 1000 )(Re 8 / ( 3 / 2 5 . 0 − + − = b b b b f f Nu

(適用於 Re=2300~1000000 ，Pr=0.5~2000，turbulent flow)

上二式皆是適用於熱力和流力性質變化不大的情況，但應用在二 氧化碳熱傳的預測上並不準確。原因在於，紐塞數受到雷諾數

(Reynolds number)和普蘭特爾數(Prandtl number)的影響，也就是與密 度以及熱傳係數的影響，又因二氧化碳在超臨界狀態時，其熱傳性能 以及流體性質會有很大的差異，如圖 1.8。

以圖 1.8 為例，圖中的點代表實驗數據，而線代表 Gnielinski Correlation 預測的數值，可以發現在低質量通率下預測較準，而在質 量通率漸漸增加時誤差逐漸加大，尤其在準臨界(Pseudo-critical)點附

近的誤差最大。 圖 1.8 可知二氧化碳在溫度與壓力分別為 40°C 與 9MPa，屬於界 流體，而其熱力性質以及流體特性由圖 1.9-圖 1.12 可以發現，其性 質在準臨界狀態皆有很大的變化。 二氧化碳在超臨界的狀態下，若以等壓的方式升降溫，當溫度到 達某一點的時候，等壓比熱將會達到一個極大值，在這點的性質稱之 為準臨界性質(Pseudocritical Property)，如準臨界溫度、準臨界密度等。 在此點附近，熱傳係數會有明顯的上升，這可歸咎於此點普蘭特爾數 的大幅增加所導致。而對於準臨界壓力及溫度，依據 NIST Refrigerants Database REFPROP[4]所提供的熱力表可以表示成和壓力間的關係： 3 5 . 2 2 0005608 . 0 01773 . 0 1657 . 0 124 . 6 6 . 122 P P P P T_pc =− + − + − 3 2 00009153 . 0 02901 . 0 233 . 4 6 . 272 P P P pc = + − + ρ 其中壓力單位為 bar，溫度和密度單位分別為℃及 kg/m3_。

在 Yoon et al.[5]的水平圓管實驗裡，在壓力 7.5~8.8MPa 之下對超 臨界二氧化碳做測試，結果發現熱傳性能在臨界點附近受到密度變化 的影響相當大，並做出了一熱傳修正式以描述其性質： 0.69 0.66

0.14 Re

Pr

b b b

Nu

=

_if

T

b

>

T

pc 0.05 1.6 0.013 Re Pr ( pc) b b b b Nu ρ ρ − = if

T

b

<

T

pc Son et al. [6]的實驗說明了在同一管截面上，由於準臨界附近的性

質變化相當大，因此壁面和管中心性質的差異上也會影響熱傳能力， 同時也提出了修正式： , 0.55 0.23 0.15 , Re Pr ( p b ) b b b p W c Nu c = if Tb >Tpc , 0.35 1.9 1.6 3.4 , Re Pr ( b ) ( p b ) b b b W p W c Nu c ρ ρ − − = if Tb <Tpc 熱通量和質通量的影響除了上述兩項因子外， Petrov et al. [7]則描述了管壁熱通量和質量通量會造成熱傳上的 差異，並提出了 Petrov Correlation： , , (1 0.001 )( p )n W O W p W c q Nu Nu G c = − 其中 NuO,W為： , 0.5 2/3 ( / 8) Re Pr 1.07 12.7( / 8) (Pr ) W W O W W f Nu f = + 上式中 n 值的大小由熱通量和質通量的比決定： 0.66 0.0004(q) n G = − if

c

p

<

c

p w, 0.9 0.0004(q) n G = − if

c

p

>

c

p w, p

c

_{的定義則為：} b w p b w

h

h

c

T

T

−

=

−

Liao et al. [8] 在水平圓管的實驗結果說明由於二氧化碳在準臨界附 近的密度變化非常大，因此由密度差造成的浮力效應將較為顯著，因 而對熱傳造成影響，而同時，作者也提到隨著管徑的減小，浮力效應 所帶來的響影將會減弱，這是因為 bouyancy parameter ( ₂ Re_b Gr )的大小 與管徑有關。並提出了紐塞數的修正式： 0.8 0.3 0.205 0.437 0.411 2 , 0.128 Re Pr ( ) ( ) ( ) Re p b b W W b W p W c Gr Nu c ρ ρ =

圖 1.9 超臨界下二氧化碳比熱變化

圖 1.11 超臨界二氧化碳黏滯性變化

1-3.2 流場計算之文獻

M.K. Bassiouny and H. Martin[9]先建立如下圖 1.13 之模型，並分 別設定流道入口與出口之控制體積，如圖 1.14 和圖 1.15： 並分別對入口與出口做分析，如下： (i) 首先對入口部分做質量守恆 c c dW AW A U A W Z dZ ρ =ρ +ρ _ + ∆ _   _{(式 1-1)} 且將 Lt Z N ∆ = 帶入上式， 整理後可得 c c AL dW u A n dZ = − (式 1-2) (ii) 對入口做動量守恆： 2 2 w c c c dP dW PA P Z A D Z A W Z AW A u W dZ τ π ρ dZ ρ ρ     −_ + ∆ _ − ∆ = _ + ∆ _ − +     (式 1-3) 在文獻中引入平均速度比 Wc W β = ，W A_c =u A_{c c} 由於文中將β視為常數，利用β 與流速之關係做整理可以知道 c c c A W W u A = − 將上式帶入(3)式，並與剪應力 2

(

)

8

w

W

f

τ

=

ρ

_{一起帶入整理} 可得

2 1 (2 ) 0 2 dP f dW W W dZ D β dZ ρ + + − = (式 1-4) (iii) 對出口做質量守恆 * * * * * c c dW A W A W Z A u dZ ρ = ρ _ + ∆ _+ρ   (式 1-5) 與(式 1-2)比較，並由入口推得知Ｕ做比較可得： * * c c c A L dW AL dW u A n dZ A n dZ = − = − 整理後可知： * * A W W A = (式 1-6) 而其中，有*即為出口之性質。 (iv) 對出口動量守恆 2 * * * * * * * 2 * * * * * * * * c c w c dP P A P Z A D Z dZ dW A W Z A W A U W dZ τ π ρ ρ ρ   −_ + ∆ _ − ∆     = _ + ∆ _ − +   _{(式 1-7)} 出口平均速度比可表示為 * * * c W W

β

= 而出口之剪應力為： *2 * *

(

)

8

w

W

f

τ

=

ρ

將上兩式帶入(式 1-7)可得： * * * *2 * * * 1 (2 ) 0 2 dP f dW W W dZ D β dZ ρ + + − = (式 1-8)

(v) 入出口合併 將(式 1-4)和(式 1-6)做比較，可得 1 𝜌 𝑑(𝑃−𝑃∗) 𝑑𝑑 + 1 2[ 𝑓 𝐷+ 𝑓∗ 𝐷∗� 𝐴 𝐴∗� 2 ]𝑊2_{− [(2 − 𝛽}∗_{) �}𝐴 𝐴∗� 2 − (2 − 𝛽)]𝑊𝑑𝑑_𝑑𝑑 = 0 (式 1-9) 而其中 * * 1 [ ] 2 f f D+D 在文獻中被視為是可以忽略的項。 (vi) 利用壓降公式 2 2 2 * 2 c t c c c u AL dW P P A n dZ ζ ρ ζ ρ    − = = − _{   }     (式 1-10) 將式 1-9 中 * P P− _{的部分代換掉} 其中 _c 1 _{t c to} e L C f C D ζ = + + + (vii) 無因次化 2 0 0 0 , , c , c u P W Z p w z W W

υ

W L

ρ

= = = = 可得 * * * ( ) [(2 ) (2 )] 0 d p p A dw w dz β A β dz − _{+ − − − =} (式 1-11) 將上式移項後可得二階常微分方程式 2 2 2 0 d w m w dz − = (式 1-12) 其中 _* * 2 (2 ) 2 2 2 [ ( ) 1] ( ) (2 ) c c nA m A A A β β β ζ − − = − − 為特徵值 這份文獻是假設在平均速度比β已知的條件下，並忽略出入口導 管摩擦損耗項，利用以上兩個假設簡化方程式。然而在真實情況下，

平均速度比會隨著分流區、匯流區位置的不同改變，在此假設為定值 並不恰當。此外，此結果有一前提假設即為流體密度在流場中為常數， 然而二氧化碳在超臨界狀態其密度變化相當劇烈，因此無法直接使用 此文獻來作為計算流場的工具。

圖 1.13 流體進入流道示意圖[9]

圖 1.14 入口控制體積示意圖[9]

1-3.3 熱傳計算之文獻

B. Prabhakara and P. KrishnaKumar [10]等人首先為了計算方便， 做了下列七個假設： (i) 熱力性質與壓力及溫度無關。 (ii) 每個流道的截面積是相同的。 (iii) 溫度轉換只包含有流道與流道之間，在出入口區並無發生溫度轉 換。 (iv) 與環境絕熱。 (v) 流體在流道中屬於塊狀流且均勻分布於流道中。 (vi) 在流道間不均勻的流體分佈可以帶入。 (vii) 平板很薄可假設成無側向熱傳。 而本文獻首先建立如圖 1.16 之模型： 圖 1.16 熱傳控制體積示意圖[10] 接著將控制體積設定在流道中，對第一個流道做熱傳分析，可得

到如下之方程式： For fluid 1： ( )₁ ( ) ( ₁ ) 2 2 i i w i w p wi i wi i dT h A h A mC T T T T dx = L − + L + −  (式 1-13) For fluid 2： 1 2 1 ( ) ( 1) ( ) ( ) 2 2 i i i i p wi i wi i dT h A h A mC T T T T dx L L − +   = − ×_ − + − _    (式 1-14) 而在板片的熱傳部分即為： 1 1 ( ) ( ) 0 2 2 i w i w i wi i wi h A h A T T T T L L − − − + − = (式 1-15) 將方程式(14)無因次化並整理可得到： 1 1 1 1 1 1 1 1 1 ( 1) 2 i i i i i i i i i i i i i i i i i h t h h dt NTU dX _{h h h} t t h h h h h h − − − − + + − + +    +  ₊       = − _{×  }     _{ + −  + }  _ + + _{ } + _    _{(式 1-16)} 其中 2( ) w p uniform hA NTU mC = _且 1, 2, 1, in in in T T t T T − = − 而邊界條件為： X=0, ti=0 i=1, 3, 5,...,N for x=0,ti =0 X=1, ti=1 i=2, 4, 6,...,N − for 1 x=1,ti =1 此文獻為在假設流體熱力性質與溫度無關的情形下，所推導出來 的結果，然而在真實情況下工作流體的熱力性質將會隨著溫度的變化

第二章

物理模式

2-1 板式熱交換器幾何尺寸

板式熱交換器的幾何圖形如圖 2.1 所示。

圖 2.1 常見板式熱交換器幾何參數其及定義 幾何參數定義如下：

(ii) 板片長度，𝐿_𝑝 (Vertical distance between ports) (iii) 板片寬度，w (Plate width)

(iv) 平均流道間距，b (Mean flow channel gap) (v) 傾斜角，β (Chevron Angle)

傾斜角β定義為山型紋波紋方向與主流方向的夾角。一般而言，

當傾斜角越大時，其熱傳效果及壓降也會隨之增加。 (vi) 流道截面積，𝐴_𝑐 (Channel flow area)

流道截面積𝐴_𝑐的定義為：

𝐴𝑐=b*w

(vii) 流道等效直徑，𝐷_𝑒 (Channel equivalent diameter)

在無因次分析中，通常使用水力直徑（𝐷_ℎ）來表示。水力直 徑的定義為： 𝐷ℎ =4 ∗ 流道截面積 濕周長 而代入板片之幾何尺寸(如圖)得： 𝐷ℎ =_{2 ∗ (b + w) =}4 ∗ 𝐴𝑐 _{2 ∗ (b + w)}4 ∗ 𝑤 ∗ 𝑏 對一般分析而言，由於 W 遠大於 b(W » b)，所以 𝐷ℎ ≒ 4 ∗ 𝑤 ∗ 𝑏_{2 ∗ (𝑤) = 2𝑏} 為了避免將其與一般水力直徑混淆， 我們特地將它稱為等效直 徑(𝐷_𝑒, equivalent diameter)。

2-2 分析假設與統御方程式

(i) 假設為連續流體。

(ii) 超臨界狀態的二氧化碳為牛頓流體(Newtonian fluid)，且黏

滯係數為等方向性。 (iii) 工作流體在板式熱交換器中，流體密度會隨溫度與壓力而 改變。 (iv) 假設為一維流動；即在同一截面下為均勻流。 (v) 由於板片的厚度很薄，因此假設板片在側向不發生熱傳 導。 (vi) 由於入口/出口端長度相對於板片長度來甚小，且管壁摩擦 係數亦遠小於板片通道，因此在此不考慮入/出口端摩擦壓 降。 (vii) 忽略重力效應。 (viii) 板片與外界無發生熱傳效果；即邊界設為絕熱。

統御方程式分別為連續方程式(Continuity equation)、動量方程式 (Momentum equation)及能量方程式(Energy equation)。

連續方程式(Continuity equation)： 𝜕𝜕 𝜕𝜕 + 𝜕 𝜕𝑥𝑗�𝜕𝑢𝑗� = 0 動量方程式(Momentum equation)： 𝜕𝜕𝑢𝑖 𝜕𝜕 + 𝜕 𝜕𝑥𝑗�𝜕𝑢𝑗𝑢𝑖� = − 𝜕𝜕 𝜕𝑥𝑖 + 𝜕𝜕 𝜕𝑥 𝜎𝑖𝑗 能量方程式(Energy equation)： 𝜕𝜕𝜕 𝜕𝜕 + 𝜕 𝜕𝑥𝑗�𝜕𝑢𝑗𝜕 − 𝜕𝑢𝑗� = − 𝜕 𝜕𝑥𝑗𝑞𝑗 + 𝜕 𝜕𝑥𝑗(𝑢𝑖𝜎𝑖𝑗)

第三章 數值計算模式與分析方法 本章將介紹及說明本研究中所使用的方法以及其步驟。第一節為 模型的介紹。將經簡化過後的板式熱交換器切割成若干個格點，以及 模型相關的何參數等。第二節為方程式的建立。參考 Bassiouny and Martin[9] 的 U 型熱交換器分流理論，並於分流區與匯流區在動量方 程式的部分做修正，並建立修正後的方程式。第三節為解題的流程。 將方程式分成流場及溫度場兩大部分，先解流場，再解溫度場，最後 將兩者的結果反覆互相迭代，以求得板式熱交換器內部真正的流動情 形。第四節為牛頓法(Newton method))的介紹。由於根據第二節所建 立的方程式，為一非線性聯立方程組，因此在此使用牛頓法來求解此 方程組。同時，針對牛頓法，介紹了牛頓法在應用上的缺點，並針對 缺點引入牛頓下降法來改良，以增加程式的收斂範圍。第五節為程式 流程。此節先分別介紹了流場與溫度場的程式流程，並於最後說明流 場及溫度場整合後的程式流程以及收斂條件等。 綜合上述，本研究在數值方面的計算過程為，首先根據使用者所 給定的板片數，包含熱通道與冷通道的數量，以及板片的幾何形狀等， 配合 Zahid [11]於 2003 所提及的板式熱交換器熱傳與壓降相關經驗式， 可以分別解出冷側及熱側流場，再將冷側及熱側流場進行熱交換，緊 接著重新更新冷側及熱側流場的速度、壓力及其它相關熱力性質，並

再次代入溫度場中，重覆步驟直到流場的速度及壓力等變數與上一次 計算完溫度場後相同時，則程式結束。 3-1 模型介紹 首先根據 2-2 節的假設，流體為連續流體且流動型式為一維，忽略重 力效應所造成的影響，且已知入口壓力及溫度並與外界絕熱。接著根 據圖 1.13U 型板熱交換流場分佈，可建立流場模型。 由於冷熱流體並不互相接觸，且行走於相異通道，為了避免混肴， 在此將冷側與熱側通道的模型配置圖分開來討論。 冷側(工作流體為水)模型如圖 3.1 所示，其中左下方為入口；由於 為 U 型配置，故出口設於左上方。在模型上方及下方分別為匯流區 以及分流區；而在中間則為板片通道，此部分佔了流場約 80%以上的 體積，而根據 2-2 假設，絕大部分的壓降以及所有的熱傳現象，也都 發生在板片通道之間。 熱側(工作流體為二氧化碳)模型如圖 3.2 所示。 由於流動方向為逆向，故熱側模型的入口在左上角而出口在左下 角，上下及下方分別為分流區及匯流區，其餘的配置則與冷測模型相 同。

圖 3.1 冷側數值模型配置圖

3-2 方程式建立

本文方程式的建立為以 Bassiouny and Martin[9] U 型板式熱交換 器分流理論為基礎。Bassiouny and Martin[9] 的所解出的理論解雖然 可以符合大部分流體在板式熱交換器的分流情形，二氧化碳在實際應 用上仍有一些較不方便的地方： (i) [9]在分流區(匯流區亦同)動量方程式的部分，設定平均速度 比𝛽為定值(如式 3-1)，由於固定𝛽代表分流區及匯流區的損 失係數(loss coefficient)為定值，然而在真實情況下，損失係 數會因為流量不同而有所改變，因此並不能假設為定值。 方程式推導如下： PA − �P +𝑑𝑃_𝑑𝑑∙ ∆Z� A = ρA �W +𝑑𝑑_𝑑𝑑 ∙ ∆Z�2− 𝜕𝐴𝑊2_{+ 𝜕𝐴𝑢} 𝑐β𝑊 (式 3-1) 其中 β = 𝑑𝐶 𝑑 = 𝑢𝑐𝐴𝑐 𝐴 𝑑 𝑢𝑐 =𝛽𝑑𝐴_𝐴𝑐 若將(式 3-1)加入離散的概念，整理後可以得到 𝑃𝑖A − 𝑃𝑖+1𝐴 = 𝜕𝐴𝑊𝑖+12− 𝜕𝐴𝑊𝑖2+ 𝜕𝐴𝑐𝑢𝑐β𝑊𝑖 (式 3-2) 移項可得 𝑃𝑖A − 𝑃𝑖+1𝐴 = 𝜕𝐴𝑊𝑖+12− 𝜕𝐴(1 − 𝛽2)𝑊𝑖2 (式 3-3)

因此β2_{項相當於分流區損失係數(loss coefficient)，在文獻中作者} 將β假設為定值，但實際應用上，假設為定值並不恰當，可能導致無 法反應出二氧化碳板式熱交換器內部真實流動情形。 (ii) 流動的性質如密度、黏滯係數等在文獻中皆假設為常數。 但當二氧化碳位於超臨界區時，性質在溫度及壓力變化非 常劇烈，假設為常數則會容易造成相當大的誤差。 以冷側為例，針對以上兩點，建立並修正後的方程式及控制體積 如下：

圖 3.3 分流區控制體積示意圖

圖 3.4 匯流區控制體積示意圖

圖 3.5 流道控制體積示意圖

方程式如下： (1) 分流區： (i) 質量守恆 1 1 ,1 i i i i in

AV

in

A u

c i in

AV

in

ρ

=

ρ

+

ρ

_{+ +} (ii) 動量守恆 1 1 2 2 2 i i i i i in in in in L in

P A

P

A

ρ

AV

ρ

AV

K

ρ

AV

+ +

−

=

−

+

(2) 匯流區： (i) 質量守恆 1 1 , i i i

out

AV

out

A u

c i N out

AV

out

ρ

=

ρ

+

ρ

_{+ +}

(ii) 動量守恆

1 1 1

2 2 2

i i i i i

out out out out c out

P

A

P

A

ρ

AV

ρ

AV

K

ρ

AV

+

−

=

−

+

+

+ (3) 流道： (i) 質量守恆 , , , 1 , 1

0

i j

A

c

u

i j i j

A u

c i j

ρ

−

ρ

_{+ +}

=

(ii) 動量守恆 2 2 , , 1 , , 1

0

i j c i j c w c i j c i j

P A

−

P

₊

A

−

τ

P

l

−

ρ

A

u

+

ρ

A u

₊

=

其中 , 2 ,

f

ρ(

)

8

i j i w c j

u

τ

=

除此之外，根據[9]，本研究還考慮轉入(出)流道的能量損失(Turning loss)如下： (i) _,1

(

)

2 ,1 ,1

A

1

0

i i in ch c t i i

P A

−

P

+ +

C

ρ

Au

=

(ii) _,

(

)

2 , ,

1

0

i N i ch c out to i N i N

P

A

−

P

A

+ +

C

ρ

Au

=

另外在熱場的部分，由於從入口到流道之間的距離非常短，因此 假設熱傳僅發生在板片通道的地方，至於方程式的建立，則參考[10] 的模型，所建立的控制體積如下圖所示：

圖 3.6 冷側熱場控制體積示意圖

方程式如下： (i) 水側板片通道熱傳 m𝐶𝑝𝑝 𝑖�𝑇𝑝 𝑖,𝑗+1 − 𝑇𝑝 𝑖,𝑗� = �_𝐴 1 𝑤ℎ𝑝 𝑖,𝑗 + ∆𝑥 𝑘𝐴𝑤 + 1 𝐴𝑤ℎ𝐶 𝑖,𝑗� −1 𝐿𝐿𝑇𝐷 + �_𝐴 1 𝑤ℎ𝑝 𝑖,𝑗 + ∆𝑥 𝑘𝐴𝑤 + 1 𝐴𝑤ℎ𝐶 𝑖+1,𝑗� −1 𝐿𝐿𝑇𝐷 (ii) 二氧化碳側板片通道熱傳 m𝐶𝑝𝐶 𝑖�𝑇𝐶 𝑖,𝑗+1 − 𝑇𝐶 𝑖,𝑗� = �_𝐴 1 𝑤ℎ𝑝 𝑖,𝑗 + ∆𝑥 𝑘𝐴𝑤 + 1 𝐴𝑤ℎ𝐶 𝑖,𝑗� −1 𝐿𝐿𝑇𝐷 + �_𝐴 1 𝑤ℎ𝑝 𝑖+1,𝑗 + ∆𝑥 𝑘𝐴𝑤 + 1 𝐴𝑤ℎ𝐶 𝑖,𝑗� −1 𝐿𝐿𝑇𝐷

3-3 解題流程 將 3-2 節提到的方程式聯立後，分別可建立流場場及溫度場矩陣 如下： 流場：

⎣

⎢

⎢

⎢

⎢

⎢

⎡

{1} {1} 0 0 0 0

₀

_{{1} 0 0 0 0}

0 {1} {1} 0 0 0

{𝑉

𝑖𝑖

} 0 0 {1} 0 0

0 {𝑈

𝑐ℎ

} 0 0 {1} 0

0 0 {𝑉

_𝑜𝑜𝑜

} 0 0 {1} ⎦

⎥

⎥

⎥

⎥

⎥

⎤

×

⎣

⎢

⎢

⎢

⎢

⎢

⎡

_{𝑈

{𝑉

𝑖𝑖

}

𝑐ℎ

}

{𝑉

𝑜𝑜𝑜

}

{𝑃

𝑖𝑖

}

{𝑃

𝑐ℎ

}

{𝑃

𝑜𝑜𝑜

}⎦

⎥

⎥

⎥

⎥

⎥

⎤

=

⎣

⎢

⎢

⎢

⎢

⎡

𝑉

₀

𝑖𝑖

0

𝑃

_𝑖𝑖

0

0 ⎦

⎥

⎥

⎥

⎥

⎤

溫度場：

�

�𝑚

𝑝

𝐶

𝑃𝐻

� �𝑈

(𝑝+𝐶)

𝐴

𝑑

�

�𝑚

𝐶

𝐶

𝑃_𝐶

� �𝑈

(𝑝+𝐶)

𝐴

𝑑

�

� × �

{𝑇

_𝑝

}

{𝑇

𝐶

}� = �

𝑇

𝑝_in

𝑇

𝐶_in

�

值得一提的是，速度場及溫度場為一非線性聯立方程組，無法藉 由直接求解法來求解，本文在此利用牛頓法來求解，而關於牛頓法將 在下一節有更進一步的介紹。 註：{𝑉_𝑖𝑖}代表由𝑉_𝑖𝑖1, 𝑉_𝑖𝑖2, 𝑉_𝑖𝑖3…𝑉_𝑖𝑖𝑛等所構成的集合，其他參數亦同。 而{1}僅表示{𝑉_𝑖𝑖}集合中𝑉_𝑖𝑖𝑖所對應的係數不為零。 至於在等號右邊V_in、P_in、T_Hin、T_Cin等為入口速度、壓力溫度等 已知條件。

3-4 牛頓法的介紹 牛頓法[12],[13](Newton method)一般又稱為切線法，為非線性方 程式相當重要的一種方法，許多解非線性方程式所使用的迭代法，皆 以牛頓法為基礎而發展而來，以下將對牛頓法有約略的介紹。 假設有一非線性方程式f(𝑥) = 0, 若存在𝑥∗使得f(𝑥∗) = 0，則稱𝑥∗ 為f(𝑥) = 0的一解。 今假設𝑥₀為f(𝑥)的一近似解，將f(𝑥)在𝑥₀處做泰勒展開(Taylor expension)，可得 𝑓(𝑥) = f(𝑥0) + (𝑥 − 𝑥0)𝑓′(𝑥0) +(𝑥 − 𝑥0) 2 2! 𝑓′′(𝑥0) … 若忽略高次項，則上式可改寫為 𝑓(𝑥) = f(𝑥0) + (𝑥 − 𝑥0)𝑓′(𝑥0) 移項可得 x = 𝑥0−_𝑓′(𝑥𝑓(𝑥0) 0) 此為牛頓法的迭代式。

3-4-1 牛頓-拉弗森法(Newton-Raphson method) 將牛頓法推廣到多維度，則稱為牛頓-拉弗森法(Newton-Raphson method)。 假設有一非線性方程組 � 𝑓1(𝑥1, 𝑥2, … 𝑥𝑖) = 0 𝑓2(𝑥1, 𝑥2, … 𝑥𝑖) = 0 ⋮ 𝑓𝑖(𝑥1, 𝑥2, … 𝑥𝑖) = 0 令�𝑥₁(𝑘), 𝑥₂(𝑘), … 𝑥_𝑖(𝑘)�為方程組的一組近似解，將原方程組對泰勒 級數展開並忽略高次項，則可得一近似方程組 ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧𝑓1�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� + � 𝜕𝑓1�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� 𝜕𝑥_𝑖 �𝑥𝑖 − 𝑥𝑖𝑘� 𝑖 𝑖=1 = 0 𝑓2�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� + � 𝜕𝑓2�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� 𝜕𝑥𝑖 �𝑥𝑖 − 𝑥𝑖 𝑘_� 𝑖 𝑖=1 = 0 ⋮ 𝑓𝑖�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� + � 𝜕𝑓𝑖�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� 𝜕𝑥𝑖 �𝑥𝑖 − 𝑥𝑖 𝑘_� 𝑖 𝑖=1 = 0 將其寫為矩陣型式，則可改寫為 ⎣ ⎢ ⎢ ⎢ ⎢ ⎡𝑓1𝑥1 𝑓1𝑥2 𝑓1𝑥3 … 𝑓_2𝑥1 𝑓2𝑥2 𝑓2𝑥3 … 𝑓3𝑥1 𝑓3𝑥2 𝑓3𝑥3 … ⋮ ⋮ ⋮ … 𝑓𝑖𝑥1 𝑓𝑖𝑥2 𝑓𝑖𝑥3 …⎦ ⎥ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎢ ⎡∆𝑥1𝑘 ∆𝑥2𝑘 ∆𝑥3𝑘 ⋮ ∆𝑥_𝑖𝑘_⎦⎥ ⎥ ⎥ ⎥ ⎤ = − ⎣ ⎢ ⎢ ⎢ ⎡_𝑓𝑓1 2 𝑓3 ⋮ 𝑓𝑖⎦ ⎥ ⎥ ⎥ ⎤

其中等號右側 𝑓1 = 𝑓1�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)�, 𝑓2 = 𝑓2�𝑥1(𝑘), 𝑥2(𝑘), … 𝑥𝑖(𝑘)� … 等。 等號左側𝑓_1𝑥1 = 𝜕𝑓1 𝜕𝑥1, 𝑓1𝑥2 = 𝜕𝑓1 𝜕𝑥2… … 𝑓𝑖𝑥𝑗 = 𝜕𝑓𝑖 𝜕𝑥𝑗…。此一矩陣又稱為 Jacobian 矩陣。 而𝑥_𝑖𝑘表示𝑥_𝑖在第 k 次迭代時的猜值，而∆𝑥_𝑖𝑘即為修正值，即 𝑥₁(𝑘+1) = 𝑥₁𝑘 _{+ ∆𝑥} 1𝑘 , 𝑥2(𝑘+1) = 𝑥2𝑘 + ∆𝑥2𝑘 , 𝑥3(𝑘+1) = 𝑥3𝑘 + ∆𝑥3𝑘… … 而在矩陣中間偏微分的部分本研究則採用中央插分法來表示： (𝑓1)𝑥1 = 𝑓1(𝑥1+ 𝛿, 𝑥2, 𝑥3, … ) − 𝑓1(𝑥1− 𝛿, 𝑥2, 𝑥3, … ) 𝛿 (𝑓1)𝑥2 = 𝑓1(𝑥1, 𝑥2+ 𝛿, 𝑥3, … ) − 𝑓1(𝑥1, 𝑥2−𝛿, 𝑥3, … ) 𝛿 ⋮ (𝑓𝑖)𝑥𝑗 = 𝑓_𝑖�𝑥1, 𝑥2, 𝑥3, … 𝑥𝑗 + 𝛿, … � − 𝑓𝑖�𝑥1, 𝑥2−𝛿, 𝑥3, … 𝑥𝑗 + 𝛿, … � 𝛿 由於牛頓法為局部收斂，當初始猜值在根附近時，收斂速度相當 快速，而當猜值與根的距離相對較遠時，則容易會有無法收斂的情形 發生，而實際上，並無法保證初始猜值的位置，這使得牛頓法應用在 解較複雜的非線性方程組時，有相當大的困難。

3-4-2 牛頓下降法 為了解決牛頓法過於依賴猜值所造成的問題，牛頓下降法便是針 對未能給定較佳的初始猜值時，所對牛頓法的一種修正。此法乃是藉 由增加牛頓法的收斂範圍來降低牛頓法對猜值的依賴性，此法在應用 上非常的重要。 若將迭代式 𝑥�(𝑘+1) _{= 𝑥�}(𝑘)_{− �∇𝑓̃}(𝑘)_�−1_{𝑓̃�𝑥�}(𝑘)_� 修正為 𝑥�(𝑘+1) _{= 𝑥�}(𝑘)_{− 𝑤�∇𝑓̃}(𝑘)_�−1_{𝑓̃�𝑥�}(𝑘)_{� , 0 < 𝑤 ≤ 1} 其中 w 為一參數，選擇適當的 w 滿足 �𝑓̃�𝑥�(𝑘+1)_{�� < �𝑓̃�𝑥�}(𝑘)_�� 一般而言，w 依次為 1 ,1 2 , 1 22,…，此稱為牛頓下山法。 其中 𝑓̃ = (𝑓1, 𝑓2, … 𝑓𝑖) 𝑥� = (𝑥1, 𝑥2, … 𝑥𝑖)

3-5 程式流程 在程式開始之前，首先必須把水及二氧化碳的熱力性質繪出 [4](如圖 1.9-圖 1.12)。接著找出可以描述此一方程式曲線並建構在程 式中，以便程式在迭代時更新熱力性質使用。另外，根據 3-2 節所述， 程式開始之前，也必須將分流區及匯流區不同截面積及流量之下所對 應的損失係數[14] (如圖 3.8、圖 3.9)建立於程式中。 圖 3.8 分流區不同流速下主流方向所對應的損耗係數[14]

流程圖(流場)： 圖 3.10 流場程式流程圖 ② 輸入板片內速度 及壓力(猜值) ⑥利用牛頓下降法計 算出新的速度及壓力 值代回③中 ⑦當⑤符合時，表示連續及動量 方程式均收斂，猜值即為解。 程式結束 ① 程式開始 ④由輸入的板 片數量建立方 程式(f1~fn) ⑤_{∑ |𝑓}𝑖 _𝑖| 𝑖=1 ≤收斂標準? NO ③由速度及壓力計算內 部流體的性質 YES

流程圖(溫度場)： 圖 3.11 熱場程式流程圖 ② 輸入板片內熱通 道及冷通道流速 ⑤利用牛頓法求解溫 度場 ⑧程式結束 ① 程式開始 ④由熱場方程式建 立矩陣 ③由速度計算 Nusselt Number 等熱力性質，輸 入溫度猜值 ⑥判斷是否滿足能量方程式 ⑦利用牛頓下降法計 算出新的溫度猜值代 回③中

②輸入板片內速 度、壓力及初始溫 度，並得到流體性質 ③計算水的流場 ④計算二氧化碳的流場 ⑤計算熱交換器內之溫度場 ⑥判斷流場速度與壓力與上 次計算，其差值是否收斂標準 ⑧程式結束 ① 程式開始 流程圖(流場及熱場整合) 圖 3.12 整合後程式流程圖 NO ⑦更新熱場溫度並代入②中 重覆此動作直至⑥符合收斂 標準 YES

第四章

研究成果與討論

本章將針對程式所模擬出來的結果分為兩大部分進行討論。第一 部分為模擬結果的正確性；藉由對照的方式，將水對水板式熱交換器 模擬出來的結果，針對流場與熱場的部分與現有文獻比較，以驗證程 式結果的準確性。第二部分則為二氧化碳板式熱交換器模擬結果；由 改變模擬參數的方式(如入口溫度、速度、壓力、板片數量等)，來探 討各個參數對熱傳結果之效應。

4-1 格點測試 由於格點的數量將會直接影響模擬結果的準確性以及計算時間， 因此在比較不同的參數前，必須先確保格點的數量並不會對模擬結果 造成太大的影響，因此在這個部分先分別就水對水以及水對二氧化碳 板式熱交換器的程式模擬結果做測試。 所模擬板片幾何參數及結果如下： 幾何及入口條件 冷側流道數(n) 10 熱側流道數(n) 9 板片間距(mm) 2.9 板片厚度(mm) 0.8 板片長度(mm) 600 板片寬度(mm) 218 入口直徑(mm) 70 流道特效直徑(mm) 5.8

熱容比(heat capacity ratio) 1 水對水板式熱交換器入口條件 冷側入口流速(m/s) 0.535 熱側入口流速(m/s) 0.522 水對二氧化碳板式熱交器入口條件 冷側(水)入口流速(m/s) 熱側(二氧化碳)入口流速(m/s) 12.816 熱側(二氧化碳)入口壓力(MPa) 8 經驗式，ref[16] 摩擦因子 0.206 1.441 Re f = ⋅ − Nusselt Number 0.663 0.333 0.17 0.3 Re Pr ( / _m) Nu= ⋅ ⋅ ⋅ m m

圖 4.1 及圖 4.2 為水對水以及水對二氧化碳板式熱交換器在熱容 比為 1 時，所得之結果，其中橫軸為單一板片所切割的格點數，而縱 軸當前格點數之出口溫度與格點數為 25 時，出口溫度的差異。 , 25, n out out

T

T

T

∆ =

−

其中Tn out, 格點數為 n 時熱側出口溫度，T25,out為格點數為 25 時， 熱側出口溫度。

而 C 為熱容比(heat capacity ratio) C=(mc_p)_min / (mc_p)_max

由於本研究皆令 Cmin在熱側，因此，C 又可表示為 ( p) / (H p C) C= mc mc 其中在此需特別注意，二氧化碳比熱為溫度與壓力的函數，不為 定值，因此本研究將二氧化碳的熱容 2 CO C 定義為二氧化碳入口狀態的 熱容值。 若將圖 4.1 及圖 4.2 相比，則可以發現水對二氧化碳板式熱交換 器在格點數較少時，出口溫度的差異相較於水對水板式熱交換器來得 大，此原因乃是因為當格點數較少時，由於二氧化碳的性質變化較大， 因此較少的格點數較無法取得板式熱交換器內部真正的流體性質，因 此誤差較大。而另一方面，當格點數增加時，則會相對的增加程式的 運算時間(當格點數增加為 2 倍時，計算時間約增加 23_{倍)，因此，在} 雙重在考量之下，本研究以 N=15 為最佳的格點數。

5 10 15 20 25 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 ∆T Number of grid C=1 water to water PHX 圖 4.1 水對水熱交換器格點測試結果 5 10 15 20 25 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 ∆T Number of grid C=1 water to CO₂ PHX 圖 4.2 水對 CO2熱交換器格點測試結果

4-2 模擬驗證 4-2.1 流場驗證 由於目前針對二氧化碳應用在板式熱交換器的相關文獻較缺乏 的緣故，在此引用水對水板式熱交換器的相關文獻結果，來驗證程式 結果的準確性，接著將改變工作流體為二氧化碳，以進行二氧化碳板 式熱交換器的設計。 首先在流場的部分，引用 Rao et al. [15]的的實驗結果來驗證，所 使用的板片以及幾何參數如圖 4.3。而圖 4.4 為實驗所得的數據以及 回歸後摩擦因子與對應雷諾數之間的關係圖。 根據[15]，所建立模型其他相關幾何參如下： 水流道個數：10 流道格點數：15 板片間距：2.4 mm 板片厚度：0.6 mm 流道等效直徑：4.8 mm

圖 4.3 文獻所使用板片幾何圖形及參數[15]

圖 4.5 壓降與流道編號關係圖

圖 4.5 為在相同板片下，平均壓降相同時，Rao el al.[15]實驗結果 與程式結果比較圖。圖 4.5 顯示模擬結果在靠近入口處的壓降值較大， 而在下游處則會小，整體而言，壓降分佈曲線較為陡峭。此結果亦顯 示了在各個流道中的流速分佈較不均勻。

圖 4.6 無因次化後速度與流道編號關係圖 其中 𝑣_𝑐 = 𝑉𝑐ℎ

𝑉𝑐ℎ𝑚

圖 4.6 中實心點為 Rao el al.[15] 根據圖 4.5 的壓降結果，並參考 Bassiouny and Martin [9] 的理論模型，所推得的速度分佈，而空心點 為則本研究程式模擬結果。圖 4.6 的結果顯示程式結果在板片中流速 的預測上，誤差約為 10%左右，而此結果符合圖 4.5 壓降分佈的預測 結果。

在流場的部分，由壓降與流速的模擬結果顯示，本研究開發出的 程式，在預測板式熱交換器的流場上，具有一定的準確性。

4-2.2 溫度場驗證 在水對水板式熱交換器溫度場的部分，本文在此引用 Gherasim et al. [16]的水對水板式熱交換器模擬結果來進行驗證。 幾何及入口條件 冷側流道數(n) 10 熱側流道數(n) 9 板片間距(mm) 2.9 板片厚度(mm) 0.8 板片長度(mm) 600 板片寬度(mm) 218 入口直徑(mm) 70 流道特效直徑(mm) 5.8

熱容比(heat capacity ratio) 1 水對水板式熱交換器入口條件 冷側入口流速(m/s) 0.535 熱側入口流速(m/s) 0.522 熱水入口溫度(K) 353 冷水入口溫度(K) 303 經驗式，ref[16] 摩擦因子 0.206 1.441 Re f = ⋅ − Nusselt Number _{0.3 Re}0.663 _Pr0.333 _{( / )}0.17 m Nu= ⋅ ⋅ ⋅ m m

圖 4.7 為 Gherasim et al. [16]的模擬結果與本研究程式的模擬結果 比較圖。其中橫軸為流道位置。 x z L = 其中 x=0 為流道入口處，而 x=L 則為流道出口贏。 縱軸為無因次化後的溫度。 C H C T T T T θ = − − Gherasim et al. [16]的結果指出，在冷側流道的部分，溫度分佈大 致可以分為三個溫度較接近的區間，分別為 1,10 流道 (冷側第一個以 及最後一個流道)、2,9 (冷側第二個以及倒數第二個流道)、以及其他 的流道。而熱側則可分為 1,9 流道(熱側第一個及最後一個流道)及其 他流道。 本研究模擬結果同樣也符合文獻中所提及的溫度分佈情形，在熱 側出口處溫度差異在 1.5 度以內(約 3%)。 由流場以及溫度場驗證結果顯示，本研究所建立之模型，可應用 於板式熱交換器的模擬中。

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

θ

z

Hot Channel - Gherasim et al. Cold Channel- Gherasim et al. Hot Channel - present research Cold Channel- present research

圖 4.7 溫度場結果與文獻比較

1, 10 channel 2, 9channel 1, 9 channel

Other hot channels

4-3 流場結果 首先，在此先考慮流場為等溫的情形(即無發生熱交換)，在不同的入 口雷諾數下，各個流道的分佈情形。 在討論流速分佈之前在此先定義一無因次化的標準差Φ來做為指 標。Φ越大，表示流場越不均勻。Φ定義如下 2 1 1 ( ) n i i n = β β Φ =

∑

− 圖 4.8 及圖 4.9 為討論水在等溫時，不同的入口雷諾數對板片間 流動不均勻性的影響。其中圖 4.8 流道數量為 10，而圖 4.9 流道數量 為 40。 由圖 489 可以發現當入口雷諾數增加時，板片間的不均勻性將隨 之增加，此結果乃是因為當流速增加時，在分流區以及匯流區的摩擦 阻力亦隨之增加，因此造成越往下游處，其速度分量將會越小。 圖 4.9 為將流道數量由 10 增加為 40 時，流場所得之結果。若將 圖 4.8 與 4.9 比較，可以發現當流道數量增加時，其流動不均性將會 增加，而且其幅度相當快速。

0 1 2 3 4 5 6 7 8 9 10 0.088 0.090 0.092 0.094 0.096 0.098 0.100 0.102 0.104 0.106 0.108 0.110 0.112 0.114 0.116 0.118 0.120 β_i=m_i/m Number of channel Re500 Φ=0.003425 Re1000 Φ=0.003776 Re2000 Φ=0.004161 Re5000 Φ=0.004729 Re10000 Φ=0.005208 Re100000 Φ=0.007151 average 圖 4.8 等溫時不同入口雷諾數對流速分佈的影響(n=10) 0 5 10 15 20 25 30 35 40 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 0.065 0.070 0.075 β_i=m_i/m Number of channel Re500 Φ=0.010486 Re1000 Φ=0.0112899 Re2000 Φ=0.012133 Re5000 Φ=0.013303 Re10000 Φ=0.014236 Re100000 Φ=0.017561 average 圖 4.9 等溫時不同入口雷諾數對流速分佈的影響(n=40)

圖 4.10、圖 4.11 為在水相同的質量流率下，當改變入口溫度且無 發生熱交換時時，流場在不同流道數量下之分佈情形。 由圖可知，當質量流率固定時，若提高入口溫度(流動過程仍為等溫)， 則流道間的流動不均勻性也會增加，造成此結果的原因為當入口溫度 升高時，則流體黏滯力降低的緣故(當溫度由 298K 提升至 358K 時， 流體黏滯力由 0.000893 Pa-s 降低至 0.000334 Pa-s)，而黏滯力降低則 造成雷諾數的增加，因而略微增加了流場分佈的不均勻性。但整體而 言，水因入口溫度不同而對流場所造成的影響有限。

0 1 2 3 4 5 6 7 8 9 10 0.088 0.090 0.092 0.094 0.096 0.098 0.100 0.102 0.104 0.106 0.108 0.110 0.112 0.114 0.116 0.118 0.120 β_i=m_i/m Number of channel T=298K Re2000 Φ=0.004161 T=328K Re3537 Φ=0.004506 T=358K Re5352 Φ=0.004774 average 圖 4.10 質量流率相同、溫度不同流場分佈情形(n=10) 0 5 10 15 20 25 30 35 40 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 0.065 0.070 0.075 β_i=m_i/m Number of channel T=298K Re2000 Φ=0.012288 T=328K Re3537 Φ=0.013019 T=358K Re5352 Φ=0.013566 average 圖 4.11 質量流率相同、溫度不同流場分佈情形(n=40)

4-4 溫度場結果 4-4.1 溫度場結果分析 在溫度場首先以模擬水-二氧化碳在二氧化碳入口壓力為 8MPa， C=1 時的流場及溫度場結果。其結果與其它參數如下： 水側入口流速：0.518 m/s (流量為 2L/s) 水入口溫度：283 K 二氧化碳入口溫度：370 K 二氧化碳出口壓降：179.3 kPa 所使用的板片模型與 Gherasim et al. [16]相同。

圖 4.12-圖 4.17 為當水(冷側)流量為 2 L/s 時、C=1 時，二氧化碳 入口壓力為 8 MPa 時，程式模擬所得流場及熱交換結果。在圖 4.12 中，水在第一個流道以及最後一個流道處，其溫度約為 312 K、而在 中間流道則約為 332 K，中間相差了 20 K，而造成此現象的原因為水 側在第一個板片與最後一個板片處，與二氧化碳流道僅有一面接觸， 熱傳效果較差，因而出口溫度較低。而在二氧化碳側(圖 4.13) 可以 發現在各個流道的溫度分佈皆相當平均，唯在第一個通道與第九個通 道處較平均溫度低了兩度，雖然在第一個與第九個通道處，二氧化碳 兩側皆與水接觸，但由於在最外側水溫度較低，溫差較大而造成熱側 通道出口溫度較低(熱傳較佳)。此結果亦可以參考圖 4.16 及圖 4.17 之熱傳量分佈曲線。 在圖 4.14 及 4.15 中，中間流道的速度分佈方式皆為隨著流道編 號的增加而遞減，此部分與 4.3 節流場結果類似。在圖 4.14 中，水在 第一個流道與最後一個流道處，由於溫度較低，黏滯係數較大的影響， 因此速度相對相鄰流道將會偏低。然後在圖 4.15 中，二氧化碳雖然 在第一個流道以及最後一個流道處溫度皆偏低，但其入口流速反而相 鄰流道來的高，此現象可用準臨界性質來解釋。在二氧化碳第一個以 及最後一個流道中，由於流道下游處溫度偏低的原因，造成二氧化碳 的密度增加且速度降低，因而減少了阻力，因此通過第一個及最後一

個流道的二氧化碳流量大於相鄰流道。另外，由圖 1.11 可以發現當 二氧化碳壓力位於 8 MPa 時，溫度由 315 K 開始，當溫度提升時，黏 滯力隨溫度變化非常微小(幾乎不改變)，因此整體而言，二氧化碳的 流速分佈較水來得均勻。 圖 4.16 及圖 4.17 為水側及二氧化碳側熱傳量在各個流道的分佈 情形。在中間流道處，由於溫度及流速在中間流道並無太大差異，因 此，在中間流道的熱傳量分佈情形大致相同，而在水側第一個流道以 及最後一個流道處，由於接觸面僅有一面的影響，因此熱傳量皆偏低。 整體而言，在各個流道的熱傳量分佈處勢與速度分佈趨勢相近。

0 1 2 3 4 5 6 7 8 9 10 304 306 308 310 312 314 316 318 320 322 324 326 328 330 332 334 336 338 340 T(K) Number of channel T_H 2o,in =283K T_H 2O,out =328.8K 圖 4.12 水側出口溫度分佈 0 1 2 3 4 5 6 7 8 9 325 326 327 328 329 330 331 332 333 T(K) Number of channel T_CO 2in =370K T_CO 2out =329.9K 圖 4.13 二氧化碳側出口溫度分佈

0 1 2 3 4 5 6 7 8 9 10 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 υ_c=u_i/u_chm Number of channel H₂O u _chm=0.316 m/s u*_chm=0.320 m/s 圖 4.14 水側無因次化流速分佈 0 1 2 3 4 5 6 7 8 9 0.94 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 1.06 υ_c=u_i/u_chm Number of channel CO₂ u_chm=7.541 m/s u*_chm=5.662 m/s 圖 4.15 二氧化碳側無因次化流速分佈

圖 4.16 水側熱傳量與流道關係圖

4-4.2 改變冷媒流量 在本節中，將以固定水側流量，改變二氧化碳流量的方法，來探討流 量對熱傳量的影響。其中，由於水側流量固定，因此水側熱容量流率 固定，改變二氧化碳流量即為改變熱容比。 圖 4.18-圖 4.22 為在不同的 C 值的情況下( 2 / 2 CO H O C C )，水-二氧化 碳板式熱交換器的結果。其中， 2 CO C 為二氧化碳在入口處的熱容量流 率。 圖 4.18 為無因次化後水側速度分佈的分佈圖，由於水側質量流率 皆相同，因此圖 4.18 的速度分佈情形皆為溫度變化所造成，而整體 而言，不同的熱傳量對速度分佈的影響十分有限。 圖 4.19 為二氧化碳在不同 C 值時的速度分佈情形。由圖可以發現 二氧化碳的流動模式隨著 C 值的減少，而從遞減型流速分佈轉變為 遞增型，此現象可歸究於入出口的壓力以及速度變化所造成(如圖 4.20、圖 4.21)。在一般的流速分佈中(流體性質為常數者皆屬於這一 型)，流體在匯流區的阻力皆大於分流區，而造成流速分佈為遞減型。 而在二氧化碳流道中，在通道出口由於速度相對入口處較低，因而造 成匯流區的阻抗小於分流區，而造成流速分佈沿著流道編號而遞增的 情形。 圖 4.22 為不同 C 值熱傳結果，其中 Case1~5 分別為 C=0.1

、0.25、0.5、0.75、1。由圖可知總熱傳量隨著 C 值的增加而上升， 不過若將 C=0.1 與 C=0.5 相比，可以發現流量增加了 5 倍，而熱傳量 僅上增加了 1 倍，此現象是因為二氧化碳在 C=0.1 處出口溫度低於臨 界溫度，而二氧化碳在接近臨界點時比熱有一個相當大幅度的上升所 導致。

0 1 2 3 4 5 6 7 8 9 10 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 υ∗ c=ui/uchm Number of channel C=0.1 u*_chm=0.3165 C=0.25 u*_chm=0.3172 C=0.5 u*_chm=0.3184 C=0.75 u*_chm=0.3194 C=1 u*_chm=0.3203 u_chm=0.3156 圖 4.18 水側速度分佈圖 0 1 2 3 4 5 6 7 8 9 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 υ∗ c=ui/uchm Number of channel C=0.1 , u*_chm=0.1577 ,u_chm=0.7546 C=0.25 , u*_chm=0.7844 ,u_chm=1.886 C=0.5 , u*_chm=2.269 ,u_chm=3.772 C=0.75 , u*_chm=30891 ,u_chm=5.658 C=1 , u*_chm=5.662 ,u_chm=7.541 圖 4.19 二氧化碳側速度分佈圖

圖 4.20 C=0.1 時，分流區及匯流區壓力變化情形