图解法绘塔板图并求塔板数的matlab程序

function N=funmain(a1,b1,a2,b2)

%%a1,b1为精馏段方程y=a1x+b1

%%a2,b2为提馏段方程y=a2x+b2

%%N为总塔板数(不包括再沸器)

syms X Y

x=0:0.0001:1;

x=x';

y1=fun1(x);%%气液平衡方程,可新建函数脚本

y2= a1*x+b1;%%精馏段方程

y3= a2*x+b2;%%提馏段方程

y4=x;%%参照线

x_d=0.7788;%%塔顶浓度

x_f=0.1740;%%进料浓度

x_w=0.0079;%%塔釜浓度

plot(x,y1,'LineWidth',1.5);hold on;plot(x,y2,'r:','LineWidth',0.75);hold on;plot(x,y3,'r:','LineWidth',0.75);hold on;plot(x,y4,'LineWidth',1);

axis([0 1 0 1]);hold on;

[X_1,Y_1]=solve(a1*X+b1-Y==0,Y-X==0,X,Y); %%精馏段方程

[X_2,Y_2]=solve(a2*X+b2-Y==0,Y-X==0,X,Y); %%提馏段方程

[X_3,Y_3]=solve(a1*X+b1-Y==0,a1*X+b1-Y==0,X,Y);

A=[X_1 X_2 X_3;Y_1 Y_2 Y_3];

A=eval(A);

scatter(A(1,:),zeros(1,size(A,2)),'r*');hold on; n=100;%%迭代次数

Y0=zeros(n,1);

X0=zeros(n,1);

Y0(1,1)=x_d;

X0(1,1)=fun1(Y0(1,1));

for i=1:n

X0(i,1)=fun1(Y0(i,1));

Y0(i+1,1)=a1*(X0(i,1))+b1;

if(X0(i,1)

break;

end

end

X0=X0(find(X0~=0));Y0=Y0(find(Y0~=0));

N1=length(find(X0~=0));

x0=X0(N1,1);

A0=zeros(n,1);

B0=zeros(n,1);

A0(1,1)=x0;

B0(1,1)=a2*(x0)+b2;

for i=1:n

B0(i+1,1)=a2*(A0(i,1))+b2;

A0(i+1,1)=fun1(B0(i+1,1));

if(A0(i+1,1)

break;

end

end

A0=A0(find(A0~=0));B0=B0(find(B0~=0)); N2=length(find(A0~=0));

%t=0.0001;%%间距

M=[x_d;X0];

N=[Y0([1:length(X0)],1);B0(1,1)];

O=[B0;x_w];

for i=1:N1

p=linspace(M(i+1),M(i),1000);

q=Y0(i)*ones(1,length(p));

plot(p,q,'LineWidth',2);hold on;

r=linspace(N(i),N(i+1),1000);

s=M(i+1)*ones(1,length(r));

plot(s,r,'LineWidth',2);hold on;

end

for j=1:N2-1

e=linspace(A0(j+1),A0(j),1000);

f=B0(j+1)*ones(1,length(e));

plot(e,f,'LineWidth',2);hold on;

g=linspace(O(j+1),O(j+2),1000);

h=A0(j+1)*ones(1,length(g));

plot(h,g,'LineWidth',2);hold on;

end

plot(x_w*ones(1,1000),linspace(0,0.8,1000),'k:','LineWidth',0.75);hold on; plot(x_d*ones(1,1000),linspace(0,0.8,1000),'k:','LineWidth',0.75);hold on; plot(x_f*ones(1,1000),linspace(0,0.8,1000),'k:','LineWidth',0.75);

xlabel('x');ylabel('y');

N=N1+N2;

end;