黃英虎，農(nóng)正東，邢靚，郭效瑛，李素萍，馬理，李偉

（百色學(xué)院化學(xué)與環(huán)境工程學(xué)院，廣西 百色 533000）

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

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

1 計(jì)算原理

1.1 飽和蒸氣壓的計(jì)算

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

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

1.2 相對揮發(fā)度

當(dāng)完全理想系達(dá)到相平衡時(shí)，對組分m，有：

式中，ym、xm分別為溫度T（K）下組分m于飽和氣液相中的摩爾分?jǐn)?shù)；p（mmHg）為體系的總壓，常壓時(shí)，p=760 mmHg。

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

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

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

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

式中，xD和xW分別為塔頂和塔底產(chǎn)品組成，xD為給定參數(shù)，xW可由全塔物料衡算求解；TD（K）為塔頂?shù)牟僮鳒囟?；xLD為自塔頂?shù)谝粔K理論板下降的飽和液相組成；TW（K）為再沸器的操作溫度；yW為再沸器產(chǎn)生的飽和氣相組成。

1.2.2 塔頂和塔底的相對揮發(fā)度

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

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

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

1.2.3 平均相對揮發(fā)度

幾何平均法采用式（5）計(jì)算平均相對揮發(fā)度：

對數(shù)平均法采用式（6）計(jì)算平均相對揮發(fā)度αˉlog：

算術(shù)平均法采用式（7）計(jì)算平均相對揮發(fā)度αˉari：

三種方法因都采用了相對揮發(fā)度，所以統(tǒng)稱為相對揮發(fā)度法。在本文案例條件下，三種相對揮發(fā)度的最終計(jì)算結(jié)果分別為2.469、2.469和2.471。

1.3 相平衡方程式

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

1.4 泡露點(diǎn)溫度的計(jì)算

給定組成（本文選進(jìn)料組成xF）的泡露點(diǎn)溫度可以用直接計(jì)算法和相對揮發(fā)度法來求取。

1.4.1 直接計(jì)算法

將式（9）與式（1）聯(lián)立即可求出泡點(diǎn)溫度Tb（K）

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

1.4.2 相對揮發(fā)度法

相對揮發(fā)度法由式（10）計(jì)算Tb（K）：

相對揮發(fā)度法由式（11）列出的方程組計(jì)算Td（K）：

1.5 飽和氣液相組成的計(jì)算

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

1.5.1 直接計(jì)算法

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

1.5.2 相對揮發(fā)度法

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

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

2 計(jì)算示例與程序

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

表1 苯和甲苯的Antonie參數(shù)[9]

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

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

新建一個(gè)名為IDVLEpy的函數(shù)，該函數(shù)可用于以直接計(jì)算法計(jì)算塔頂液相組成和塔頂?shù)牟僮鳒囟?，編輯代碼如下：

function xt= IDVLEpy（xx）

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

%xx 的兩個(gè)元素依次為飽和液相組成xLD 和塔頂?shù)牟僮鳒囟萾D（K）

global CA P xD

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

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

Ps=exp（lnPs）；

%以下兩步分別給出塔頂飽和氣液相中輕重兩組分的摩爾分?jǐn)?shù)列向量y和x

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

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

end

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

function yt = IDVLEpx（xx）

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

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

global CA P xW

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

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

Ps=exp（lnPs）；

%以下兩步分別給出再沸器飽和氣液相中輕重兩組分的摩爾分?jǐn)?shù)列向量y和x

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

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

end

2.2 三種平均相對揮發(fā)度的計(jì)算

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

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

% rvmean 函數(shù)以給定的方法計(jì)算全塔輕組分對重組分的平均相對揮發(fā)度alpha

%xD和xW分別為塔頂產(chǎn)品和塔底產(chǎn)品組成

%xLD 為塔頂溫度下飽和液相中輕重組分的摩爾分?jǐn)?shù)列向量；yW 為再沸溫度下飽和氣相中輕重組分的摩爾分?jǐn)?shù)列向量

%str為方法判斷字符串

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

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

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

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′）%對數(shù)平均值法

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

elseif strcmp（str，′ari′） %算術(shù)平均值法

alpha=（alpha_D+alpha_W）/2；

end

end

2.3 泡露點(diǎn)溫度的計(jì)算

新建一個(gè)函數(shù)命名為BPT，可為直接計(jì)算法求給定組成的泡點(diǎn)溫度提供方程，編輯代碼如下：

function fun = BPT（t）

%BPT函數(shù)返回直接計(jì)算法求泡點(diǎn)溫度的方程

% t為待求解的給定組成的泡點(diǎn)溫度，K

global CA P xF

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

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

x=[xF；1-xF]；%給出飽和液相中輕重兩組分的摩爾分?jǐn)?shù)列向量

fun=sum（x.*Ps）-P；%建立泡點(diǎn)溫度求解方程，取一側(cè)化為0后的表達(dá)式

end

新建一個(gè)函數(shù)命名為DPT，可為直接計(jì)算法求解給定組成的露點(diǎn)溫度提供方程組，編輯代碼如下：

function fun = DPT（xx）

%DPT 函數(shù)返回直接計(jì)算法求解露點(diǎn)的溫度和露點(diǎn)液相組成的方程組

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

global CA P xF

y=[xF；1-xF]；%給出飽和氣相中輕重兩組分的摩爾分?jǐn)?shù)列向量

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

