數(shù)值分析答辯

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

2、強(qiáng)大的圖形演示功能，能圖文并茂的展示化工過程，特別是分離過程中各組分從開始到穩(wěn)定的變化過程。更直觀、更高效、更穩(wěn)定的模擬化工過程。 下面以連續(xù)多組分精餾塔從開車到穩(wěn)定的動態(tài)濃度計算為例，簡要說明MATLAB在化工過程中的應(yīng)用。問題描述：連續(xù)多組分精餾在一個精餾塔中進(jìn)行對組分液體混合物的分離，該混合物含有苯、苯乙烯和甲苯三個組分。已知進(jìn)料流量F=40kmol/hr，回流比R=5，進(jìn)料組成（摩爾分?jǐn)?shù)）x1=0.6，x2=0.25，塔頂冷凝器中的滯液量M1=75kmol，塔板滯液量M=10kmol，塔釜中的滯液量MN=150kmol，進(jìn)料狀態(tài)q=1(飽和)，相對揮發(fā)度：a1=2.75,a2=1,a3

3、=0.4,從塔釜蒸發(fā)上來的蒸汽流量V=150 kmol/h，塔板總數(shù)Nt=10（包括塔頂冷凝器和塔釜），進(jìn)料板位置Nf=5。1. 求塔頂和塔釜產(chǎn)品從進(jìn)料開始直到穩(wěn)態(tài)的動態(tài)濃度曲線。2. 畫出穩(wěn)態(tài)時精餾塔各塔板上的濃度曲線。3. 研究操作變量（進(jìn)料流量、進(jìn)料組成和回流比）的變化對精餾的影響。數(shù)學(xué)模型對塔頂冷凝器（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)對進(jìn)料板（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)若飽和液體進(jìn)料（泡點(diǎn)進(jìn)料），則q=1。氣液相平衡關(guān)系 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在各個塔板上的液相摩爾分?jǐn)?shù)），故只需對組分1和2的有關(guān)動態(tài)方程進(jìn)行求解。x1 x2初值的選取依據(jù)：開車時塔內(nèi)所有板上的x1和x2分別與進(jìn)料的z1和z2相同，故初值x1 x2中向量x1的各個元素都取為z1，

6、向量x2的各個元素都取為z2。matlab程序：Function ConDistill%連續(xù)多組分（三元）精餾塔的模擬計算clear allclcglobal F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 W 定義全局函數(shù)F=40； %進(jìn)料流量，kmol/hrR=5; %回流比Z1=0.6； %苯的進(jìn)料組成（摩爾分率）Z2=0.25; %甲苯的進(jìn)料組成（摩爾分率）Z3=1-Z2-Z1； %苯乙烯的進(jìn)料組成（摩爾分率）M1=75； %塔頂冷凝器中的滯液量（kmol）M=10； %塔板上的滯液量（kmol）MN=150； %

7、塔釜中的滯液量（kmol）q=1; %飽和進(jìn)料tf=35; 截止dt=1; 步長%相對揮發(fā)度alpha1=2.75; alpha2=1; alpha3=0.4;Nt=10; %塔板總數(shù)Nf=5； %進(jìn)料位置V1=150； %從塔釜蒸發(fā)上來的蒸汽流量（kmol/hr）%精餾段V=V1-(1-q)*F； 氣相的量D=V/(R+1)； 塔頂產(chǎn)品的量L=V-D； 液相的量%提餾段L1=L+F； 液相的量W=L1-V1； 塔釜產(chǎn)品的量%初始化x1和x2開車時塔內(nèi)所有板上的x1和x2分別于進(jìn)料的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算法%輸出結(jié)果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(塔頂和塔釜產(chǎn)品從進(jìn)料開始直至穩(wěn)態(tài)的動態(tài)濃度曲線)legend(x_1_1,x_1_2,x_1_3,x_N_1,x_N_2,x_N_3)%穩(wěn)態(tài)圖figureplate = 1:Nt;plot(plate,x1(end,:),r.-,plate,x2(end,:),k.-,plate,x3(end,:),b.:)xlabel(塔板)ylabel(穩(wěn)態(tài)時苯，甲苯，苯乙烯的組成)title(穩(wěn)態(tài)時精餾的濃度曲線)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%進(jìn)料板（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個循環(huán)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;計算結(jié)果：（1）塔頂和塔釜產(chǎn)品從進(jìn)料開始直至穩(wěn)定的動態(tài)濃度曲線如圖1所示（經(jīng)過30h基本達(dá)到穩(wěn)態(tài)）；（2）穩(wěn)態(tài)時各塔板上的濃度曲線示于圖2中；符號說