符號運算軟體在化工熱力學應用之教材製作

行政院國家科學委員會專題研究計畫 成果報告

符號運算軟體在化工熱力學應用之教材製作

計畫類別： 個別型計畫 計畫編號： NSC94-2522-S-151-001- 執行期間： 94 年 08 月 01 日至 95 年 07 月 31 日 執行單位： 國立高雄應用科技大學化學工程與材料工程系 計畫主持人： 陳和瑞 報告類型： 精簡報告 處理方式： 本計畫涉及專利或其他智慧財產權，1 年後可公開查詢

中 華 民 國 95 年 9 月 18 日

符 號 運 算 軟 體 在 化 工 熱 力 學 應 用 之 教 材 製 作 計 劃 編 號 : NSC 94-2522-S-151-001 陳 和 瑞 國立高雄應用科技大學 化學工程與材料工程系

摘要

本計畫旨在利用電腦的運算能力，配合運用符號運算軟體(mathematica)程式的符號運算 能力，建構一套可提升科技大學化工系學生之熱力學性質計算能力與觀念建立的教材，內容 包含非理想純物質熱力學性質、多成分混合物熱力學性質、多成分混合物相平衡的計算。

簡介

化工熱力學為化工學門的基礎科目，其基礎理論在其他領域上的應用也非常的廣泛，例 如材料科學、儲能技術、電化學等。然而化工熱力學抽象的觀念對初學者在觀念上的建立比 較困難，而且對於真實條件的熱力學性質的計算也比較繁複，勢必要依賴電腦與適當的軟體 工具輔助工具始能竟其功。近年來微電子技術的進步，製造快速的微處理機，以及許多商業 上發行功能強大的運算軟體，例如 Mathematica、Mapple、Matlab、Mathcad 等。其中 Mathematica 符號運算軟體，它有符號輸入方式與繪圖功能，以及計算核心系統[1]，可讓繁複的計算步驟 程式碼符號化，因而可獲得簡化程式碼的攥寫，而且其又可很容易的輸出 postscript 的圖形， 所以有很好的視覺化輸出。因為 Mathematica 有這些優點，所以在學術上的研究常被作為計 算工具[2,3,4]，也普遍被作為教材的輔助工具，但在化學工程教材上的運用尚不多見。工欲善 其事必先利其器，本計畫使用微軟 OFFICE 系統及借助於 Mathematica 符號運算軟體的優點， 編制一套可提升科技大學化工系學生之熱力學性質計算能力與觀念建立的教材，其內容必須 涵蓋化工熱力學多成分相平衡計算的章節，以期建立學生良好的理論基礎。

教材內容

本計畫執行重點在於化工熱力教學，融入 mathematica 程式設計的教材製作，以期對學 者有較佳的學習效果，教材內容的設計分成六個章節，依次為： 1. 簡介 2. 能量平衡與熵平衡的應用 3. 真實物質熱力學性質 4. 純物質相平衡 5. 多成分混合物的熱力學性質 6. 混合物成分 Gibbs 自由能與逸壓估算 每一章節內容概述如下： 一、簡介

在此章節將描述化工熱力學常用的名詞定義，尤其對於可逆過程與不可逆過程的定義將 明確詳述。此外在下列章節中使用到的 mathematica 函數的語法在此作敘述並列舉使用範例。 二、能量平衡與熵平衡的應用 往 昔 常 見 化 工 熱 力 學 教 學 直 接 以 密 閉 系 統 ( U Q W  ) 及 機 械 能 方 程 式 ( 2 constant 2 u P gZ     )之公式作為計算問題的切入點，而少有強調公式使用之條件，遇到以 開放系統架構去解決的問題則不知如何著手。因此本章節一開始就引導建立開放系統能量均 衡式，即      





2 2 1 1 2 1 2 1 { ( / 2 )} ( / 2 ) ( ) {( ) / 2 } / 2 K K k k s k k k k K k k s k K k s k k d U M v dt dV M U v Q W P M PV dt dV M U PV v Q W P dt dV M H v Q W P dt                              









          並舉兩個典型的例題引導作題，這些例題在尚未介紹真實氣體之前均以理想氣體模型或以查 表方式來解題。例如： 【例題 1】 有一絕熱良好的箱子中間被隔膜分割成空間相 等的兩部分(參考右圖)，左邊內裝 1 kg 氣體右為真空，若將 隔膜移開，請計算最後的溫度與壓力。一開始(1) 氣體為溫 度 500 K；壓力 1 bar 的理想氣體。(2) 氣體為溫度 400℃； 壓力 60 bar 的蒸氣。 【例題 2】 如右圖之兩理想氣體進行混合，假設 混合槽為絕熱。最初槽壓為 0 bar，當槽溫達 195 ℃時 A 氣體之莫耳分率若干？ 此題目導入開放系統公式，並刪除不存在的項， 經推導之後的公式為



