数值分析答辩[稻谷书苑].doc
常微分方程方法计算连续多组分精馏问题背景：精馏是化工中的重要的分离过程，但连续多组分精馏的分析比较复杂，需要通过解偏微分方程组求得。如采用一般的编程软件解决这个问题，首先编写或查找相应的 FORTRAN、C和BASIC等语言的程序包，再进行编程、调试，最后得到结果。这种方法要求计算者对 FORTRAN、C和BASIC等传统的计算语言具有一定的编程能力，调试过程中花费很多时间和精力。也可以使用商业软件如ASPEN PLUS、HYSIS等，但这些集成软件不需要建模过程，不能利于用户深层理解，另外当输入的参数不准确时，很难收敛。MATLAB相对于以上传统的计算语言有明显的优越性，并且MATLAB有强大的图形演示功能，能图文并茂的展示化工过程，特别是分离过程中各组分从开始到稳定的变化过程。更直观、更高效、更稳定的模拟化工过程。 下面以连续多组分精馏塔从开车到稳定的动态浓度计算为例，简要说明MATLAB在化工过程中的应用。问题描述：连续多组分精馏在一个精馏塔中进行对组分液体混合物的分离，该混合物含有苯、苯乙烯和甲苯三个组分。已知进料流量F=40kmol/hr，回流比R=5，进料组成（摩尔分数）x1=0.6，x2=0.25，塔顶冷凝器中的滞液量M1=75kmol，塔板滞液量M=10kmol，塔釜中的滞液量MN=150kmol，进料状态q=1(饱和)，相对挥发度：a1=2.75,a2=1,a3=0.4,从塔釜蒸发上来的蒸汽流量V=150 kmol/h，塔板总数Nt=10（包括塔顶冷凝器和塔釜），进料板位置Nf=5。1. 求塔顶和塔釜产品从进料开始直到稳态的动态浓度曲线。2. 画出稳态时精馏塔各塔板上的浓度曲线。3. 研究操作变量（进料流量、进料组成和回流比）的变化对精馏的影响。数学模型对塔顶冷凝器（reflux drum）的任意组分j（这里j=1,2,n）,Midx1,fdt=Vy2,j-L+Dx1,j i=1 (1)对精馏塔(column enriching section)第i块塔板的任意组分j,Mdxi,1dt=Lxi-1,j-xi,j+Vyi-1,j-yi,j i=2,3,Nf-1 (2)对进料板（column feedplate）的任意组分j,Mdxi,jdt=Lxi-1,j-Lxi,j+Vyi+1,j-Vyi,j+Fzi,j i=Nf (3)对提馏段（column stripping section）第i块板的任意组分j,Mdxi,jdt=Lxi-1,j-xi,j+Vyi+a,j-yi,j i=Nf+1,Nf+2,Nt-1(4)对塔釜（column reboiler）的任意组分j,MBdxB,jdt=LxN,j-WxB,j-VyB,j i=Nt (5)在精馏段中 L=RD (6) V=(R+1)D (7)在提馏段中 L=L+ q F (8) V=V-(1-q)F (9)若饱和液体进料（泡点进料），则q=1。气液相平衡关系 yi,j=i,jxi,jj=1ni,jxi,j (10)式中，i为塔板序号（i=2,Nt-1）,j为组分序号(j=1,2,n)。程序说明 用ode45()求解由DistMassBalance()定义的物料平衡方程组（1）（5）。由于是三元组分，各塔板上满足x3=1-x1-x2（其中x1、x2、x3均为向量，分别表示组分1、2和3在各个塔板上的液相摩尔分数），故只需对组分1和2的有关动态方程进行求解。x1 x2初值的选取依据：开车时塔内所有板上的x1和x2分别与进料的z1和z2相同，故初值x1 x2中向量x1的各个元素都取为z1，向量x2的各个元素都取为z2。matlab程序：Function ConDistill%连续多组分（三元）精馏塔的模拟计算clear allclcglobal F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 W 定义全局函数F=40； %进料流量，kmol/hrR=5; %回流比Z1=0.6； %苯的进料组成（摩尔分率）Z2=0.25; %甲苯的进料组成（摩尔分率）Z3=1-Z2-Z1； %苯乙烯的进料组成（摩尔分率）M1=75； %塔顶冷凝器中的滞液量（kmol）M=10； %塔板上的滞液量（kmol）MN=150； %塔釜中的滞液量（kmol）q=1; %饱和进料tf=35; 截止dt=1; 步长%相对挥发度alpha1=2.75; alpha2=1; alpha3=0.4;Nt=10; %塔板总数Nf=5； %进料位置V1=150； %从塔釜蒸发上来的蒸汽流量（kmol/hr）%精馏段V=V1-(1-q)*F； 气相的量D=V/(R+1)； 塔顶产品的量L=V-D； 液相的量%提馏段L1=L+F； 液相的量W=L1-V1； 塔釜产品的量%初始化x1和x2开车时塔内所有板上的x1和x2分别于进料的z1和z2相同x1 = z1 * ones(1,Nt);x2 = z2 * ones(1,Nt);t,y = ode45(DistMassBalances,0:dt:tf,x1 x2) 四阶-五阶Runge-Kutta算法 ode45表示采用四阶-五阶Runge-Kutta算法%输出结果X1=y(:,1:Nt); %苯的液相组成（摩尔分率）X2=y(:,Nt+1:2*Nt); %甲苯的液相组成（摩尔分率）X3=1-x1-x2; %苯乙烯的液相组成（摩尔分率）plot(t,x1(:,1),r-,t,x2(:,1),k-,t,x3(:,1),b:,t,x1(:,end),r.-,t,x2(:,end),k-.,t,x3(:,end),b.-)xlabel(Time (h)ylabel(x_1_1, x_1_2, x_1_3, x_N_1, x_N_2, x_N_3)title(塔顶和塔釜产品从进料开始直至稳态的动态浓度曲线)legend(x_1_1,x_1_2,x_1_3,x_N_1,x_N_2,x_N_3)%稳态图figureplate = 1:Nt;plot(plate,x1(end,:),r.-,plate,x2(end,:),k.-,plate,x3(end,:),b.:)xlabel(塔板)ylabel(稳态时苯，甲苯，苯乙烯的组成)title(稳态时精馏的浓度曲线)legend(苯,甲苯,苯乙烯)% - function dydt = DistMassBalances(t,y) %物料平衡方程组global F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 Wx1 = y(1:Nt); % 组分1（苯）x2 = y(Nt+1:2*Nt); % 组分2（甲苯）x3 = 1-x1-x2; % 组分3（苯乙烯）%气相平衡denom = alpha1*x1+alpha2*x2+alpha3*x3;y1 = alpha1*x1./denom;y2 = alpha2*x2./denom;%对塔顶冷凝器（i=1）i = 1;dx1dt(i) = (V*y1(i+1)-(L+D)*x1(i)/M1; dx2dt(i) = (V*y2(i+1)-(L+D)*x2(i)/M1;%精馏段（i=2Nf-1）for i=2:Nf-1dx1dt(i) = (L*(x1(i-1)-x1(i)+V*(y1(i+1)-y1(i)/M;dx2dt(i) = (L*(x2(i-1)-x2(i)+V*(y2(i+1)-y2(i)/M;end%进料板（i=Nf）i = Nf; dx1dt(i) = (F*z1+L*x1(i-1)-L1*x1(i)+V1*y1(i+1)-V*y1(i)/M;dx2dt(i) = (F*z2+L*x2(i-1)-L1*x2(i)+V1*y2(i+1)-V*y2(i)/M;%提馏段（Nf+1Nt-1）for i=Nf+1:Nt-1 2个循环dx1dt(i) = (L1*(x1(i-1)-x1(i)+V1*(y1(i+1)-y1(i)/M; dx2dt(i) = (L1*(x2(i-1)-x2(i)+V1*(y2(i+1)-y2(i)/M; end%塔釜（i=Nt）i = Nt;dx1dt(i) = (L1*x1(i-1)-V1*y1(i)-W*x1(i)/MN;dx2dt(i) = (L1*x2(i-1)-V1*y2(i)-W*x2(i)/MN;dydt = dx1dt dx2dt;计算结果：（1）塔顶和塔釜产品从进料开始直至稳定的动态浓度曲线如图1所示（经过30h基本达到稳态）；（2）稳态时各塔板上的浓度曲线示于图2中；符号说明：D 塔顶馏出液流量，kmol/h R 回流比F 进料流量，kmol/h V 蒸汽流量，kmol/hL 液体流量，kmol/h x 液相摩尔分数M 液体滞留量，kmol/h y 气相摩尔分数N 塔板数 z 进料组成（摩尔分数）Q 进料状态 相对挥发度下标 提馏段 I 组分i1 塔顶冷凝器 j 组分数B 再沸器 N 塔板总数图1F 进料板 t总数（包括塔顶冷凝器和塔釜）图25思维
