matlab-物料衡算

1、工程应用数学作业5材研1401 642080502011 黄惠洪题目：1.用不同温度下固体铅的定压热容数据，求298K时的绝对熵。2.气液平衡计算。第一题：实验数据的数值积分问题用不同温度下固体铅的定压热容数据，求298K时的绝对熵。数据见下表：由热力学方法可知，绝对熵与定压热容有如下关系：本题可用多种方法来解决，如Simpson数值积分法，或采用的形式，或先用数据点拟合方程得回归参数，然后对回归方程进行积分（回归方程可用，也可在此基础上探索更好的回归方程形式）。1方法与思路采用变型公式：求解给定区间数值积分问题。2. 过程与结果1）积分区间的确定：根据网上教程，matlab对于(-, +)的

2、定义是x=-inf: inf，但是正常情况下是无需这样定义的，因系统默认变量的取值是从负无穷到正无穷的。2 ) 编写程序Entropy_Int1.m并保存：这里将最小二乘B样条拟合法和三次样条插值函数计算作比较，命令行程序如下：function Entropy_Intclear all, clc% 读入数据T = 5. 10. 15. 20. 25. 30. 50. 70. 100. 150. 200. 250. 298.;Cp = 0.305 2.8 7.0 10.8 14.1 16.5 21.4 23.3 24.5 25.4 25.8 26.2 26.5;lnT=log(T);% 用最小二

3、乘B样条拟合法和三次样条插值函数计算lnTi = linspace(lnT(1),lnT(end),299);knots = 3; K = 3; % 三次B样条sp = spap2(knots,K,lnT,Cp);cs = csapi(lnT,Cp); % 生成三次样条插值函数cs % 绘制浓度拟合曲线Ti = linspace(lnT(1),lnT(end),299);Cp_B = fnval(sp,Ti);Cp_c = fnval(cs,Ti);plot(lnT,Cp,'ro',Ti,Cp_B,'b-',Ti,Cp_c,'Bl-.')xlab

4、el('LnT')ylabel('Cp')legend('实验值','B样条拟合','立方样条拟合')% 进行数值积分pp1= fnint(sp);s1= fnval(pp1,lnT)pp2= fnint(cs);s2= fnval(pp2,lnT)end命令窗口输入Entropy_Int，结果如下。s1 = Columns 1 through 11 0 1.1069 3.3382 5.9047 8.5613 11.2145 20.9711 28.7027 37.3706 47.4362 54.7432 Colum

5、ns 12 through 13 60.5284 65.1646s2 = Columns 1 through 11 0 0.6604 2.5965 5.1400 7.9193 10.7146 20.5046 28.0510 36.5924 46.7244 54.0893 Columns 12 through 13 59.8897 64.51873.结果分析1）查阅CRC化学和物理手册，可得到固体铅298K时的绝对熵为64.8，与计算结果之间的差别很小。s1=(s1-s)/s=0.5%，s2=(s2-s)/s=0.4%。2）对比B样条拟合和立方样条拟合，可以发现，立方样条拟合和实验值的匹配度更高

6、。最后拟合出的结果也是立方样条更接近绝对熵的理论值。但是B样条拟合可以适当的曾大结点数knots来提高其拟合程度。3）其中，T=0K时的实验数据由于等于-，如果采用-inf带入求解会出现如下警告：Warning: All data points with NaN or Inf as their site will be ignored. 说明-将被系统忽略。因此直接将此点舍去来让程序可以正常运行。便可计算出结果。第二题：气液平衡计算准确的汽液平衡数据可以用来计算液相活度系数和Gibbs过量自由能。对下表的双组分汽液平衡数据，计算液相活度系数，并对结果进行热力学的一致性检验，假定气相是理想的，即

7、 ，其中 为纯组分的饱和蒸汽压。数据来源：S. Weissman and S. E. Wood, "Vapor-Liquid Equilibrium of Benzene-2,2,4-trimethylpentane Mixtures," J. Chem. Phys., vol.32, 1960, p.1153表51汽液平衡数据(液、汽相摩尔分率为组分1：苯)实验点计算值液相摩尔分率汽相摩尔分率平衡总压P(mmHg)-0.00.0178.081.4411.0-0.3650.08190.1869202.741.4151.008-0.3386-0.41180.21920.406

8、5236.861.3431.011-0.2840-0.48370.35840.5509266.041.2501.046-0.1787-0.54210.38310.5748270.731.2421.048-0.1670-0.56950.52560.6786293.361.1581.116-0.0367-0.65180.84780.8741324.661.0231.5080.38760.57150.98720.9863327.391.0001.9680.67670.69471.01.0327.051.02.0260.706 对低压汽液平衡，假定气相为理想是合理的，此时有：由上表数据和公式，结合由数

