基于Python語言輔助化工熱力學教學*

張建蘭，崔 健，蔡江濤，汪廣恒，張亞婷

(西安科技大學化學與化工學院，陜西 西安 710054)

化工熱力學是以化學熱力學和工程熱力學為基礎，在解決化工過程中實際問題而逐漸發展形成的一門實用工程學科，是化學工程的一個重要分支，也是化學工程與工藝專業的核心基礎課程[1]。化工熱力學的主要研究內容包括物系的熱力學性質計算和預測、相平衡和反應平衡計算、能量轉化和有效利用等實際問題，是化工過程研究、開發和設計的理論基礎，在培養學生工程觀點，提高學生工程能力方面起著舉足輕重的作用[2-3]。

化工熱力學前承學生們大一大二的高等數學、物理化學和計算機語言等課程，后啟大三大四的分離工程、反應工程和化工設計等專業課，是一門理論性、邏輯性非常強的學科。由于化工熱力學的研究對象是實際的化工過程，因此不能再沿用物理化學中的理想氣體或者理想溶液模型來進行工程計算，而需要更接近真實情況、準確性更好、更復雜的非理想模型，比如：描述真實流體p-V-T關系的RK，SRK、PR等立方型狀態方程，或者virial、BWR、MH等多參數狀態方程；而真實溶液則采用Margules、van Laar、Wilson等活度系數模型，或者是NRTL、UNIQUAC、UNIFAC等更復雜的活度系數模型[4-5]。這些復雜的模型使得化工熱力學這門課程涉及的公式多且抽象，計算量非常大，通常需要迭代試差，且采用手工計算較為困難，所以往往讓學生覺得枯燥無味，心生畏懼，從而限制了學生應用化工熱力學分析和解決實際問題的能力發展[6]，因此如何提高學生的學習興趣，使學生在掌握理論知識的同時，能夠理論聯系實際應用，培養學生解決復雜工程問題的計算能力，這是化工熱力學教學中的關鍵和難點。

隨著科學技術的高速發展以及計算機的廣泛應用，現代計算機技術為工程技術人員提供了十分豐富的程序語言和應用軟件腳本語言，如一般的高級程序語言(C/C++，VB，Java等)，也有集成強大數學運算能力的商業應用軟件(Matlab，Mathematica，Maple等)，以及專業的化工流程模擬軟件Aspen Plus，其中Matlab，Aspen Plus等軟件在化工熱力學的教學中得到了廣泛的應用[7-9]。這些軟件大大降低了模型的計算難度并有利于培養學生的應用能力，很大程度提高了學生的學習興趣和積極性。但是筆者在教學過程中發現，學生在使用這些軟件的過程中只是按照指令選取計算方法，對程序內部的計算過程和原理并不清楚，這顯然不利于培養學生以后解決未知復雜工程問題的能力。通過將計算原理進行語言編程，然后解決化工熱力學中的模型求解，不僅能將學生大一所學習的計算機知識進行應用，促進學科交叉融合，而且可以培養學生將理論知識應用于工程實際問題的能力，為培養優秀的工程技術人才奠定良好的基礎。近幾年來Python語言由于簡單易學、開源、可移植、可擴展、可嵌入以及面向對象等特點而被廣泛應用，因此筆者將結合化工熱力學的典型計算，重點介紹Python語言在流體p-V-T關系計算中的應用。

1 Python語言介紹

Python語言是由Guido van Rossum(吉多·范羅蘇姆)于1989年開發的一種腳本語言，是一種面向對象、解釋型計算機程序設計語言。第一個公開發行版本發布于1991年，經過三十多年的發展，Python語言憑借其語法精煉、可讀性強、兼容性和可移植性、強大而豐富的第三方庫等優點而成為主流的程序設計語言之一。Python語言是目前最接近自然語言的通用編程語言，憑借其強大的功能廣泛應用于GUI、Web框架、系統編程，網頁爬蟲，數據挖掘，科學計算在內幾乎所有的程序設計領域，是全球范圍內最大的單一語言編程社區。從解決計算問題角度，相比較C語言、VB語言和Java語言等高級語言，Python語言更適合以應用為主的非計算機專業學生，這是因為該語言只關心計算問題的求解，其輕量級的語法和高層次的語言表示表達了應用計算機解決問題的計算思維理念[10-12]。

2 Python語言應用實例