x=[xx（1）；1-xx（1）]；%給出露點(diǎn)溫度下飽和液相中輕重兩組分的摩爾分?jǐn)?shù)列向量

fun=x.*Ps-y*P；%建立露點(diǎn)溫度求解方程組，取一側(cè)化為0后的表達(dá)式

end

編輯一個(gè)名為Tdx_rv的函數(shù)，為相對揮發(fā)度法計(jì)算露點(diǎn)溫度提供方程組

function f= Tdx_rv（xx）

% Tdx_rv 函數(shù)返回相對揮發(fā)度法求露點(diǎn)溫度和露點(diǎn)液相組成的方程組

% 解向量xx 的兩個(gè)元素依次是露點(diǎn)液相組成xd和露點(diǎn)溫度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

編輯一個(gè)名為BDPT_2M 的函數(shù)，可以給定的方法和參數(shù)計(jì)算泡露點(diǎn)溫度

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

%BPT_2M函數(shù)可以用兩類方法計(jì)算泡點(diǎn)溫度，tb、td分別為待求解的給定組成的泡點(diǎn)及露點(diǎn)溫度

global CA P

if strcmp（str，′dc′）

%以下兩步用以求解給定組成的泡點(diǎn)溫度

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

tb=fzero（@BPT，t0）-273.15；%調(diào)用fzero 函數(shù)求解tb（℃）

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

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

r=fsolve（@DPT，xx0）；%調(diào)用fsolve 函數(shù)解DPT 函數(shù)所定義的方程組

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 飽和氣液相組成的計(jì)算

編輯一個(gè)名為sxy 的函數(shù)，用給定的方法和參數(shù)計(jì)算飽和氣液相組成。

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

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

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

global CA P

if strcmp（str，′dc′）%直接計(jì)算法

Psf=exp（CA（:，1）-CA（:，2）./（CA（:，3）+TF））；%計(jì)算出給定溫度下輕重兩組分的飽和蒸氣壓列向量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 主程序

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

clc，clear

global CA P xF xD xW alpha

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

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

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

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

p=input（′請輸入原料液和塔頂產(chǎn)品摩爾流量（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））；%根據(jù)物料衡算式求塔釜產(chǎn)品組成xW

%以下數(shù)步用以計(jì)算xLD

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

xt=fsolve（@IDVLEpy，xt0）；%調(diào)用fsolve 函數(shù)求解IDVLEpy定義的方程組

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

%以下數(shù)步用以計(jì)算yW

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

yt=fsolve（@IDVLEpx，yt0）；%調(diào)用fsolve函數(shù)解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（′請輸入進(jìn)料溫度（介于泡露點(diǎn)之間）：（℃）′）；

TF=tF+273.15；%進(jìn)料溫度單位換算（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 結(jié)果與討論

運(yùn)行tbtdxeye腳本，命令窗口顯示：

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

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

請輸入進(jìn)料溫度（介于泡露點(diǎn)之間）：（℃）

輸入95 并回車，計(jì)算完成后，在Workspace 中出現(xiàn)全部計(jì)算數(shù)據(jù)。

從Workspace 中復(fù)制相關(guān)計(jì)算結(jié)果，可得表2。從表2 可以看出，當(dāng)以直接計(jì)算法為基準(zhǔn)時(shí)，幾何平均法計(jì)算的tb、td、xe和ye整體偏差最?。ㄏ鄬ζ钜来螢?.10%、0.34%、0.03%和0.03%）；對數(shù)平均法的整體偏差也較小（相對偏差依次為0.10%、0.35%、0.06%和0.06%）；算術(shù)平均法的偏差整體相對較大（相對偏差依次為0.09%、0.36%、0.12%和0.12%）。

表2 四種計(jì)算方法的計(jì)算結(jié)果

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

4 結(jié)論

本文創(chuàng)建了八個(gè)函數(shù)用于計(jì)算氣液相平衡關(guān)系：IDVLEpy、IDVLEpx、rvmean、BPT、DPT、Tdx_rv、BDPT_2M、sxy。其中，IDVLEpy、IDVLEpx、BPT、DPT、Tdx_rv返回的均是需要求解的方程或方程組：IDVLEpy、IDVLEpx專用于以直接計(jì)算法計(jì)算塔頂和塔底再沸器的飽和液相或氣相組成及其操作溫度，而BPT、DPT、Tdx_rv 均為BDPT_2M 返回待求解方程組；rvmean 函數(shù)可根據(jù)給定的方法計(jì)算三種相對揮發(fā)度；BDPT_2M、sxy 可以兩類方法（共四種方法）分別計(jì)算給定組成的泡露點(diǎn)溫度、給定溫度下的飽和氣液相組成。以苯-甲苯體系為例，通過主程序tbtdxeye 調(diào)用相關(guān)函數(shù)可實(shí)現(xiàn)以四種方法計(jì)算泡露點(diǎn)溫度及飽和氣液相組成。通過相對偏差分析發(fā)現(xiàn)，以直接計(jì)算法為基準(zhǔn)，三種相對揮發(fā)度法的整體精度順序?yàn)椋簬缀纹骄?對數(shù)平均法>算術(shù)平均法。三種方法均可計(jì)算二元完全理想系相平衡關(guān)系，可進(jìn)一步用于后續(xù)的精餾計(jì)算中。