1





1



A VA A VB PA A A PB B A x C  x C T C T x C T x     再實際導入 mathematica 程式，使用內建的 Solve 函數，解出線性方程式的根，如下： 1 mol 氣體 真空 A 氣體、3 bar、25℃ CP= 20.9 J/mol·K B 氣體、3 bar、25℃ CP= 29.3 J/mol·K

對於熱力學初學者，熵是非常抽象的名詞。而熵為何要定義為S Q T  ？在本章節中所作的解 釋如下： 以卡諾引擎一個循環所得到 1 2 1 2 0 Q Q T  T 結果，假設一個可逆的循環過程可以由許多極微小的 卡諾循環來取代，則 1 0 M i i i Q T  



 dQ 0 T 





上式 dQ 0 T 





表示其積分與路徑無關，代表dQ T 是狀態函數，所以定義 dQ dS T  。 在熱力學系統，常常是能量均衡式與熵的均衡式必須合併，即 







1 2 2 1 { ( / 2 )} / 2 K k k gen k K k s k k dS Q M S S dt T d dV U M v M H v Q W P dt   dt           _{       }  





      上兩式忽略動能與未能項及W，共同銷去 Q 得_s





1 1 ˆ ˆ ˆ K K k k k k k k dU TdS dM H TS PdV TdS PdV dM G    



    



如以密閉系統的理想氣體可逆過程為例，消去不必要項之後可導出 * 2 2 * 2 * 2 2 2 1 1 1 1 1 1 ( , ) ( , ) _Pln T ln P _Pln V _Vln P S T P S T P C R C C T P V P           _{ } _{ } _{ } _{ }        關係式，過程若為絕熱可 逆則可得 * 2 2 2 1 1 2 1 1 ( , ) ( , ) 0 P R C P S T P S T P T T P      _{  }   上式的理想氣體應用例題如下： 【例題 3】 理想氣體在 1bar、290K 之條件進行絕熱壓縮至 10 bar，最後溫度幾度？ (C_V* 21 J/mol K) 程式碼：

而利用熵的均衡式檢驗程序若為可行，則必須滿足S_gen 0，本章列舉例子加以說明。 三、真實物質熱力學性質 在本章節中率先導出馬克思威爾關係式(Maxwell Relations)，並說明馬克思威爾關係式的用途 為：將方程式中出現之 T S V     _  換成 V P T     _  （或 T S P     _  換成 P V T     _  ），也就是將 S 轉換為 PVT 關係。 為使學者了解理想氣體與真實氣體的熱力學性質差異，在本章節中有詳細推導，綜合結果如 下表所示： 結果 項目 理想氣體 真實氣體模型 T P H V V T P T    _{ }  _  _  _     0 2 3  2 RT V a V b a P V V              

( van der Waal equation)

V T U P T P V T   _   _{ }  _  _{ }   0 2 RT a P

Vb V ( van der Waal equation)

P P T _P C V T V P T T         _{ }  _{ }     _  _  _{ }       0 ( 1 PV B RT  V ) V T T V C U V T V            _  _{ } 

     0 (Peng Robinson equation)

對於真實物質性質計算以 Peng Robinson 方程式為模型，藉由 mathematica 程式程式設計繪 出其相關性質與 PT 的關係(以 O2為例)如下：

本章沿用前第二章所舉之例題，但套入真實氣體狀態方程式並推導之，教材內詳述何將導出 之方程式攥寫成 mathematica 程式碼，所舉之相關例題如下： 【例題 4】 使用 Peng-Robinson 狀態方程式重複例題 2 的計算，其中 A 為 N2氣體；B 為 He 氣體。相關資料如下： N2：TC= 126.2 K；PC= 33.94 bar；0.04；CP 27.2 4.2 10 3T     He：TC= 5.19 K；PC= 2.27 bar；0.387；CP 27.2 4.2 10 3T     程式碼： 【例題 5】 使用 Peng-Robinson 狀態方程式，若為 N2氣體重複例題 3 的計算。 程式碼： -80 -60 -40 -20 S,J molK 0.1 0.5 1 5 10 50 100 P,bar T140 K 圖 2 Peng Robinson 方程式壓力對熵等溫線

四、純物質相平衡

本章首先推導相平衡與系統穩定條件的關係式，由於此兩條件使得以 Peng Robinson 方程式為 模型所繪出之 PVT 圖形，如圖 4 之 PVT 關係由原本的 12345 路徑曲線，改變為實際存在的 135 曲線。

根據平衡條件，汽、液之逸壓相等的關係可推導出 vap sat( vap liq)

liq PdV P V V



，此關係式相當 於圖 4 中 1231 所圍之面積等於 3453 所圍之面積，所以 0 ( ) sat PdV P V V   _   



