Matlab在二元完全理想系相平衡關系計算中的應用

黃英虎，農正東，邢靚，郭效瑛，李素萍，馬理，李偉

（百色學院化學與環境工程學院，廣西 百色 533000）

精餾作為一種典型的分離單元操作，其設計計算在化工計算中尤為重要。在精餾塔的逐板計算中，相平衡關系的計算因涉及非線性方程（組）的求解，計算繁瑣，是全過程中最難穿越的“荊棘”，導致本就復雜的計算過程愈加困難[1-2]。因此，精餾計算常需依托軟件編程來完成。在眾多的編程軟件中，Matlab因其數值計算功能卓異、庫函數浩博、語言簡易、語法靈活、易于可視化等優勢而脫穎而出，為化工計算提供了極大的便利[3-5]。

在所有的氣液兩相體系中，二元完全理想系的相平衡計算最為簡單。在完全理想系的相平衡關系中，相對揮發度是關鍵參數。在恒壓精餾過程中，因其是溫度的函數，所以其值會隨塔板溫度的變化而變化[6]。因此，確定合適的相對揮發度對精餾計算非常重要。若有相平衡實驗數據，可以擬合出相對揮發度[7]；反之，則需要選取合適的平均相對揮發度[8]。本文通過幾何平均法、對數平均法和算術平均法求得平均相對揮發度，并將三者用于給定組成的泡露點溫度、給定溫度的飽和氣液相組成的計算。最后以直接計算法的結果作為基準，通過相對偏差分析對三種方法進行評價。本文以Matlab 函數或腳本實現全部計算過程。

1 計算原理

1.1 飽和蒸氣壓的計算

組分m的飽和蒸氣壓可由Antoine方程[9]計算：

式中Am、Bm和Cm為組分m 的Antoine 參數；(T)為組分m于溫度T（K）下的飽和蒸氣壓（mmHg）。

1.2 相對揮發度

當完全理想系達到相平衡時，對組分m，有：

式中，ym、xm分別為溫度T（K）下組分m于飽和氣液相中的摩爾分數；p（mmHg）為體系的總壓，常壓時，p=760 mmHg。

1.2.1 塔頂和再沸溫度及另一相組成

本節采用直接計算法計算塔頂和再沸器的操作溫度及未知的另一相組成。

在溫度T（K）下，將式（2）應用于輕、重兩組分，建立一側歸零化的方程組：

式中，y1、x1分別為飽和氣液相中輕組分的摩爾分數，輕重兩組分的飽和蒸氣壓與溫度T的關系滿足式（1）。為使敘述簡練，下文將輕組分的摩爾分數簡稱為組成。

式中，xD和xW分別為塔頂和塔底產品組成，xD為給定參數，xW可由全塔物料衡算求解；TD（K）為塔頂的操作溫度；xLD為自塔頂第一塊理論板下降的飽和液相組成；TW（K）為再沸器的操作溫度；yW為再沸器產生的飽和氣相組成。

1.2.2 塔頂和塔底的相對揮發度

塔頂和再沸器操作溫度下的相對揮發度可用式（4）計算：

式中，k1(T)，k2(T)分別為溫度T（K）下輕組分和重組分的氣液平衡常數；y(T)，x(T)分別為該溫度下的飽和氣液相組成。

取T=TD和T=TW，則可由1.2.1提及或計算出的xD[y（TD）]、xLD[x（TD）]、yW[y（TW）]、xW[x（TW）]對應計算出α(TD)和α(TW)，兩者分別為塔頂和再沸器操作溫度下的相對揮發度。

1.2.3 平均相對揮發度

幾何平均法采用式（5）計算平均相對揮發度：

對數平均法采用式（6）計算平均相對揮發度αˉlog：

算術平均法采用式（7）計算平均相對揮發度αˉari：

三種方法因都采用了相對揮發度，所以統稱為相對揮發度法。在本文案例條件下，三種相對揮發度的最終計算結果分別為2.469、2.469和2.471。

1.3 相平衡方程式

式中，x為飽和液相組成；α為某一種平均相對揮發度，可由式（5）（6）或（7）求得；ye為與x同溫的飽和氣相組成。

1.4 泡露點溫度的計算

給定組成（本文選進料組成xF）的泡露點溫度可以用直接計算法和相對揮發度法來求取。

1.4.1 直接計算法

將式（9）與式（1）聯立即可求出泡點溫度Tb（K）

將T=Td和y1=xF代入式（3），并聯立式（1）即可求出Td（K）。

1.4.2 相對揮發度法

相對揮發度法由式（10）計算Tb（K）：

相對揮發度法由式（11）列出的方程組計算Td（K）：

1.5 飽和氣液相組成的計算

