基于Octave圖解計(jì)算精餾塔理論塔板數(shù)技術(shù)研究

葉世貴，李鵬飛，盧宏濤，文 芳，吳孟琪

(中石化廣元天然氣凈化有限公司，四川 廣元 628400)

為從多組分體系中分離獲得高純度的產(chǎn)品，通常可以采用多級平衡閃蒸的精餾方式，該方式也是化工單元操作中最重要的一種操作方法[1]。多組分精餾塔可分為簡單塔和復(fù)雜塔2種，前者僅有一股進(jìn)料且無側(cè)線出料和中間換熱設(shè)備，在生產(chǎn)中較為常見。對于二元物料精餾體系，經(jīng)常面臨的一類問題是：已知條件為產(chǎn)品濃度，待求量為理論塔板數(shù)。對于二元組分混合物連續(xù)精餾，由于僅含有2種組分，一般可通過逐板計(jì)算法或圖解法求解其理論塔板數(shù)[2]。

1 圖解法計(jì)算原理和過程

1.1 相關(guān)方程

對于連續(xù)精餾塔(見圖1)，塔頂設(shè)置全凝器，泡點(diǎn)回流。其氣液平衡方程可以表示為相對揮發(fā)度α的函數(shù)式：

(1)

根據(jù)物料衡算可以得到精餾段和提餾段塔板的下降液相組成xn(xm)與進(jìn)入該層塔板的氣相組成yn+1(ym+1)關(guān)系式：

(2)

(3)

精餾段操作方程式(2)與提餾段操作方程式3相交得到1個(gè)交點(diǎn)，令其坐標(biāo)為(xq,yq)，則有

(4)

(5)

式4與式5消去參數(shù)xD即可得到交點(diǎn)的軌跡方程，稱為q線方程，其關(guān)系式表示如下：

(6)

1.2 相關(guān)數(shù)據(jù)

用一連續(xù)精餾塔分離乙醇-水混合液，原料液中含有乙醇的摩爾分?jǐn)?shù)xF為0.40，要求塔頂餾出液含乙醇摩爾分?jǐn)?shù)xD為0.78，塔釜含乙醇摩爾分?jǐn)?shù)xW為0.02。塔頂設(shè)全凝器，泡點(diǎn)回流，所用回流比R為3，試用圖解法求進(jìn)料液相分率q為1.103時(shí)的理論板數(shù)和加料板位置。已知乙醇-水的汽-液平衡數(shù)據(jù)見表1。

表1 乙醇-水的汽-液平衡數(shù)據(jù)

1.3 圖解手動計(jì)算過程

1)相平衡線繪制。按照表1所給的平衡數(shù)據(jù)可以做出平衡曲線圖(見圖2)。

4)提餾段操作線繪制。由xW=0.02在對角線上確定點(diǎn)c，連接c、d可得提餾段操作線。

5)階梯線繪制。從a點(diǎn)開始，在平衡線與精餾段操作線之間做梯級并計(jì)數(shù)，當(dāng)梯級與平衡線的交點(diǎn)橫坐標(biāo)小于兩操作線的交點(diǎn)d的橫坐標(biāo)時(shí)，在平衡線與提餾段之間做梯級并繼續(xù)計(jì)數(shù)。當(dāng)最后一個(gè)交點(diǎn)橫坐標(biāo)翻越xW后完成計(jì)數(shù)，由圖2可以得到總理論板數(shù)為7，在第5塊板進(jìn)行加料。

2 Octave自動解算及繪圖過程

2.1 解算及繪圖步驟

1)羅列出所給固定參數(shù)條件：xD=0.78，xW=0.02，xF=0.40，R=3，q=1.103。汽液平衡數(shù)據(jù)賦值給矩陣參數(shù)Eqdata。當(dāng)需要引用所有x值時(shí)，可以用Eqdata(1,:)表示，當(dāng)需要引用所有y值時(shí)，可以用Eqdata(2,:)表示。

2)在直角坐標(biāo)圖上使用plot(xD,xD)、plot(xW,xW)定出塔頂組成a點(diǎn)、塔底組成c點(diǎn)。

3)在直角坐標(biāo)圖上使用函數(shù)plot(Eqdata(1,:),Eqdata(2,:),'*')繪制汽液平衡點(diǎn)，點(diǎn)的呈現(xiàn)形式為*。

