1、.常微分方程方法計算連續多組分精餾問題背景：精餾是化工中的重要的分離過程，但連續多組分精餾的分析比較復雜，需要通過解偏微分方程組求得。如采用一般的編程軟件解決這個問題，首先編寫或查找相應的 FORTRAN、C和BASIC等語言的程序包，再進行編程、調試，最后得到結果。這種方法要求計算者對 FORTRAN、C和BASIC等傳統的計算語言具有一定的編程能力，調試過程中花費很多時間和精力。也可以使用商業軟件如ASPEN PLUS、HYSIS等，但這些集成軟件不需要建模過程，不能利于用戶深層理解，另外當輸入的參數不準確時，很難收斂。MATLAB相對于以上傳統的計算語言有明顯的優越性，并且MATLAB有

2、強大的圖形演示功能，能圖文并茂的展示化工過程，特別是分離過程中各組分從開始到穩定的變化過程。更直觀、更高效、更穩定的模擬化工過程。 下面以連續多組分精餾塔從開車到穩定的動態濃度計算為例，簡要說明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

3、=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,

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

5、)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，

6、向量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； %

7、塔釜中的滯液量（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 = o

8、de45(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_

9、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) %物

10、料平衡方程組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)

11、= (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）fo

12、r 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中；符號說