給定溫度（本文選擇進料溫度TF，要求Tb≤TF≤Td）下的飽和氣液相組成亦可由上述兩類方法計算。

1.5.1 直接計算法

進料溫度下的飽和液氣相組成xe、ye可分別由式（12）和（13）計算：

1.5.2 相對揮發度法

相對揮發度法由式（14）列出的方程組求xe、ye

式中，(TF)為進料溫度下輕組分的飽和蒸氣壓（mmHg），可由式（1）計算。

2 計算示例與程序

計算案例：已知苯-甲苯精餾塔的原料流量為10 kmol/h，其中苯含量為50%（摩爾分數，下同），在1 atm下進行精餾，塔頂餾出液的流量為5 kmol/h，苯含量為95%。在操作條件下，苯-甲苯可按完全理想體系進行計算。苯（1）、甲苯（2）的飽和蒸氣壓可由式（1）計算，其對應的Antonie參數見表1。

表1 苯和甲苯的Antonie參數[9]

試求：①分別以三種方法計算進料組成的泡露點溫度tb和td；②在兩溫度間選取進料溫度tF，分別以三種方法計算進料溫度下的飽和氣液相組成。

2.1 塔頂、再沸器的操作溫度及另一相組成的計算

新建一個名為IDVLEpy的函數，該函數可用于以直接計算法計算塔頂液相組成和塔頂的操作溫度，編輯代碼如下：

function xt= IDVLEpy（xx）

%IDVLEpy 函數在已知壓力P，塔頂氣相組成xD下，返回關于飽和液相組成和塔頂操作溫度的方程組

%xx 的兩個元素依次為飽和液相組成xLD 和塔頂的操作溫度tD（K）

global CA P xD

%以下兩步以Antoine 方程計算塔頂溫度下輕重兩組分的飽和蒸氣壓列向量Ps（mmHg）

lnPs=CA（:，1）-CA（:，2）./（CA（:，3）+xx（2））；

Ps=exp（lnPs）；

%以下兩步分別給出塔頂飽和氣液相中輕重兩組分的摩爾分數列向量y和x

y=[xD；1-xD]；x=[xx（1）；1-xx（1）]；

xt=y*P-x.*Ps；%建立求解方程組，取一側化為0后的表達式

end

新建一個名為IDVLEpx的函數，該函數可用于以直接計算法計算再沸器氣相組成和再沸器的操作溫度，編輯代碼如下：

function yt = IDVLEpx（xx）

%IDVLEpx 函數在已知壓力P，塔底再沸器液相組成xW下，返回關于飽和氣相組成和再沸器的操作溫度的方程組

%xx 的兩個元素依次為飽和氣相組成和塔底再沸器的操作溫度，K

global CA P xW

%以下兩步以Antoine方程計算塔底再沸器溫度xx（2）（K）下輕重兩組分的飽和蒸氣壓列向量，mmHg

lnPs=CA（:，1）-CA（:，2）./（CA（:，3）+xx（2））；

Ps=exp（lnPs）；

%以下兩步分別給出再沸器飽和氣液相中輕重兩組分的摩爾分數列向量y和x

y=[xx（1）；1-xx（1）]；x=[xW；1-xW]；

yt=y*P-x.*Ps；%建立求解方程組，取一側化為0后的表達式

end

2.2 三種平均相對揮發度的計算

新建一個名為rvmean 的函數，可以選擇幾何平均法、對數平均法或算術平均法之一計算全塔輕組分對重組分的平均相對揮發度alpha，編輯代碼如下：

function alpha= rvmean（xD，xW，xLD，yW，str）

% rvmean 函數以給定的方法計算全塔輕組分對重組分的平均相對揮發度alpha

%xD和xW分別為塔頂產品和塔底產品組成

%xLD 為塔頂溫度下飽和液相中輕重組分的摩爾分數列向量；yW 為再沸溫度下飽和氣相中輕重組分的摩爾分數列向量

%str為方法判斷字符串

%以下兩步分別計算塔頂及再沸器操作溫度下的氣液平衡常數kD和kW

kD=[xD；1-xD]./xLD；kW=yW./[xW；1-xW]；

%以下兩步分別計算塔頂及再沸器中輕組分對重組分的相對揮發度

alpha_D=kD（1）/kD（2）；alpha_W=kW（1）/kW（2）；

if strcmp（str，′geo′） %幾何平均值法

alpha=sqrt（alpha_D*alpha_W）；

elseif strcmp（str，′log′）%對數平均值法

alpha=（alpha_D-alpha_W）/log（alpha_D/alpha_W）；

elseif strcmp（str，′ari′） %算術平均值法

alpha=（alpha_D+alpha_W）/2；

end

end

2.3 泡露點溫度的計算

新建一個函數命名為BPT，可為直接計算法求給定組成的泡點溫度提供方程，編輯代碼如下：

