空分装置多股流换热器的模拟与性能分析

ABSTRACT-1g）目标函数（4-1a）是由两个部分构成，第一部分由设备两端出口温度的计算值与测量值的误差总和构成，第二部分由设备中间层温度的计算值和同一层有测量点的测量值误差总和构成。由于MHEX复杂的结构特点决定了温度分布形式，在同一个横截面上（图4-1）各热流股之间温度接近，各冷流股之间温度接近，热流股与冷流股之间的温差接近最小温差，为了约束这些温度，基于现场数据观察，本文假定：热流股之间温差为0℃；冷流股之间温差为0℃；最小温差为1℃。目标函数第二部分就是由这些约束构成。公式（4-1b～g）构成约束方程，（4-1b）表示任一段b（Bundle）的回归模型，其中β表示被估参数。传递方程（4-1c）和（4-1d）表示上一层模型的出口温度传递给下一层模型的入口温度。方程（4-1e）和（4-1f）保证同一截面上最小温差1℃。最后，方程（4-1g）表示被估参数的结构特征，例如当热流股流量增加其出口温度随之增大，冷流股流量增加其出口温度减小等。复杂MHEX模型由多个Bundle的回归模型联立构成，每个Bundle由一个简单MHEX模型控制（模型推导见2.2），单个Bundle的回归模型只需将冷热流股数带入公式（2-5）扩展得到，三个Bundle的具体回归模型见公式（4-2），（4-3），（4-4），其中,表示热流股i在第b个Bundle的流量，表示冷流股j在第b个Bundle的流量，表示热流股i在第b个Bundle的进口温度，表示冷流股j在第b个Bundle的进口温度，β表示被估参数的集。表示任意流股i/j在第b个Bundle的流量，表示任意流股i/j在第b个Bundle的出口温度。Bundle1的回归模型由五股流股28个结构参数组成变量组，回归结构如下， （4-2）Bundle2的回归模型由六股流股39个结构参数组成变量组，回归结构如下， （4-3）Bundle3的回归模型由五股流股28个结构参数组成变量组，回归结构如下， （4-4）4.2多层出口的复杂多股流换热器模型求解与结果4.2.1多层出口的复杂多股流换热器模型求解按上述方法对空分装置主换热器和过冷器建模，建立一个多Bundle联立回归模型来模拟冷量回收的过程。复杂MHEX在GeneralAlgebraicModelingSystem(GAMS)中建立模型，采用联立方程法求解。回归问题被再构成一个二次优化问题，搜索一组最优的结构参数，使得出口温度的计算与测量值误差最小。采用CONOPT求解器求解非线性问题，具体的求解程序见第四章末尾的程序附录。我们需要一组结构参数的初值去搜索联立方程的最优解。初值很大程度上影响计算时间，而且差的初值很可能导致局部最优解而非全局最优解。因此，初值的求解是模型求解的第一步工作，这是一个迭代求解的过程见图4-2，全部在Matlab@R2011b进行。图4-SEQ图4-\*ARABIC2复杂模型初值计算过程首先将每一个Bundle的训练数据做简化处理，假定中间层热流股温度相同，冷流股温度相同，且热流股与冷流股温差为1℃，利用Matlab数据库中的Regress函数先对Bundle1进行最小二乘回归得到一组结构参数，用结构参数和入口温度重新计算Bundle1出口温度，将热流股温度传递给Bundle2的入口温度，同样的方法回归Bundle2结构，计算Bundle2热流股出口温度传递热流给Bundle3，回归Bundle3结构参数计算冷流股出口温度传递回Bundle2，依次按照图4-2的顺序迭代。值得一提的是与其让迭代达到局部最优解，不如仅仅只迭代最主要的5次，得到三个Bundle模型的结构参数作为联立方程求解的初值。在求解初值时，对于没有测量的数据，根据MHEX传热结构的特征，假定在MHEX同一截面上没有测量点的热流股温度近似等于有测量点的热流股温度，同理冷流股亦是，热流股的温度近似比冷流股高1℃，流股流量由一个Bundle传递到另一个Bundle,只要流股没有变化，流量也等于前一个Bundle的流量。初值求解的结果包括每一个Bundle中每一股流股的结构参数，进口温度，出口温度和流量。4.2.2多层出口的复杂多股流换热器模型求解结果与讨论将初值传递给GAMS，开始联立模型的求解，主换热器的回归模型计算如表4-1所示，复杂模型包含了39,423个变量，100,501个约束（其中等式约束为37,500个），用时253分钟在1500训练样本中搜索到最优的结构参数，Bundle1由两股热流股和三股冷流股组成，回归得到5组，每组28结构参数。Bundle2由三股热流股和三股冷流股组成，回归得到6组，每组39结构参数。Bundle1由两股热流股和三股冷流股组成，回归得到5组，每组28结构参数。表4-SEQ表4-\*ARABIC1模型计算统计变量约束回归时间计算时间计算精度多层出口复杂MHEX模型39,423100,501253分钟17秒2%简单MHEX模型5070.12秒0.00138秒1%然后再用结构参数求解1000组测试样本的出口温度，此过程亦是一个联立计算的过程，用时17秒。将出口温度的计算值与测试值比较来检验模型估算的准确性。对于训练样本，最优结构参数在多Bundle模型中被固定去测试出口温度的计算值和测量值。当结构参数在测试模型中被固定，则模型自由度为0，变成一个线性模拟问题。虽然多层出口的复杂模型在计算上确实比第二章的简单模型要费时的多，在计算进度上也不如简单模型，但是多层出口的复杂模型更合乎流股换热的内在关系，更能反映实际换热过程，且能够处理某些流股缺乏仪表测量的情况。用CPLEX求解器只需要几秒钟就能求解出线性模型的解，训练样本和测试样本的误差见图4-3。图4-3画出了七个流股训练数据出口温度的回归误差与测试数据出口温度的模拟误差，图中上层实线表示出口温度的误差，对应左边Y轴坐标，下层实线表示主流空气H3的进口温度（记TMinH3），对应右边Y轴坐标，图中虚线左边是训练样本，右边是测试样本。从图中可以看出，训练样本的误差基本分布在-2℃到2℃之间，个别误差大于2℃，最大最对误差小于3℃。误差较为均匀的分布在0轴的两侧，表明回归模型得到了很好的训练。其次我们发现训练样本的大误差大多发生在主流空气H3的进口温度剧烈波动时，其原因同样是分子筛周期性切换引起。当空分换热系统偏离正常稳态操作的一段时间内或者说面临工况切换时，会导致了回归模型的不准确性。如果只考虑稳态模拟，就不应该将这些变负荷工况动态切换过程的操作数据纳入回归范围，这将明显提高多层出口的复杂MHEX模型的回归精度，若是考虑变负荷操作动态切换过程，则出于对测量数据不准确性的考虑，应该对测量数据做标定预处理。图4-SEQ图4-\*ARABIC3回归结果误差分析复杂MHEX模型回归和测试的平均相对误差和最大绝对误差见表4-2。从图4-3和表4-2中可以发现，虽然考虑多层出口的复杂MHEX模型训练结果和测试结果的最大误差和平均误差都比第二章考虑单层出口的MHEX模型要大一些，一方面是由于目前模型的约束条件做了简化处理还不够完善，缺乏多股流传热结构的约束，另一方面是模型缺乏一个更高级的算法。但是考虑多层出口的复杂MHEX模型应用性更普遍，是对只考虑单层出口模型的简单MHEX模型的扩展，若模型只有单层出口对应复杂模型只有一个Bundle的情况，若模型流股数众多，有多层出口，且某些流股缺乏测量数据时，考虑单层出口模型的简单MHEX模型将不再适用，但是由于考虑多层出口的复杂MHEX模型是一个优化命题，可以通过条件约束来解决。表4-SEQ表4-\*ARABIC2训练与测试误差分析H1H2H3C1C2C3C4TrainingdataMean*0.2910.0900.1330.1350.1470.3680.135Max*1.9220.5280.6330.9460.9182.9721.199CheckingdataMean0.4990.1550.2180.1450.1410.5870.172Max2.5250.5480.7231.1470.9362.9851.502*Intable1,“Mean”referstomeanabsoluteerror,“Max”referstomaxabsoluteerror.4.3小结本章考虑流股常常从MEHX中间进入或离开，改进建立了考虑多层出口的复杂MHEX联立捷算模型，这个模型不再用传统处理黑箱的方法将整个MHEX看成一个黑箱，而是将MHEX在中间进出的位置分成多段，每一段记一个Bundle，每个Bundle是一个黑箱由一个模型控制。复杂MHEX模型由多个Bundle的回归模型联立构成，由于有中间变量的存在，不再是输入变量的显函数。模型的预测平均误差不大于2%，平均绝对误差小于0.6℃。此模型更符合设备结构的传热模式且能解决某些流股缺乏测量数据的问题。程序附录*MHEXregressionSETSIflow/17,18,19,21,23,25,32/BBundle/1*3/INOUT1(I)in(out)flowofeachBundle/17*19,21,23/INOUT2(I)in(out)flowofeachBundle/17*19,21,25,32/INOUT3(I)in(out)flowofeachBundle/18,19,21,25,32/HOThotflow/17,18,25/COLDcoldflow/19,21,23,32/Ndata/1*120/M1numberregressionstructure1/1*29/M2numberregressionstructure2/1*40/M3numberregressionstructure3/1*29/;PARAMETERSTMin(I,B,N)inlettemperaturemeasurementTMout(I,B,N)outlettemperaturemeasurementFM(I,N)flowratemeasurementVARIABLESTCin(I,B,N)calculatedtemperatureTCout(I,B,N)calculatedtemperatureoutK1(INOUT1,M1)regressioncoefficientsforeachBundleK2(INOUT2,M2)regressioncoefficientsforeachBundleK3(INOUT3,M3)regressioncoefficientsforeachBundleZdummyobjectiveEQUATIONSREGRESSION1(INOUT1,N)regressionfunctionsREGRESSION2(INOUT2,N)regressionfunctionsREGRESSION3(INOUT3,N)regressionfunctionsCONNECTIONS1(N)connectionequationsCONNECTIONS2(N)connectionequationsCONNECTIONS3(N)connectionequationsCONNECTIONS4(N)connectionequationsCONNECTIONS5(N)connectionequationsCONNECTIONS6(N)connectionequationsCONNECTIONS7(N)connectionequationsCONNECTIONS8(N)connectionequationsCONNECTIONS9(N)connectionequationsOBJdummyobjectivefunciton;OBJ..z=e=sum(N,sqr(TCout("19","1",N)/TMout("19","1",N)-1)+sqr(TCout("21","1",N)/TMout("21","1",N)-1)+sqr(TCout("23","1",N)/TMout("23","1",N)-1)+sqr(TCout("17","2",N)/TMout("17","2",N)-1)+sqr(TCout("32","2",N)/TMout("32","2",N)-1)+sqr(TCout("18","3",N)/TMout("18","3",N)-1)+sqr(TCout("25","3",N)/TMout("25","3",N)-1));REGRESSION1(INOUT1,N)..Fm(INOUT1,N)*TCout(INOUT1,"1",N)=E=K1(INOUT1,"1")+K1(INOUT1,"2")*FM("17",N)*FM("17",N)+K1(INOUT1,"3")*FM("18",N)*FM("18",N)+K1(INOUT1,"4")*FM("19",N)*FM("19",N)+K1(INOUT1,"5")*FM("21",N)*FM("21",N)+K1(INOUT1,"6")*FM("23",N)*FM("23",N)+K1(INOUT1,"7")*FM("17",N)*FM("19",N)+K1(INOUT1,"8")*FM("17",N)*FM("21",N)+K1(INOUT1,"9")*FM("17",N)*FM("23",N)+K1(INOUT1,"10")*FM("18",N)*FM("19",N)+K1(INOUT1,"11")*FM("18",N)*FM("21",N)+K1(INOUT1,"12")*FM("18",N)*FM("23",N)+K1(INOUT1,"13")*FM("17",N)*TMin("17","1",N)+K1(INOUT1,"14")*FM("18",N)*TMin("18","1",N)+K1(INOUT1,"15")*FM("19",N)*TCin("19","1",N)+K1(INOUT1,"16")*FM("21",N)*TCin("21","1",N)+K1(INOUT1,"17")*FM("23",N)*TMin("23","1",N)+K1(INOUT1,"18")*FM("17",N)*TCin("19","1",N)+K1(INOUT1,"19")*FM("17",N)*TCin("21","1",N)+K1(INOUT1,"20")*FM("17",N)*TMin("23","1",N)+K1(INOUT1,"21")*FM("18",N)*TCin("19","1",N)+K1(INOUT1,"22")*FM("18",N)*TCin("21","1",N)+K1(INOUT1,"23")*FM("18",N)*TMin("23","1",N)+K1(INOUT1,"24")*FM("19",N)*TMin("17","1",N)+K1(INOUT1,"25")*FM("19",N)*TMin("18","1",N)+K1(INOUT1,"26")*FM("21",N)*TMin("17","1",N)+K1(INOUT1,"27")*FM("21",N)*TMin("18","1",N)+K1(INOUT1,"28")*FM("23",N)*TMin("17","1",N)+K1(INOUT1,"29")*FM("23",N)*TMin("18","1",N);REGRESSION2(INOUT2,N)..Fm(INOUT2,N)*TCout(INOUT2,"2",N)=E=K2(INOUT2,"1")+K2(INOUT2,"2")*FM("17",N)*FM("17",N)+K2(INOUT2,"3")*FM("18",N)*FM("18",N)+K2(INOUT2,"4")*FM("25",N)*FM("25",N)+K2(INOUT2,"5")*FM("19",N)*FM("19",N)+K2(INOUT2,"6")*FM("21",N)*FM("21",N)+K2(INOUT2,"7")*FM("32",N)*FM("32",N)+K2(INOUT2,"8")*FM("17",N)*FM("19",N)+K2(INOUT2,"9")*FM("17",N)*FM("21",N)+K2(INOUT2,"10")*FM("17",N)*FM("32",N)+K2(INOUT2,"11")*FM("18",N)*FM("19",N)+K2(INOUT2,"12")*FM("18",N)*FM("21",N)+K2(INOUT2,"13")*FM("18",N)*FM("32",N)+K2(INOUT2,"14")*FM("25",N)*FM("19",N)+K2(INOUT2,"15")*FM("25",N)*FM("21",N)+K2(INOUT2,"16")*FM("25",N)*FM("32",N)+K2(INOUT2,"17")*FM("17",N)*TCin("17","2",N)+K2(INOUT2,"18")*FM("18",N)*TCin("18","2",N)+K2(INOUT2,"19")*FM("25",N)*TMin("25","2",N)+K2(INOUT2,"20")*FM("19",N)*TCin("19","2",N)+K2(INOUT2,"21")*FM("21",N)*TCin("21","2",N)+K2(INOUT2,"22")*FM("32",N)*TCin("32","2",N)+K2(INOUT2,"23")*FM("17",N)*TCin("19","2",N)+K2(INOUT2,"24")*FM("17",N)*TCin("21","2",N)+K2(INOUT2,"25")*FM("17",N)*TCin("32","2",N)+K2(INOUT2,"26")*FM("18",N)*TCin("19","2",N)+K2(INOUT2,"27")*FM("18",N)*TCin("21","2",N)+K2(INOUT2,"28")*FM("18",N)*TCin("32","2",N)+K2(INOUT2,"29")*FM("25",N)*TCin("19","2",N)+K2(INOUT2,"30")*FM("25",N)*TCin("21","2",N)+K2(INOUT2,"31")*FM("25",N)*TCin("32","2",N)+K2(INOUT2,"32")*FM("19",N)*TCin("17","2",N)+K2(INOUT2,"33")*FM("19",N)*TCin("18","2",N)+K2(INOUT2,"34")*FM("19",N)*TMin("25","2",N)+K2(INOUT2,"35")*FM("21",N)*TCin("17","2",N)+K2(INOUT2,"36")*FM("21",N)*TCin("18","2",N)+K2(INOUT2,"37")*FM("21",N)*TMin("25","2",N)+K2(INOUT2,"38")*FM("32",N)*TCin("17","2",N)+K2(INOUT2,"39")*FM("32",N)*TCin("18","2",N)+K2(INOUT2,"40")*FM("32",N)*TMin("25","2",N);REGRESSION3(INOUT3,N)..Fm(INOUT3,N)*TCout(INOUT3,"3",N)=E=K3(INOUT3,"1")+K3(INOUT3,"2")*FM("18",N)*FM("18",N)+K3(INOUT3,"3")*FM("25",N)*FM("25",N)+K3(INOUT3,"4")*FM("19",N)*FM("19",N)+K3(INOUT3,"5")*FM("21",N)*FM("21",N)+K3(INOUT3,"6")*FM("32",N)*FM("32",N)+K3(INOUT3,"7")*FM("18",N)*FM("19",N)+K3(INOUT3,"8")*FM("18",N)*FM("21",N)+K3(INOUT3,"9")*FM("18",N)*FM("32",N)+K3(INOUT3,"10")*FM("25",N)*FM("19",N)+K3(INOUT3,"11")*FM("25",N)*FM("21",N)+K3(INOUT3,"12")*FM("25",N)*FM("32",N)+K3(INOUT3,"13")*FM("18",N)*TCin("18","3",N)+K3(INOUT3,"14")*FM("25",N)*TCin("25","3",N)+K3(INOUT3,"15")*FM("19",N)*TMin("19","3",N)+K3(INOUT3,"16")*FM("21",N)*TMin("21","3",N)+K3(INOUT3,"17")*FM("32",N)*TMin("32","3",N)+K3(INOUT3,"18")*FM("18",N)*TMin("19","3",N)+K3(INOUT3,"19")*FM("18",N)*TMin("21","3",N)+K3(INOUT3,"20")*FM("18",N)*TMin("32","3",N)+K3(INOUT3,"21")*FM("25",N)*TMin("19","3",N)+K3(INOUT3,"22")*FM("25",N)*TMin("21","3",N)+K3(INOUT3,"23")*FM("25",N)*TMin("32","3",N)+K3(INOUT3,"24")*FM("19",N)*TCin("18","3",N)+K3(INOUT3,"25")*FM("19",N)*TCin("25","3",N)+K3(INOUT3,"26")*FM("21",N)*TCin("18","3",N)+K3(INOUT3,"27")*FM("21",N)*TCin("25","3",N)+K3(INOUT3,"28")*FM("32",N)*TCin("18","3",N)+K3(INOUT3,"29")*FM("32",N)*TCin("25","3",N);CONNECTIONS1(N)..TCout("17","1",N)=E=TCin("17","2",N);CONNECTIONS2(N)..TCout("18","1",N)=E=TCin("18","2",N);CONNECTIONS3(N)..TCin("19","1",N)=E=TCout("19","2",N);CONNECTIONS4(N)..TCin("21","1",N)=E=TCout("21","2",N);CONNECTIONS5(N)..TCin("19","2",N)=E=TCout("19","3",N);CONNECTIONS6(N)..TCin("21","2",N)=E=TCout("21","3",N);CONNECTIONS7(N)..TCin("32","2",N)=E=TCout("32","3",N);CONNECTIONS8(N)..TCout("18","2",N)=E=TCin("18","3",N);CONNECTIONS9(N)..TCout("25","2",N)=E=TCin("25","3",N);ModelMHEX/all/;*===ImportdatafromExcel*TMin$CALLGDXXRW.EXEktest.xlspar=TMin1rng=TMin!B3:I123ParameterTMin1(N,I);$GDXINktest.gdx$LOADTMin1$GDXIN$CALLGDXXRW.EXEktest.xlspar=TMin2rng=TMin!J3:Q123ParameterTMin2(N,I);$GDXINktest.gdx$LOADTMin2$GDXIN$CALLGDXXRW.EXEktest.xlspar=TMin3rng=TMin!R3:Y123ParameterTMin3(N,I);$GDXINktest.gdx$LOADTMin3$GDXIN*===SetTMin(I,B,N)asTMin1,2,3TMin(I,"1",N)=TMin1(N,I);TMin(I,"2",N)=TMin2(N,I);TMin(I,"3",N)=TMin3(N,I);*TMout$CALLGDXXRW.EXEktest.xlspar=TMout1rng=TMout!B2:I122ParameterTMout1(N,I);$GDXINktest.gdx$LOADTMout1$GDXIN$CALLGDXXRW.EXEktest.xlspar=TMout2rng=TMout!J2:Q122ParameterTMout2(N,I);$GDXINktest.gdx$LOADTMout2$GDXIN$CALLGDXXRW.EXEktest.xlspar=TMout3rng=TMout!R2:Y122ParameterTMout3(N,I);$GDXINktest.gdx$LOADTMout3$GDXIN*===SetTMout(I,B,N)asTMout1,2,3TMout(I,"1",N)=TMout1(N,I);TMout(I,"2",N)=TMout2(N,I);TMout(I,"3",N)=TMout3(N,I);*FM(I,B,N)$CALLGDXXRW.EXEktest.xlspar=FM1rng=FM!B1:I121ParameterFM1(N,I);$GDXINktest.gdx$LOADFM1$GDXINFM(I,N)=FM1(N,I);*K1,K2,K3$CALLGDXXRW.EXEktest.xlspar=K10rng=K!B1:AE6ParameterK10(INOUT1,M1);$GDXINktest.gdx$LOADK10$GDXINK1.l(INOUT1,M1)=K10(INOUT1,M1);$CALLGDXXRW.EXEktest.xlspar=K20rng=K!B7:AP13ParameterK20(INOUT2,M2);$GDXINktest.gdx$LOADK20$GDXINK2.l(INOUT2,M2)=K20(INOUT2,M2);$CALLGDXXRW.EXEktest.xlspar=K30rng=K!B14:AE19ParameterK30(INOUT3,M3);$GDXINktest.gdx$LOADK30$GDXINK3.l(INOUT3,M3)=K30(INOUT3,M3);*optionsNLP=IPOPT;optionsreslim=3600;solveMHEXusingNLPminimizingz;displayK1.l,K2.l,K3.l*savesolutionexecute_unload"k.gdx"K1.lK2.lK3.lexecute'gdxxrw.exek.gdxvar=K1.lrng=B1:AE6'execute'gdxxrw.exek.gdxvar=K2.lrng=B7:AP13'execute'gdxxrw.exek.gdxvar=K3.lrng=B14:AE19'第五章总结与展望5.1总结在实际操作中寻找多股流换热器的最佳工作点是空分装置节能的重要手段，也是过程实时优化的研究方向。本文以20000Nm3/h内压缩空分装置的MHEX为代表，对结构复杂，换热温差小，耦合程度高的MHEX进行了建模、求解和传热性能研究，并对AspenPlus平台上的空分流程模拟进行了补充完善。本文的研究工作和成果包括：1）本文结合传热机理和数据回归模型给出模型自变量的选取规则，基于数据回归方法，建立MHEX的操作数据模型。模型利用进出口温度、和流量回归出口温度和进口温度、流量的关系，得到一组回归参数。用于模拟的模型只有进出变量，没有中间变量，是输入变量的显函数。这模型求解出口温度时计算速度开，收敛性好，计算得到的出口温度与实际出口温度相对误差小于1%，主换热器平均绝对误差小于0.35℃，过冷器平均绝对误差小于0.05℃。2）本文采用的机理分析从几种评价方式中选取复热不足冷损做MHEX传热性能的表征方式，从数据角度验证了选择的正确性。关于原料空气量做一维标度化分析，复热不足冷损比冷量回收率关于空气总流量有更为明显的灵敏度关系。采用MHEX模型和PCA做高维变量的数据处理，PCA通过特征提取，将多股流换热器的高维变量转化成相互独立的低维变量，实现变量降维的目的。PCA结果表明主元一和主元二包含了原变量94%的信息。空分装置现场控制时，复热不足冷损在工程上直观表现为热端温差，MHEX的相关变量应尽量控制在主元一较小的区域。3）本文采用Fortran语言编译了MHEX用户自定义模型（Custommodular），实现了AspenPlus平台的MHEX单元操作模拟。然后替换内压缩空分流程模拟中的MHeatX，通过灵敏度分析调和误差大的测量数据，实现全流程的操作模拟，主要参数的模拟误差基本小于5%，为空分变负荷优化奠定了先决条件。4）