化工计算机模拟读书报告

化工过程计算机模拟读书报告题目：利用Matlab确定乙烯加氢制乙烷反应动力学模型参数组员姓名：肖竹钱2111301197杨利2011301202指导教师：夏陆岳专业：化学工程学院：化学工程与材料学院提交日期2013年12月1.绪论乙烷可以在冷冻设施中作为致冷剂使用。在科学研究中液态的乙烷在电子显微技术中被用来使得含水量高的样本透明化。薄的水层在快速沉浸到设施-150度或更低的液乙烷中时迅速冻结，不会形成晶体。这样的迅速冻结不会像冰结晶时那样会破坏液态水中软物质的结构。乙烷(C2H6)是一种低级的烷烃,在生产和生活中有很重要的应用，它能发生很多烷烃的典型反应。乙烷的制备方法之一就是乙烯（C2H4）与氢气(H2)的加成反应，即乙烯和氢气在加热加压催化剂的条件下，乙烯中的碳碳双键断裂变单键，每个碳原子再结合一个氢，生成乙烷．在乙烯与氢气的反应过程求解中,动力学参数的求解往往很关键，但也很复杂。在化学反应的研究中，反应速率随着一些条件的变化而变化，利用一系列变化的数据去分析化学反应往往是非常繁琐的工作，然而利用MATLAB的二乘法函数可以很快捷的估算这些参数，进而求出反应速率，为研究带来了很大的方便。2.问题叙述乙烯（E）加氢（H）制乙烷（EA）反应如下：H2+C2H4→C2H6，反应是在钴钼催化剂上进行的。在微分反应器中测得的实验数据如下表所示。表2-1乙烯加氢动力学实验数据实验号反应速率/[mol/(kgcat·s)]PE/atmPEA/atmPH/atm11.04001.00001.00001.000023.13001.00001.00003.000035.21001.00001.00005.000043.82003.00001.00003.000054.19005.00001.00003.000062.39100.50001.00003.000073.86700.50000.50005.000082.19900.50003.00003.000090.75000.50005.00001.0000根据表中数据确定乙烯加氢制乙烷动力学模型中的k和K.3.问题求解3.1动力学模型乙烯加氢制乙烷反应速率方程为：(1)3.2求解思路先进行线性回归：取（1）式的倒数，并把PE、PH移向左侧，则(2)其中，(3)在b1、b2中含有未知参数k和KE，所以它们成为新的的未知量。按表2-1中的数据，用线性最小二乘法估计（2）式的b1、b2，再由（3）式求出k和KE。然后，由上述线性回归得到的k和KE作为非线性回归的初值，利用非线性最小二乘法的函数lsqnonlin（）直接估计非线性模型（1）中的参数k和KE，并计算其置信区间，即函数图像中哪个区间得到的数据比较可靠。3.3计算结果计算结果：k=3.187±0.288，KE=2.101±0.264；残差平方和为0.042.拟合残差见下图3-2.图3-1程序运行结果图3-2拟合残差4.结果分析通过matlab的运算，最终得到k=3.187±0.288,KE=2.101±0.264;数据点(离散点)分布在回归线两侧，这些点是随机分布的，反应速率一样时，由于压力等其他原因，数据点的分布也会有所不同。这些数据点距离回归线很近，其残差平方和为0.042。由n=9,查表Za/2=1.96,σ=0.2，由置信区间公式求得95%的置信区间为[1.970,2.232]。KE的95%置信区间表明如果进行100次试验，KE的计算值有95次在[1.970,2.232]范围内，即KE=2.101±0.264.由残差平方和越小，拟合程度越好可知，此拟合效果比较满意。对于残差平方和数值，通过matlab程序的优化，可以使其更小，使结果更加的精确。5.结论通过运算，得到了参数k和KE的值。由可知乙烯加氢制乙烷反应速率的计算变的相当简单。利用非线性最小二乘法的函数lsqnonlin（）直接估计非线性模型（1）中的参数k和KE，不仅方便快捷，而且具有很高的精确性。对于matlab的编程而言，也并不复杂。犹如此类的反应都可以用这种方法求解，使化学研究变得快捷高效。附录functionKineticsEst2%动力学参数估计clearallclcrA=[1.043.135.213.824.192.3913.8672.1990.75];Pe=[111350.50.50.50.5];Ph=[135333531];%线性拟合RA=Pe.*Ph./rA;y=RA';X=[ones(size(y))Pe'];b=X\y;k=1/b(1);KE=b(2)*k;%非线性拟合beta0=[kKE][beta,resnorm,residual,exitflag,output,lambda,jacobian]=...lsqnonlin(@ObjFunc,beta0,[],[],[],rA,Pe,Ph)ci=nlparci(beta,residual,jacobian)%残差关于拟合值的残差图rAc=RateA(beta,Pe,Ph)plot(rAc,residual,'*')xlabel('反应速率（-rA）拟合值,mol/(kgcat.s)')ylabel('残差R,mol/(kgcat.s)')refline(0,0)%参数辨识结果fprintf('EstimatedParameters:\n')fprintf('\tk=%.3f±%.3f\n',beta(1),ci(1,2)-beta(1))fprintf('\tKE=%.3f±%.3f\n',beta(2),ci(2,2)-beta(2))fprintf('\tThesumofthesquaresis:%.3f',resnorm)%------------------------------------------------------------------functionf=ObjFunc(beta,rA,Pe,Ph)f=rA-RateA(