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

《MATLAB图解精馏塔理论塔板数程序代码.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、MATLAB图解精馏塔理论塔板数程序代码functiondistillation%文件名“distillation”可以更改%输入计算所需参数q=1;%输入进料热状况参数R=1.5;%输入回流比xD=0.95;%输入塔顶轻组分摩尔分数xW=0.04;%输入塔底轻组分摩尔分数xF=0.52;%输入进料轻组分摩尔分数%以下输入相平衡数据x0=[00.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0

2、.0.1];y0=[00.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.700330.0.0.0.762990.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.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

3、,yQ,xD];yfit=linspace(0,1,1001);xfit=interp1(y0,x0,yfit,'pchip');%%绘制图形holdonboxonplot([01],[01],'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(

4、1)],[yn(1),yn(1)],'b-')text(xn(1),yn(1),num2str(1),...'HorizontalAlignment','center','VerticalAlignment','bottom')whilexn(k)>xWyn(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

5、),yn(k)],'b-')text(xn(k),yn(k),num2str(k),...'HorizontalAlignment','center','VerticalAlignment','bottom')endN=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

6、Feed(1);text(0.5,0.5,{strcat('所需理论板：',num2str(N)),...strcat('进料板位置：',num2str(N_Feed))},...'HorizontalAlignment','left','VerticalAlignment','top')%以下代码是为了去掉顶端和右边坐标轴的刻度boxoffax2=axes('Position',get(gca,'Position'),...'XAxisLocation','top',...'YAxisLocation','right',...'C

7、olor','none',...'XColor','k','YColor','k');set(ax2,'YTick',[]);set(ax2,'XTick',[]);boxon运行结果如图1所示：图1图解苯-甲苯精馏塔理论塔板数和进料位置