MATLAB图解精馏塔理论塔板数程序代码

MATLAB图解精馏塔理论塔板数程序代码

function distillation %文件名“distillation”可以更改

% 输入计算所需参数

q=1;%输入进料热状况参数

R=1.5;%输入回流比

xD=0.95;%输入塔顶轻组分摩尔分数

xW=0.04;%输入塔底轻组分摩尔分数

xF=0.52;%输入进料轻组分摩尔分数

%以下输入相平衡数据

x0=[0

0.0196078

0.0392156

0.0588235

0.0784313

0.0980392

0.1176471

0.1372549

0.1568627

0.1764706

0.1960784

0.2156863

0.2352941

0.254902

0.2745098

0.2941176

0.3137255

0.3333333

0.3529412

0.372549

0.3921569

0.4313725 0.4509804 0.4705882 0.4901961 0.5098039 0.5294118 0.5490196 0.5686275 0.5882353 0.6078431 0.627451 0.6470588 0.6666667 0.6862745 0.7058824 0.7254902 0.745098 0.7647059 0.7843137 0.8039216 0.8235294 0.8431373 0.8627451 0.8823529 0.9019608 0.9215686 0.9411765 0.9607843

0.9803922

1];

y0=[0 0.0437029

0.1258286 0.1643911 0.2013788 0.2368595 0.2708994 0.303563 0.3349129 0.3650094 0.3939109 0.4216732 0.4483501 0.4739928 0.4986506 0.5223702 0.5451963 0.5671715 0.5883362 0.6087289 0.6283862 0.6473428 0.6656317 0.6832842 0.70033 0.7167974 0.7327131 0.7481026 0.76299 0.7773982 0.791349 0.8048631 0.8179601 0.8306587

0.8549313

0.8665382

0.877813

0.8887702

0.899424

0.9097874

0.9198734

0.9296939

0.9392607

0.9485847

0.9576768

0.966547

0.9752052

0.9836608

0.9919228

1];

Yr=@(x)R/(R+1).*x+xD/(R+1);%精馏段操作线

fun=@(x)(q-1)*(R/(R+1).*x+xD/(R+1))-(q*(x-xF)+(q-1)*xF); xQ=fzero(fun,0.5);%求操作点

yQ=Yr(xQ);

xOP=[xW,xQ,xD];

yOP=[xW,yQ,xD];

yfit=linspace(0,1,1001);

xfit=interp1(y0,x0,yfit,'pchip');

%%绘制图形

hold on

box on

plot([0 1],[0 1],'k');

xlabel('x')

ylabel('y')

plot(x0,y0,'r')

plot(xfit,yfit,'r-')

plot(xF,xF,'b*')

plot(xQ,yQ,'bo')

plot(xOP,yOP,'b-')

k=1;

yn(1)=xD;

xn(1)=interp1(y0,x0,yn(1),'pchip');

plot([xD,xn(1)],[yn(1),yn(1)],'b-')

text(xn(1),yn(1),num2str(1),...

'HorizontalAlignment','center','VerticalAlignment','bottom') while xn(k)>xW

yn(k+1)=interp1(xOP,yOP,xn(k));

k=k+1;

xn(k)=interp1(y0,x0,yn(k),'pchip');

plot([xn(k-1),xn(k-1)],[yn(k-1),yn(k)],'b-')

plot([xn(k-1),xn(k)],[yn(k),yn(k)],'b-')

text(xn(k),yn(k),num2str(k),...

'HorizontalAlignment','center','VerticalAlignment','bottom' ) end

N=k;

plot([xn(N),xn(N)],[yn(N),xn(N)],'b-')

text(xn(N),yn(N),num2str(N),...

'HorizontalAlignment','center','VerticalAlignment','bottom' )

N_Feed=find(xn

N_Feed=N_Feed(1);

text(0.5,0.5,{strcat('所需理论板：',num2str(N)),...

strcat('进料板位置：',num2str(N_Feed))},...

'HorizontalAlignment','left','VerticalAlignment','top')

%以下代码是为了去掉顶端和右边坐标轴的刻度

box off