_{。可利用此關係式} 計算出系統處於某一溫度下之飽和蒸氣壓，所舉的例子如下： 【例題 6】 用 Peng_Robinson 狀態方程式計算 O2在 140 K 之飽和蒸氣壓，TC= 154.6 K、PC= 50.46 bar、0.021。 1 2 3 4 5 6 7 sat P 圖 4 Peng Robinson 方程式模型體積等溫線

程式碼： In[3]:= 0.021; Tc154.6; Pc50.46; T140.; PRvpT,, Tc, Pc Out[4]= 28.0781 其中 PRvp[T,,Tc,Pc]為模組化後計算飽和蒸氣壓之函數，28.0781 bar 為 O2在 140 K 之飽和蒸 氣壓。 汽、液平衡曲線可藉由一系列的不同溫度(小於 Tc)計算出相對應的飽和蒸氣壓即可繪出，例如 圖 5 為 O2之 PV 汽液平衡曲線。其他的 PH、PS 之平衡曲線可經由算出之飽和蒸氣壓代入計算 式得之。 五、多成分混合物的熱力學性質 不同物質混合後因非理想條件故混合後熱力學性質產生改變，這種改變可從成分之部分莫耳 性質經由計算式估計之，主要項目為體積與焓的改變。本章是以硫酸與水的混合為例，資料 如下表，測量不同莫耳比混合溶液之比重，經由





_

_





_

_

2 ₂ 1 1 , 1 ₁ 2 2 _, mix mix T P mix mix T P V V x V V x V V x V V x        _         _  式可算出 各成分的部分莫耳體積。其中





1 _, mix T P V x   與





2 _, mix T P V x   可利用 mahtematica 內插函數計算， 0.05 0.1 0.5 1 5 V, m 3  kgmol 1 2 5 10 20 50 P,bar 圖 5 O2之 PV 汽液平衡曲線

並可繪出圖形如下：

1

x (H2O) (m3/mol)106 V_mix(m3/mol) 10 6 V_mix(m3/mol) 10 6

0.0 786.846 40.7221 0 0.1162 806.655 37.7015 -0.3883 0.2221 825.959 35.0219 -0.6688 0.2841 837.504 33.5007 -0.7855 0.3729 855.031 31.3572 -0.9174 0.4186 864.245 30.2812 -0.9581 0.5266 887.222 27.7895 -1.0032 0.6119 905.376 25.9108 -0.9496 0.7220 929.537 23.5759 -0.7904 0.8509 957.522 20.9986 -0.4476 0.9489 981.906 19.0772 -0.1490 1.0 997.047 18.0686 0 計算出的數據如下： x1 V  1V1 m3_{mol10}6 V—1 m3mol106 V  2V2 m3_{mol10}6 V—2 m3mol106 0. 3.57076 14.4976 1.010711017 40.7221 0.1162 3.07727 14.9911 0.0347172 40.6874 0.2221 2.37933 15.689 0.180477 40.5416 0.2841 1.98592 16.0824 0.309122 40.413 0.3729 1.62604 16.4423 0.495941 40.2261 0.4186 1.37903 16.6893 0.655065 40.067 0.5266 0.96665 17.1017 1.04388 39.6782 0.6119 0.552473 17.5159 1.57577 39.1463 0.722 0.243768 17.8246 2.20995 38.5121 0.8509 0.00049725 18.0689 3.00505 37.717 0.9489 0.0046855 18.073 3.00291 37.7192 1. 0. 18.0684 2.79986 37.9222 比照部分莫耳體積計算模式，硫酸－水系統之部分莫耳焓的數據如下：

六、混合物成分 Gibbs 自由能與逸壓估算 處於混合物平衡系統，各成分之逸壓相等(或逸壓係數)。由此條件可發展出對稱型與非對稱 型的汽液平衡關係式，即 L i i i V i _i L i i i i V i _i y K x y f K x _P                對稱型 非對稱型

本章中列舉混合物之起泡點(bubble point)之計算為例，非對稱型以 van Laar 及 Wilson 方程式 為例，對稱型之逸壓係數(、L_i  )以 Peng Robinson 方程式計算，其中起泡點計算範例如下：V_i 對稱型範例： 【例題 7】 己烷、庚烷、辛烷各佔 33%、37%、30%莫耳分率之液體，以對稱型計算在 1.2 bar 氣壓下此混合物之起泡點。相關資料如下表： kij _{C6H14 C7H16 C8H18 Tc, K Pc, bar}  C6H14 0 -0.008 0.003 507.4 29.69 0.296 C7H16 -0.008 0 0 540.2 27.36 0.351 C8H18 0.003 0 0 568.8 24.82 0.394 程式碼：

非對稱型範例：

【例題 8】 已知乙醇(1)－正己烷(2)在 331.25 K 時各純質成份的蒸氣壓為：

1 323.5