4)通過曲線擬合函數(shù)polyfit將散點(diǎn)函數(shù)化，其調(diào)用格式為：p=polyfit(x,y,n)。該函數(shù)采用最小二乘法來對已知數(shù)據(jù)x、y進(jìn)行擬合，擬合多項(xiàng)式階數(shù)為n。返回值p為一行向量，長度為n+1，包含了冪次降序的多項(xiàng)式系數(shù)，即。擬合完成后，可以使用函數(shù)polyval計(jì)算多項(xiàng)式的擬合值，其調(diào)用格式為y=polyval(p,x)。該函數(shù)輸入上面產(chǎn)生的一維向量p，返回值y為在x處的多項(xiàng)式計(jì)算值。如果輸入的自變量x為一矩陣或向量，則返回函數(shù)計(jì)算多項(xiàng)式在x每一個(gè)元素處的擬合值。

5)使用函數(shù)plot(x,y)可以繪制出x在0～1范圍內(nèi)，依據(jù)擬合產(chǎn)生的汽液平衡曲線。

6)使用函數(shù)line繪制對角線。

7)使用fplot函數(shù)繪制精餾段方程所表示的操作線，函數(shù)以x為自變量，R，xD為固定參數(shù)。

8)使用fplot函數(shù)繪制q線方程所表示的操作線，函數(shù)以x為自變量，q，xF為固定參數(shù)。

9)采用數(shù)值解法求精餾段操作線與進(jìn)料線的交點(diǎn)d，調(diào)用的函數(shù)為fzero，其格式為x=fzero(fun,x0)。其中fun為待求根函數(shù)，其具體表現(xiàn)形式為精餾段函數(shù)與q線函數(shù)相減。fun函數(shù)以x為自變量，R，xD，q，xF為固定參數(shù)，x0為求解初值。

10)求解出交點(diǎn)后，使用line繪制交點(diǎn)d與塔底組成c點(diǎn)間的直線，即可得到提餾段操作線。

11)從塔頂組成a點(diǎn)開始，在平衡線和精餾段操作線之間繪制梯級。在梯級跨過交點(diǎn)d時(shí)，轉(zhuǎn)為在平衡線和提餾段操作線之間繪制梯級，直至梯級跨過c點(diǎn)。此過程中形成的一個(gè)三角形梯級表示一塊理論塔板，所繪制的梯級總數(shù)即為總理論塔板數(shù)，過d點(diǎn)的梯級為加料板，最后一個(gè)梯級為塔底再沸器。程序化繪制梯級水平線，其原理是為給定y值求解其平衡x值，通過函數(shù)fzero解平衡多項(xiàng)式即可得到。繪制梯級的垂直線，為給定x值求解操作線上的y值，通過直接將x帶入對應(yīng)的操作線方程求解即可。

2.2 程序代碼

1)設(shè)置主函數(shù)并儲存給定值：

function distillation

clear all;clc;clf;

global xD=0.78 xW=0.02 xF=0.4 R=3 q=1.103;

Eqdata=[

1 0.95 0.894 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.05 0

1 0.942 0.894 0.82 0.755 0.698 0.657 0.614 0.575 0.52 0.43 0.31 0];

2)繪制塔頂組成a點(diǎn)及塔底組成c點(diǎn)：

plot(xD,xD)

text(xD+0.02,xD,'a')

hold on

plot(xW,xW)

text(xW+0.02,xW,'c')

hold on

3)繪制汽液平衡點(diǎn)：

plot(Eqdata(1,:),Eqdata(2,:),'*')

hold on

4)對汽液平衡數(shù)據(jù)進(jìn)行擬合：

p=polyfit(Eqdata(1,:),Eqdata(2,:),9);

5)根據(jù)汽液平衡擬合函數(shù)繪制擬合曲線：

x1=[0.0:0.01:1.0];

y1=polyval(p,x1);

plot(x1,y1)

hold on

6)繪制對角線：

line([0.0,1.0],[0.0,1.0])

7)繪制精餾段操作線：

legend('hide')

recope2=@(x)recope(x,R,xD);

fplot(recope2,[0 0.78],'Color','black')

hold on

8)繪制q線：

qline2=@(x)qline(x,q,xF);

fplot(qline2,[0.4 1],'Color','black')

text(xF+0.02,xF,'f');

hold on

9)求精餾段與進(jìn)料線交點(diǎn)d：

cross2=@(x)cross(x,R,xD,q,xF);

xc=fzero(cross2,0.5);

yc=recope(xc,R,xD);

text(xc+0.02,yc,'d');

10)繪制提餾段操作線：

line([xc,xW],[yc,xW],'Color','black')

11)繪制梯級并顯示總理論塔板數(shù)：