流體的p-V-T關系是化工熱力學的基石，是化工過程開發和設計、安全操作和科學研究必不可少的基礎數據，也是化工熱力學教學中的開篇之章，具有非常重要的地位。因此下面我們將重點介紹Python語言在流體p-V-T關系計算中的應用。

例：用RK方程計算異丁烷：(1)420 K，2.0 MPa時摩爾體積(實驗值為V=1 411.2 cm3/mol)；(2)在380 K時的飽和氣、液相摩爾體積，已知該溫度下的蒸汽壓是2.25 MPa(實驗值分別是866.1 cm3/mol和140.8 cm3/mol)。

這是陳新志版《化工熱力學》第二章中的一個實例[3]，在求解過程中課本直接給出了p-V相圖上380 K，420 K的兩條等溫線，用圖解法給出了方程的根。由于沒有具體的計算過程，學生們并不清楚p-V相圖和方程的根是如何得到的，因此工程計算能力并沒有得到鍛煉，無法解決其它狀況的熱力學性質計算，達不到培養學生們工程素養的目標。接下來我們看一下如何應用Python語言解決此類問題。

2.1 Python語言繪制p-V相圖

前面介紹到Python 語言含有豐富的第三方庫，因此我們可以直接使用Python語言中的matplotlib繪圖庫。matplotlib可以輕松地將數據圖形化，并且提供多樣化的輸出格式，很方便地將數據可視化并進行對比分析，在圖像的美化方面也比較完善。matplotlib中應用最廣的是matplotlib.pyplot模塊，下面是針對上述例題所編寫的程序：

>import matplotlib.pyplot as plt

>import numpy as np

>tc=408.1

>pc=3.648

>w=0.176

>T1=380

>T2=420

>a=0.42748*(8.314**2)*(tc**2.5)/pc

>b=0.08664*8.314*tc/pc

>x=np.arange(b+1，2000，0.01)

>y1=(8.314*T1)/(x-b)-(a/(T1**0.5))/(x*(x+b))

>y2=(8.314*T2)/(x-b)-(a/(T2**0.5))/(x*(x+b))

>plt.xlabel('v/(cm3/mol)')

>plt.ylabel('p/Mpa')

>plt.ylim(-1，10)