P mmHg P₂ 537.1 mmHg

van Laar 及 Wilson 係數為：2.409 1.970 及  ₁₂ 0.0952  ₂₁ 0.2713，利用非對 稱式計算在 331.25 K 及 x1= 0.412 時之泡點壓力與汽相組成。 假設總壓力不大 ( , , ) ( ) L sat sat i i i i i f T P x x P T  及液壓係數定義 v v i i i f y P    _v v i i i f y P 汽液平衡： V L i i f f  v_iy P_i _ix P_{i i}sat( )T _isat 因壓力不大， v sat i i   所以 ( ) sat i i i i y Px P T  P 



_ix P_{i i}sat( )T

a. van Laar equation

從 _{1 2} 1 2 ln 1 x x             、 _{2 2} 2 1 ln 1 x x             式

1 2 2 1 2 2.409 ln 0.699 2.409 0.411 1 1 1.970 0.588 x x          _     _ _   1 2.011  2 2 2 1 1.970 ln 0.4195 1.970 0.588 1 1 x _{2.409 0.411} x            _      _{ }   2 1.521   2.011 0.412 323.5 1.521 0.588 537.1 748.3 P        mmHg 1 2.011 0.412 323.5 0.358 748.3 y     、y₂ 0.642 b. Wilson equation 從 12 1 1 12 2 2 1 12 2 2 21 1 21 12 2 2 21 1 1 2 21 1 1 12 2 21 ln ln( ) ln ln( ) x x x x x x x x x x x x x x          _  _            _  _     兩式 12 1 1 12 2 2 1 12 2 2 21 1 21 ln ln(x x ) x 0.6509 x x x x    _    _     1 1.917 21 12 2 2 21 1 1 2 21 1 1 12 2 ln ln(x x ) x 0.4329 x x x x    _    _      2 1.542 742.4 P  mmHg、y₁0.344 上面所提之 van Laar 方程式之係數可從實驗數據以非線性回歸方式計算出，所舉例子如下： 【例題 8】苯與 2,2,4 三甲基戊烷在 55℃的汽液平衡數據如下表[5]，以非線性回歸計算出 van Laar 方程式中之、 之參數。苯與 2,2,4 三甲基戊烷(TMP)在 55℃的飽和蒸氣壓分別為 0.43596 及 0.23738 bar。 xB yB P(Bar) ( 1) ( 1 / 1 1 ) vap P x P y Ln Ln   ( 2) ( 2 / 2 2 ) vap P x P y Ln Ln   0.0819 0.1869 0.26892 0.341938 0.0032995 0.2192 0.4065 0.31573 0.294936 0.0109434 0.3584 0.5509 0.35463 0.223428 0.0446933 0.3831 0.5748 0.36088 0.216721 0.046736 0.5256 0.6786 0.39105 0.146776 0.109809 0.8478 0.8741 0.43277 0.0232059 0.410837 0.9872 0.9863 0.43641 0.0001196 0.676871

假設氣相為理想氣體且忽略 Poynting corrections，所以 P y y P T f i V i ( , , ) vap i i i i i vap i i i P yP yP xP x   / 以上式計算上表第四與第五行數據，又 van Laar 方程式： ₂ 2 1 1 1 ln         x x     、 ₂ 1 2 2 1 ln         x x     利用上表第一行與第四行數據進行非線性回歸，結果與程式碼如下：

結論

圖表 6 苯與 TMP 烷溶液以 van Laar 方程式非線性回歸 之活性係數與組成關係曲線 B  TMP 

Mathematica 符號運算軟體，它有符號輸入方式與強勢的繪圖功能與內鍵函數庫，可讓 繁複的計算步驟程式碼符號化，因而可獲得簡化程式碼的攥寫，因此非常適合將其運用於化 工熱力學的教學。本計劃製作一套可提升科技大學化工系學生之化工熱力學基礎與觀念建立 的教材，內容包含非理想純物質熱力學性質、多成分混合物熱力學性質、多成分混合物相平 衡的計算。

參考文獻

1. Wolfram, S., The Mathematica Book 4rd Ed., New Yark: Cambridge University Press(1999)

2. Kaziro, R. W., “Mathematica-assisted learning in physical chemistry”, Journal of Chemical Information and Computer Sciences 391 Jan-Feb 1999 p96-103.

3. Mordechai S., Neima B. and Michael P., “Compare software for interactive solution of systems of nonlinear algebraic equations”, Computer chem. Engng Vol. 22, No. 1-2, pp. 323-331, 1998

4. 郭東義, 陳和瑞, 利用符號運算設計最適 PID 控制器參數, Proceeding of the 1999 CIChE Annual Meeting and Conferences Chemical Engineering Technology.

5. Stadley I. S., Chemical and engineering thermodynamics, 3rd Ed., John Wiley and Sons(1997)