axis square

axis([0 1 0 1])

xlabel('x')

ylabel('y')

x0=xD;

y0=xD;

n=0;

while 1

n=n+1;

[x0,y0]=tri(x0,y0,R,xD,p,xc,yc,xW,n);

if x0

break;

end

end

s=sprintf('Total number:%d',n);

text(0.03,0.9,s)

12)自定義精餾段函數(shù)：

function y=recope(x,R,xD)

y=R*x/(R+1)+xD/(R+1);

13)自定義進(jìn)料線函數(shù)：

function y=qline(x,q,xF)

y=q*x/(q-1)-xF/(q-1);

14)自定義交點(diǎn)函數(shù)：

function y=cross(x,R,xD,q,xF)

y=recope(x,R,xD)-qline(x,q,xF);

15)自定義平衡關(guān)系逆運(yùn)算函數(shù)：

function y=inverseop(x,p,yn)

y=polyval(p,x)-yn;

16)自定義動態(tài)繪制梯級函數(shù)：

function[x,y]=tri(x0,y0,R,xD,p,xc,yc,xW,n)

inverseop2=@(x0)inverseop(x0,p,y0);

x=fzero(inverseop2,x0);

if x>xc

y=recope(x,R,xD);

else

y=yc+(yc-xW)*(x-xc)/(xc-xW);

end

xd=x0;yd=y0;

while xd>x

xd=xd-0.01;

yd=y0;

pause(0.0005)

plot(xd,yd,'r')

hold on

end

while xd<=x & yd>y

xd=x;

yd=yd-0.01;

pause(0.0005)

plot(xd,yd,'r')

hold on

end

text(x-0.03,y0,int2str(n),'VerticalAlignment','bottom','color','red')

2.3 運(yùn)行測試

將上述代碼保存在Octave的運(yùn)行路徑下，如“D:Octavework”，文件取名為distillation.m。在Octave程序命令窗口輸入distillation，即可求解并動態(tài)繪制精餾塔理論塔板數(shù)。運(yùn)行最終結(jié)果如圖3所示，該結(jié)果與手動繪制結(jié)果一致。

2.4 改進(jìn)方向

本程序?yàn)榱搜芯糠奖?，將q值設(shè)置為大于1。由于q值可以表征進(jìn)料加熱狀態(tài)，即q值實(shí)際表示進(jìn)料中飽和液體所占的分率。若進(jìn)料經(jīng)預(yù)熱或部分汽化，其q值較小。當(dāng)沸點(diǎn)進(jìn)料時(shí)，q=1；當(dāng)氣液混合物進(jìn)料時(shí)，q=0.5；當(dāng)進(jìn)料為飽和蒸汽時(shí)，q=0。因此程序在使用q線方程前，應(yīng)該對q值進(jìn)行判斷。1)當(dāng)q=1時(shí)，q線為一條垂直于x軸的直線，其與精餾段的交點(diǎn)d橫坐標(biāo)為xF，縱坐標(biāo)為xF值帶入精餾段操作線方程得出值；2)當(dāng)q=0時(shí)，q線為一條垂直與y軸的直線，其與精餾段的交點(diǎn)縱坐標(biāo)為xF，橫坐標(biāo)為xF值帶入精餾段操作線方程逆向計(jì)算得出值；3)當(dāng)q<0時(shí)，若交點(diǎn)在平衡線上或平衡線外，則無解。

設(shè)計(jì)條件下，如果選用較小的回流比R，則精餾段操作線與提餾段操作線交點(diǎn)往汽液平衡線移動，達(dá)到指定分離程度所需的理論塔板數(shù)增多，當(dāng)回流比降至某一數(shù)值時(shí)，兩操作線的交點(diǎn)d落在平衡線上，此時(shí)也無解。因此在得出交點(diǎn)d坐標(biāo)后，應(yīng)該對該坐標(biāo)的位置進(jìn)行判斷，若交點(diǎn)在平衡線內(nèi)方程有解，否則無解。

3 結(jié)語

利用Octave開源強(qiáng)大的數(shù)據(jù)擬合、數(shù)值計(jì)算和圖表繪制等功能，解決了二元物料體系精餾理論塔板數(shù)圖解計(jì)算問題。程序可移植性強(qiáng)，計(jì)算精度高，省時(shí)省力。自定義動態(tài)繪圖功能，為方便教學(xué)理解圖解過程具有重要的意義。改進(jìn)后的程序，對進(jìn)料加熱狀態(tài)及回流比適應(yīng)性更強(qiáng)。