9、据手册上得到的饱和蒸汽压数据或计算数据，可得到数据点对应的活度系数，见上表。对汽液平衡数据进行热力学一致性检验要用Gibbs-Duhem方程。对双组分情况，Gibbs-Duhem方程的形式是：如采用面积检验法，有：积分的区间是0,1，因此还缺少两端点，即无限稀释的活度系数，需采用估计计算。无限稀释的活度系数的定义为：和 在分离过程设计和汽液平衡活度系数模型的参数估计中都需要无限稀释的活度系数，但往往很难用实验测定，因此它的计算是很有用的。这里介绍一个估计计算无限稀释活度系数的方法。首先计算和-，然后将它们和一起对作图，可得下图。分别将-和曲线外推到，然后将它们和，其中用接近零的液相摩尔分率和活

10、度系数计算，这样计算得到的两端点的也列于上表，并可作出整个区间0,1的待积分图形。1方法与思路1）首先计算、-和，对、-进行多项式拟合，得到端点值，从而确定边值条件。2）对所有数据点进行三次样条插值，对汽液平衡数据进行热力学一致性检验。2. 过程与结果编写程序VLE.m并保存：function VELclear all,clc%确定边界条件%导入数据，不含端点x1=0.0819,0.2192,0.3584,0.3831,0.5256,0.8478,0.9872;y1=0.1869,0.4065,0.5509,0.5748,0.6786,0.8741,0.9863;p=202.74,236.86

11、,266.04,270.73,293.36,324.66,327.39;p10=327.05;p20=178.08;r1=(p.*y1)./(p10.*x1);r2=(p.*(1-y1)./(p20.*(1-x1);t1=log(r2)./(x1).2);t2=-log(r1)./(1-x1).2);t3=log(r2./r1);x15=0.0819,0.2192,0.3584,0.3831,0.5256;T15=t2(1) t2(2) t2(3) t2(4) t2(5);x67=0.8478,0.9872;T67=t1(6) t1(7);%多项式拟合P1=polyfit(x15,T15,1)

12、;P2=polyfit(x67,T67,1);%端值计算t1end=polyval(P2,1-1e-100);t2int=polyval(P1,1e-100);%立方样条差值拟合t1int=log(1)./(1e-100).2);t1=t1 t1end;t1=t1int t1;t2end=-log(1)./(1e-100).2);t2=t2int t2;t2=t2 t2end;t3=t2int t3 t1end;X1=1e-100 x1 1-1e-100;cs=csapi(X1,t3); X2=1e-100 x1;X3=x1 1-1e-100;X1all=linspace(X1(1),X1(e

13、nd),100);cs_X1=fnval(cs,X1all);Y1=1e-100 x15;Z1=polyval(P1,Y1);Y2=x67 1-1e-100;Z2=polyval(P2,Y2);%端点活度系数Vint=exp(-polyval(P1,0)Vend=exp(polyval(P2,1)t3int=log(1/Vint)t3end=log(Vend)%绘制气液平衡曲线与多项式拟合曲线plot(X1all,cs_X1,'b-',X1,t3,'b*',Y1,Z1,'r-o',Y2,Z2,'g-o')xlabel('x

14、1')ylabel('%')legend('ln(2/1)立方样条插值拟合','ln(2/1)实验数据点','ln2/x12','-ln1/x22')%进行数值积分pp=fnint(cs);s=fnval(pp,X1)命令窗口输入VEL，结果如下：t1 = 0 1.2283 0.2280 0.3479 0.3185 0.3975 0.5716 0.6945 0.7058t2 = -0.3646 -0.4115 -0.4838 -0.5428 -0.5695 -0.6522 -1.0022 -0.7750 0

15、t3 = -0.3646 -0.3386 -0.2840 -0.1787 -0.1700 -0.0370 0.3876 0.6767 0.7058Vint = 1.4400Vend = 2.0255t3int = -0.3646t3end = 0.7058s = 0 -0.0287 -0.0721 -0.1039 -0.1082 -0.1242 -0.0733 0.0002 0.00903.结果分析1) 将对应计算结果填入下表表51汽液平衡数据(液、汽相摩尔分率为组分1：苯)实验点计算值液相摩尔分率汽相摩尔分率平衡总压P(mmHg)-0.00.0178.081.4411.0-0.3646-0.36460.08190.1869202.741.4151.008-0.3386-0.41151.22830.21920.4065236.861.3431.011-0.2840-0.48380.22800.35840.5509266.041.2501.046-0.1787-0.54280.34790.38310.5748270.731.2421.048-0.1700-0.56950.318