《数值分析答辩》由会员分享，可在线阅读，更多相关《数值分析答辩（5页珍藏版）》请在金锄头文库上搜索。
1、常微分方程方法计算连续多组分精馏常微分方程方法计算连续多组分精馏问题背景：问题背景： 精馏是化工中的重要的分离过程，但连续多组分精馏的分析比较复杂，需 要通过解偏微分方程组求得。如采用一般的编程软件解决这个问题，首先编写 或查找相应的 FORTRAN、C 和 BASIC 等语言的程序包，再进行编程、调试， 最后得到结果。这种方法要求计算者对 FORTRAN、C 和 BASIC 等传统的计算 语言具有一定的编程能力，调试过程中花费很多时间和精力。也可以使用商业 软件如 ASPEN PLUS、HYSIS 等，但这些集成软件不需要建模过程，不能利于 用户深层理解，另外当输入的参数不准确时，很难收敛。
2、 MATLAB 相对于以上传统的计算语言有明显的优越性，并且 MATLAB 有 强大的图形演示功能，能图文并茂的展示化工过程，特别是分离过程中各组分 从开始到稳定的变化过程。更直观、更高效、更稳定的模拟化工过程。下面以连续多组分精馏塔从开车到稳定的动态浓度计算为例，简要说明 MATLAB 在化工过程中的应用。 问题描述：连续多组分精馏问题描述：连续多组分精馏 在一个精馏塔中进行对组分液体混合物的分离，该混合物含有苯、苯乙烯 和甲苯三个组分。已知进料流量 F=40kmol/hr，回流比 R=5，进料组成（摩尔分 数）x1=0.6，x2=0.25，塔顶冷凝器中的滞液量 M1=75kmol，塔板滞液
3、量 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）,(1)1,= 2, ( + )1, = 1对精馏塔(column enric
4、hing section)第 i 块塔板的任意组分 j,(2),1= ( 1, ,)+ ( 1, ,) = 2,3, 1对进料板（column feedplate）的任意组分 j,(3),= 1, ,+ + 1, ,+ , = 对提馏段（column stripping section）第 i 块板的任意组分 j,(4),= ( 1, ,)+ ( + , ,) = + 1,+ 2, 1对塔釜（column reboiler）的任意组分 j,(5),= , , , = 在精馏段中 L=RD (6)V=(R+1)D (7) 在提馏段中 L=L+ q F (8)V=V-(1-q)F (9) 若饱和液
5、体进料（泡点进料） ，则 q=1。气液相平衡关系 (10),=, = 1,式中，i 为塔板序号（i=2,-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的各个元
6、素都取为 z1，向量 x2的 各个元素都取为 z2。 matlab 程序：程序： Function ConDistill %连续多组分（三元）精馏塔的模拟计算 clear all clc global F z1 z2 z3 R alpha1 alpha2 alpha3 M1 MN M Nt Nf V1 V D L L1 W 定义全局函数 F=40； %进料流量，kmol/hr R=5; %回流比 Z1=0.6； %苯的进料组成（摩尔分率） Z2=0.25; %甲苯的进料组成（摩尔分率） Z3=1-Z2-Z1； %苯乙烯的进料组成（摩尔分率） M1=75； %塔顶冷凝器中的滞液量（kmol） M
7、=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
8、和 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.
9、-,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) %稳态图 figure plate = 1:Nt; plot(plate,x1(end,:),r.-,plate,x2(end,:),k.-,plate,x3(end,:),b.:) xlabel(塔板) ylabel(稳态时苯，甲苯，苯乙烯的组成)
10、 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 W x1 = 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./den
11、om; 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-1 dx1dt(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
12、*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/h L 液体流量，kmol/h x 液相摩尔分数 M 液体滞留量，kmol/h y 气相摩尔分数 N 塔板数 z 进料组成（摩尔分 数） Q 进料状态 相对挥发度 下标 提馏段