>plt.plot(x，y1，color='y'，label='380k，T

>plt.plot(x，y2，color='r'，label='420k，T>TC')

>plt.legend()

>plt.show()

通過運行就可以得到p-V相圖中T=380 K和T=420 K時的兩條等溫線。

圖1 基于Python語言繪制相關p-V相圖Fig.1 The related p-V diagram drawed by python language

通過繪制RK方程的相圖我們可以清晰的看出，p-V相圖中不同溫度下的等溫線形狀并不相同，因此不同壓力下求解摩爾體積的根的數目也就不同。通過繪制相圖，學生們可以很直觀的理解立方型狀態方程關于求解體積根的不同情況。

2.2 Python語言求解RK方程

2.2.1 直接迭代法求解RK方程

學生們在高等數學中已經學習過，在立方型或者高次型方程的求解過程中，采用解析法往往比較困難，更多的是采用迭代試差的方法進行求解。通常迭代的方法有兩種：直接迭代法和牛頓迭代法。通過將迭代原理利用Python語言進行編寫，就可以得到RK方程的解。并且我們可以很直觀的對比兩種不同迭代方法的不同。首先看一下直接迭代法的編寫程序：

>import sympy as sy

>T1=380

>p1=2.25

>a=0.42748*(8.314**2)*(408.1**2.5)/3.648

>b=0.08664*8.314*408.1/3.648

>x=sy.symbols('x')

>x0=8.314*T1/p1

>x1=b

>times1=0

>times2=0

>g=1e-10

>while True：

>temp_x=x0

>x0=8.314*T1/p1+b-(a/(T1**0.5))*(x0-b)/(p1*x0*(x0+b))

>times1+=1

>gap=abs(x0-temp_x)

>if gap

>print("可得飽和氣相的摩爾體積v(vapor)="，x0，"cm3/mol")

>print("v(vapor)迭代次數"，times1)

>break

>elif times1>1 000：

>print("1 000次迭代仍無法得出結果，請重新設置迭代初始值，或方程無解")

>break

>while True：

>temp_x=x1

>x1=(p1*(x1**3)-8.314*T1*(x1**2)-a*b/(T1**0.5))/(p1*(b**2)+ b*8.314*T1-a/(T1**0.5))

>times2+=1

>gap=abs(x1-temp_x)

>if gap

>print("可得飽和液相的摩爾體積v(liquid)="，x1，"cm3/mol")

>print("v(liquid)迭代次數"，times2)

>break

>elif times2>1000：

>print("1 000次迭代仍無法得出結果，請重新設置迭代初始值，或方程無解")

>break

通過運行結果我們可以得到T=380 K，p=2.25 MPa時的飽和氣、液相摩爾體積，同時還可以輸出迭代的次數，結果如下：

飽和氣相的摩爾體積v(vapor)=916.2 245 396 473 064 cm3/mol

v(vapor)迭代次數 43

飽和液相的摩爾體積v(liquid)=174.17 686 314 303 995 cm3/mol

v(liquid)迭代次數 111

同樣的思路，通過改變輸入溫度、壓力，就可以得到T=420 K，p=2.0 MPa時的摩爾氣相體積和迭代次數：

飽和氣相的摩爾體積v(vapor)=1404.5 052 868 409 211 cm3/mol

v(vapor)迭代次數 23

2.2.2 牛頓迭代法求解RK方程

根據牛頓迭代法原理，同樣利用Python語言進行編寫求解，下面是牛頓迭代法的編寫程序：

>import sympy as sy

>T1=380

>p1=2.25

>a=0.42748*(8.314**2)*(408.1**2.5)/3.648

>b=0.08664*8.314*408.1/3.648

>x=sy.symbols('x')

>x0=8.314*T1/p1

>x1=b+1

>def fun_value(_x)：

>return (8.314*T1)/(_x-b)-(a/(T1 ** 0.5))/(_x*(_x+b))-p1

>times1=0

>times2=0

>g=1e-10

>while True：

>temp_x=x0

>x0=x0-fun_value(x0)/sy.diff(fun_value(x)).subs(x，x0)

>times1+=1

>gap=abs(x0-temp_x)

>if gap

>print("可得飽和氣相的摩爾體積v(vapor)="，x0，"cm3/mol")

>print("v(vapor)迭代次數"，times1)

>break

>elif times1>1 000：

>print("1 000次迭代仍無法得出結果")

>break

>while True：

>temp_x=x1

>x1=x1-fun_value(x1)/sy.diff(fun_value(x)).subs(x，x1)

>times2+=1

>gap=abs(x1-temp_x)

> if gap

>print("可得飽和液相的摩爾體積v(liquid)="，x1，"cm3/mol")

>print("v(liquid)迭代次數"，times2)

>break

>elif times2>1 000：

>print("1 000次迭代仍無法得出結果")

>break

運行程序就可以得到T=380 K，p=2.25 MPa時的飽和氣、液相摩爾體積，結果如下：

飽和氣相的摩爾體積v(vapor)=916.224 539 647 234 cm3/mol

v(vapor)迭代次數 6

飽和液相的摩爾體積v(liquid)=174.176 863 143 347 cm3/mol

v(liquid)迭代次數 15

改變輸入溫度和壓強，同樣可以得到T=420 K，p=2.0 MPa時的摩爾氣相體積：

飽和氣相的摩爾體積v(vapor)=1404.50 528 684 091 cm3/mol

v(vapor)迭代次數 6

通過對比我們可以發現，兩種迭代方法求解的體積根的數值基本是一致的，但是牛頓迭代法的迭代次數卻是明顯少于直接迭代法的，這就是牛頓迭代法的優勢——收斂速度更快。但是牛頓迭代法的形式相比較直接迭代法而言更為復雜，因為要涉及方程的一階求導。通過原理的理解和編程過程的求解，學生們不僅可以利用兩種迭代方法解決物性的計算問題，還可以清晰的看到兩種迭代方法的不同和優缺點，極大的培養了學生們理論聯系實際，解決工程實際問題的計算能力。

3 結 語

通過Python語言在流體p-V-T關系計算中的應用，可以清晰地看到Python語言中強大的第三方庫和簡潔的編程語言能夠有效幫助學生解決化工熱力學中的復雜計算。通過將該編程運用于化工熱力學的教學中，不僅可以加深學生對于課程內容的理解和應用，提高課堂的教學效率，而且提高了學生的學習興趣和工程計算能力，為培養優秀的工程技術人才奠定良好的基礎。