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

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

2、2549020.27450980.29411760.31372550.33333330.35294120.3725490.39215690.41176470.43137250.45098040.47058820.49019610.50980390.52941180.54901960.56862750.58823530.60784310.6274510.64705880.66666670.68627450.70588240.72549020.7450980.76470590.78431370.80392160.82352940.84313730.86274510.88235290.9019608

3、0.92156860.94117650.96078430.98039221;y0=00.04370290.08562240.12582860.16439110.20137880.23685950.27089940.3035630.33491290.36500940.39391090.42167320.44835010.47399280.49865060.52237020.54519630.56717150.58833620.60872890.62838620.64734280.66563170.68328420.700330.71679740.73271310.74810260.762990.

4、77739820.7913490.80486310.81796010.83065870.84297680.85493130.86653820.8778130.88877020.8994240.90978740.91987340.92969390.93926070.94858470.95767680.9665470.97520520.98366080.99192281;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

5、(xQ);xOP=xW,xQ,xD;yOP=xW,yQ,xD;yfit=linspace(0,1,1001);xfit=interp1(y0,x0,yfit,'pchip');%绘制图形hold onbox onplot(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-

6、9;)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

7、),'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' )endN=k;plot(xn(N),xn(N),yn(N),xn(N),'b-')text(xn(N),yn(

8、N),num2str(N),. 'HorizontalAlignment','center','VerticalAlignment','bottom' )N_Feed=find(xn<xF);N_Feed=N_Feed(1);text(0.5,0.5,strcat('所需理论板：',num2str(N),. strcat('进料板位置：',num2str(N_Feed),. 'HorizontalAlignment','left','VerticalAlignment','top')%以下代码是为了去掉顶端和右边坐标轴的刻度box offax2 = axes('Position',get(gca,'Position'),. 'XAxisLocation','top',. 'YAxisLocation','right',. 'Color','none',. 'XColor