function fun = BPT（t）

%BPT函數返回直接計算法求泡點溫度的方程

% t為待求解的給定組成的泡點溫度，K

global CA P xF

%將Antoine 參數矩陣CA、總壓P（mmHg），給定組成xF定義為全局變量

Ps=exp（CA（:，1）-CA（:，2）.（/CA（:，3）+t））；%以Antoine公式計算輕重兩組分在（tK）下的飽和蒸氣壓列向量，mmHg

x=[xF；1-xF]；%給出飽和液相中輕重兩組分的摩爾分數列向量

fun=sum（x.*Ps）-P；%建立泡點溫度求解方程，取一側化為0后的表達式

end

新建一個函數命名為DPT，可為直接計算法求解給定組成的露點溫度提供方程組，編輯代碼如下：

function fun = DPT（xx）

%DPT 函數返回直接計算法求解露點的溫度和露點液相組成的方程組

%解向量xx 的兩個元素依次為露點溫度下的飽和液相組成及給定組成的露點溫度，K

global CA P xF

y=[xF；1-xF]；%給出飽和氣相中輕重兩組分的摩爾分數列向量

Ps=exp（CA（:，1）-CA（:，2）./（CA（:，3）+xx（2）））；%計算輕重兩組分在溫度xx（2）（K）下的飽和蒸氣壓列向量，mmHg

x=[xx（1）；1-xx（1）]；%給出露點溫度下飽和液相中輕重兩組分的摩爾分數列向量

fun=x.*Ps-y*P；%建立露點溫度求解方程組，取一側化為0后的表達式

end

編輯一個名為Tdx_rv的函數，為相對揮發度法計算露點溫度提供方程組

function f= Tdx_rv（xx）

% Tdx_rv 函數返回相對揮發度法求露點溫度和露點液相組成的方程組

% 解向量xx 的兩個元素依次是露點液相組成xd和露點溫度Td（K）

global CA P alpha xF

f（1）=xx（1）*exp（CA（1，1）-CA（1，2）/（CA（1，3）+xx（2）））-xF*P；

f（2）= alpha*xx（1）/（1+（alpha-1）*xx（1））-xF；

end

編輯一個名為BDPT_2M 的函數，可以給定的方法和參數計算泡露點溫度

function [tb，td] = BDPT_2M（xF，alpha，str）

%BPT_2M函數可以用兩類方法計算泡點溫度，tb、td分別為待求解的給定組成的泡點及露點溫度

global CA P

if strcmp（str，′dc′）

%以下兩步用以求解給定組成的泡點溫度

t0=95+273.15；%給出初值t0（K）

tb=fzero（@BPT，t0）-273.15；%調用fzero 函數求解tb（℃）

%以下三步用以求解給定組成的露點溫度td（℃）

xx0=[xF，t0]；%給出初值向量xx0

r=fsolve（@DPT，xx0）；%調用fsolve 函數解DPT 函數所定義的方程組

td=r（2）-273.15；

elseif strcmp（str，′rv′）

y=@（x）alpha*x./[1+（alpha-1）*x]；

f=@（T）xF*exp［CA（1，1）-CA（1，2）/（CA（1，3）+T）］-y（xF）*P；

tb=fzero（f，373）-273.15；

xx0=[xF，373]；

r=fsolve（@Tdx_rv，xx0）；

td=r（2）-273.15；

end

end

2.4 飽和氣液相組成的計算

編輯一個名為sxy 的函數，用給定的方法和參數計算飽和氣液相組成。

function [ x，y ] = sxy（TF，xF，alpha，str）

%sxy可以用直接計算法或相對揮發度法計算給定溫度下的飽和液、氣相組成x、y

% TF（K）為給定溫度，xF為給定組成，alpha為某一種相對揮發度，str為方法字符串

global CA P

if strcmp（str，′dc′）%直接計算法

Psf=exp（CA（:，1）-CA（:，2）./（CA（:，3）+TF））；%計算出給定溫度下輕重兩組分的飽和蒸氣壓列向量Psf（mmHg）

fun=@（x）Psf（1）*x+Psf（2）*（1-x）-P；

x=fzero（fun，xF）；%求給定溫度下的飽和液相組成x

y=x*Psf（1）/P；%求給定溫度下的飽和氣相組成y elseif strcmp（str，′rv′）

y=@（x）alpha*x./（1+（alpha-1）*x）；

f=@（x）x*exp（CA（1，1）-CA（1，2）/（CA（1，3）+TF））-y（x）*P；

x=fzero（f，xF）；

y=y（x）；

end

end

2.5 主程序

新建一個名為tbtdxeye 的腳本文件，編輯代碼如下：

clc，clear

global CA P xF xD xW alpha

%將Antoine 參數矩陣CA、總壓P（mmHg）、進料液組成xF、塔頂產品組成xD 和塔釜產品組成xW、所選相對揮發度alpha定義為全局變量

CA=[15.9008 2788.51-52.36；16.0137 3096.52-53.67]；

%定義Antoine 參數矩陣CA：第一行元素依次為苯的Antoine 參數A，B，C；第二行元素依次為甲苯的Antoine參數A，B，C

P=760；%給定氣相總壓（mmHg）

p=input（′請輸入原料液和塔頂產品摩爾流量（kmol/h）、原料液和塔頂組成：′）；%本例輸入[10 5 0.5 0.95]

xF=p（3）；xD=p（4）；

xW=（p（1）*xF-p（2）*xD）/（p（1）- p（2））；%根據物料衡算式求塔釜產品組成xW

%以下數步用以計算xLD

xt0=[xD 81+273.15]；%定義初值向量xt0

xt=fsolve（@IDVLEpy，xt0）；%調用fsolve 函數求解IDVLEpy定義的方程組

xLD=xt（1）；xLD=[xLD；1-xLD]；

%以下數步用以計算yW

yt0=[xW 110+273.15]；%給出解IDVLEpx函數的初始值向量yt0

yt=fsolve（@IDVLEpx，yt0）；%調用fsolve函數解IDVLEpx定義的方程組

yW=yt（1）；yW=[yW；1-yW]；

str={′geo′，′log′，′ari′}；

a=zeros（size（str））；

for i=1∶length（str）

a（i）=rvmean（xD，xW，xLD，yW，str{i}）；

end

tb= zeros（4，1）；td=tb；

[tb（1），td（1）]= BDPT_2M（xF，a（1），′dc′）；

for i=1∶length（a）

alpha=a（i）；

[tb（i+1），td（i+1）]= BDPT_2M（xF，alpha，′rv′）；

end

tF=input（′請輸入進料溫度（介于泡露點之間）：（℃）′）；

TF=tF+273.15；%進料溫度單位換算（K）

x= zeros（4，1）；y=x；

[ x（1），y（1）] = sxy（TF，xF，a（1），′dc′）；

for i=1∶length（a）

alpha=a（i）；

[x（i+1），y（i+1）]= sxy（TF，xF，alpha，′rv′）；

end

3 結果與討論

運行tbtdxeye腳本，命令窗口顯示：

請輸入原料液和塔頂產品摩爾流量（kmol/h）、原料液和塔頂組成：

輸入[10 5 0.5 0.95 ]并回車，命令窗口顯示：

請輸入進料溫度（介于泡露點之間）：（℃）

輸入95 并回車，計算完成后，在Workspace 中出現全部計算數據。

從Workspace 中復制相關計算結果，可得表2。從表2 可以看出，當以直接計算法為基準時，幾何平均法計算的tb、td、xe和ye整體偏差最?。ㄏ鄬ζ钜来螢?.10%、0.34%、0.03%和0.03%）；對數平均法的整體偏差也較?。ㄏ鄬ζ钜来螢?.10%、0.35%、0.06%和0.06%）；算術平均法的偏差整體相對較大（相對偏差依次為0.09%、0.36%、0.12%和0.12%）。

表2 四種計算方法的計算結果

由此可見，與直接計算法結果相比，幾何平均法的整體計算效果相對最好，對數平均法次之，算術平均法最差；無論相對揮發度采用哪一種方法計算，三種相對揮發度法的各個計算結果的相對偏差均在0.2%以內，因此三種方法均可用于二元完全理想系相平衡關系的計算。

4 結論

本文創建了八個函數用于計算氣液相平衡關系：IDVLEpy、IDVLEpx、rvmean、BPT、DPT、Tdx_rv、BDPT_2M、sxy。其中，IDVLEpy、IDVLEpx、BPT、DPT、Tdx_rv返回的均是需要求解的方程或方程組：IDVLEpy、IDVLEpx專用于以直接計算法計算塔頂和塔底再沸器的飽和液相或氣相組成及其操作溫度，而BPT、DPT、Tdx_rv 均為BDPT_2M 返回待求解方程組；rvmean 函數可根據給定的方法計算三種相對揮發度；BDPT_2M、sxy 可以兩類方法（共四種方法）分別計算給定組成的泡露點溫度、給定溫度下的飽和氣液相組成。以苯-甲苯體系為例，通過主程序tbtdxeye 調用相關函數可實現以四種方法計算泡露點溫度及飽和氣液相組成。通過相對偏差分析發現，以直接計算法為基準，三種相對揮發度法的整體精度順序為：幾何平均法>對數平均法>算術平均法。三種方法均可計算二元完全理想系相平衡關系，可進一步用于后續的精餾計算中。