基于稳定性的循环物流系统流程模拟——以催化裂化反应-再生系统为例

曹森山，许锋，罗雄麟

（中国石油大学（北京）信息科学与工程学院自动化系，北京 102249）

引 言

在化工过程系统中经常会遇到循环物流，这一般是为了充分利用物料[1-3]，例如，在工艺反应后分离出未反应的或剩余的反应物而进行循环，以提高该反应物流的转化率；在吸收过程中，将吸收液中的溶剂在解吸塔中解吸出来循环回到吸收塔中，实现吸收剂的循环利用，减少了溶剂消耗；在液相或固相催化剂作用下的化学反应中，对反应混合物进行分离后，对液相或固相的催化剂再进行循环；此外，还有共沸精馏中共沸剂的循环利用和萃取过程中萃取剂的回收循环等等。因此，对含有循环物流的化工过程进行流程模拟研究是非常必要的。

在进行化工过程的流程模拟时，从计算策略的角度出发，要解决两个基本问题，一是流程的分隔，二是包含循环物流的子系统迭代。两者相比，如何处理循环物流的迭代是更困难的任务，因为循环物流系统中的循环是贯穿整个系统的，在进行迭代计算时会出现代数环[4-8]问题。代数环让输入输出同时参与迭代计算，引起逻辑混乱，使得迭代计算无法继续进行。目前，化工过程流程模拟方法主要有序贯模块法[9-12]、联立模块法[13-14]和联立方程法[15-16]，其中，应用最广泛的是序贯模块法，该方法的主要思路是将整个系统划分为若干个相互连接的不可分割模块，分别计算每个模块，使其输出相应的信号，最后联立模块计算出最终模型解。对于如何处理化工过程模拟计算中的循环物流，避免流程模拟迭代计算时出现代数环的问题，主要有择优选择断裂物流[17-19]、提供较优初始值、选择合适的加速收敛方法[20-22]、增大迭代次数[23-27]等方法，这对于流程模拟计算的收敛具有重要意义。序贯模块法在对其进行模拟计算时要将所有回路断裂，并给断裂流股变量赋初值，然后依次调用单元模块不断进行迭代计算，直至断裂流股变量的计算结果收敛。在20世纪60年代初，Kou等[28]就提出了断裂的思想。此后，Shacham 等[29]、Motard 等[30]、Naseem 等[31]分别对断裂技术的发展进行了系统的评述，给出了系统的断裂最优准则。随着流程模拟技术的不断发展，目前主流的流股断裂方法主要有Lee-Rudd(L-R)分解法、Upadhye-Grens(U-G)断裂法。而对于断裂后如何迭代计算，常用的迭代方法有直接迭代法、加权迭代法、牛顿迭代法以及拟牛顿法等。汪德友等[32]通过直接法和牛顿法对机坪管网水力流程模拟计算,证明了牛顿迭代法相对于直接迭代法收敛更快。管慧莹[33]提出在以牛顿迭代法为主函数的基础上，引入插值法与其结合，提高了迭代算法的效率。王晓锋[34]提出了两种修正的牛顿迭代法，分析了它们的稳定性和迭代速度，用数值实验证明了两方面都优于经典牛顿迭代法。Fang等[35]通过假定最后三次迭代得到信息关系来近似雅可比矩阵，建立新的拟牛顿迭代方程，具有局部超线性收敛性质，提高了迭代收敛速度。 Bublitz 等[36]基于 Dulmage-Mendelsohn 算法形成块对角雅可比矩阵，降低了拟牛顿迭代方程对子系统初始量的敏感度，确保了整个系统的收敛。这些研究成果都是对牛顿迭代法的有效性进行研究，或是提出新的修正格式，以使收敛结果更加接近真实值。但是这些方法仍然存在一些问题，比如对简单过程系统尚可，对复杂系统则因其规模庞大，很难计算收敛，还有关于循环回路的切断优化判据自相矛盾或者与某些实例不符等等，说明在如何处理含循环物流的流程模拟计算仍需要更深一步的研究。

本文针对含循环物流的化工流程模拟计算难收敛的问题，将稳态模拟计算转化为动态化的差分方程，用离散系统的稳定性理论解决收敛性问题。该方法首先把过程系统的对象模型划分为正向模型和反向模型两类，并且将循环物流中的变量定义为迭代变量和收敛变量，在每个正向模型和反向模型相邻的位置进行断裂，加入增益系数，进而得到收敛递推式。利用控制理论中的稳定性理论来确定迭代方程的增益系数，使得整个系统流程模拟计算达到收敛。根据收敛递推式的形式，采用劳斯判据确定增益系数稳定范围。以催化裂化装置反应-再生系统为例进行稳态流程模拟计算，以反应器作为正向模型，再生器作为反向模型，以再生温度和再生催化剂含碳量作为迭代变量，构造了催化裂化装置反应-再生系统流程模拟的迭代方程，模拟计算结果验证了该方法的可行性和有效性。

1 正向反向模型

化工过程系统是较为复杂和庞大的工业过程，其涉及的建模问题较为复杂，涉及的模型种类也较多。一般情况下需要对真实的化工过程系统单元进行一定的简化处理，使建立的模型能够用现有的数学工具进行求解。由于化工过程系统中的各种单元设备都可以用合适的模型方程进行描述，其模型方程都具有一定的相似性，因此在进行流程模拟计算之前，根据化工过程单元模型的计算形式，可以分为正向模型和反向模型两种。

定义X为单元设备的入口变量；Y为单元设备的出口变量。

正向模型如图1所示。

图1 正向模型Fig.1 Forward model

正向模型单元数学表达式为f(X) =Y，显式表达式，物料流从入口到出口，如实线所示，计算流也是从入口到出口，如虚线所示。这表明正向模型单元在进行流程模拟计算时，其可以从入口一步一步直接计算到出口，不需要附加方法对其处理。管式反应器为典型的正向模型。

反向模型如图2所示。

图2 反向模型Fig.2 Backward model

反向模型单元数学表达式为g(Y) =X，隐式表达式，物料流从入口到出口，如实线所示，但是计算流不能直接从入口计算到出口，计算流如虚线所示，需要对其做断裂，加入附加条件处理后才能进行迭代计算。全混流反应釜为典型的反向模型。

2 循环物流系统

一个化工过程系统可由正向模型与反向模型组合表达。如催化裂化装置的反应-再生系统，提升管是正向模型，计算流与物料流方向相同，从入口计算到出口；而密相床再生器是反向模型，物料流从入口到出口，但是从入口到出口的计算是隐式方程，不能直接从入口计算到出口，需要进行反向计算，计算流与物料流方向相反。正向模型与反向模型组合得到的化工系统如图3所示。

图3 含循环物流的信息流图Fig.3 Information flow diagram with stream circulation

图3 所示的一个简单化工流程中，物料流如实线所示，从入口到出口再到入口组成了一个循环流。但是整个系统的计算流却不能按物料流的顺序进行迭代计算，只有正向模型单元可以按物料流顺序计算，而反向模型单元无法按物料流顺序计算。因为正向模型在模拟计算时因其数学模型是显式表达式，出口变量可有入口变量直接计算得到；反向模型因其数学表达式为隐式表达式，无法从入口到出口直接计算。在模拟计算起始时，应对正向模型与反向模型直接的循环流进行流股断裂，加入增益系数模块，然后再对整个过程系统进行流程模拟计算。

对图3中正向模型和反向模型之间进行流股断裂，如图4所示，其中，X为迭代变量，Y为断裂变量，加入增益系数后有

图4 断裂后的含循环物流的信息流图Fig.4 Information flow diagram of stream circulation system after tearing stream

根据图4，对正向模型和反向模型之间进行流股断裂，加入增益系数后，解决了代数环问题。

化工过程一般为非线性系统，f()和g()为非线性函数，对其基于稳态点进行线性化，得

式(2)为循环物流系统的迭代方程(1)的线性化形式，如果式(2)是稳定的，则式(1)必然收敛，即含循环物流的系统在进行流程模拟计算时会达到计算收敛，所以应当设法通过调整增益系数阵K使式(2)稳定。

3 增益系数的稳定范围

对于迭代方程的线性化方程式(2)，要想满足稳定性，根据雅可比迭代收敛的充分必要条件，(I-KA)的特征值λ必须满足|λ| ≤1。由特征值λ和特征向量X的关系有

式(3)是一个齐次线性方程，存在非零解的充要条件是系数行列式为零，即

式(4)是一个关于λ的高阶方程，与控制理论中离散系统的闭环特征方程形式相同。当方程式(4)的根满足在单位圆内的条件时，λ的取值满足|λ| ≤1。

可以利用线性离散系统稳定性判据中的劳斯判据对式(4)方程的根是否在单位圆内进行判断，进而得出增益系数阵K的元素取值范围。

对增益系数阵K求解的流程如图5所示。

图5 求取收敛模块的流程图Fig.5 The flow chart for obtaining convergence module

线性离散系统的劳斯判据具体要点如下。

表1 劳斯表Table 1 Rolls table

如此继续，每一行的系数由前两行系数组成相应的行列式进行计算。

(4)由劳斯表判断系统稳定性。为保证系统稳定，必须有劳斯表第一列全部大于零，即满足式(6)。

只要式(6)成立，则系统是稳定的，即式(4)的根都在单位圆内，迭代式(1)必将收敛。

对于规模较小的系统(迭代变量X的维数小于等于2)，劳斯表中的元素可表示为增益系数K的函数，利用控制理论中线性离散系统的劳斯判据，根据图5 所示的计算流程可求解增益系数K的稳定范围。当增益系数K在稳定范围内时，则迭代方程式(1)必然收敛，进而含循环物流的系统在进行流程模拟计算时也会达到收敛。对于规模较大的系统，也可以利用控制理论中线性离散系统的劳斯判据，判断增益系数K是否能使迭代方程式(1)收敛，无须进行模拟计算。该方法避免了对反向单元隐式数学模型求解的复杂性，使得不用直接求解隐式数学模型就可以使得整体系统流程模拟迭代达到计算收敛。

4 实例分析

4.1 甲苯氧化反应-分离系统

甲苯氧化制苯甲酸过程是一个简单的反应分离再循环系统。在反应器中进行甲苯氧化制苯甲酸反应后，由于甲苯的注入必须过量，有部分甲苯没反应完，如果不处理这些甲苯，将导致资源的浪费，故将反应器中的物料注入到分离器（精馏塔）中，塔底采出产品苯甲酸，塔顶将未反应完的甲苯和未采出的苯甲酸一同采出再回流到反应器中，以此来不断循环利用甲苯，以节约甲苯资源，提高甲苯利用率。甲苯氧化制苯甲酸反应-分离循环的工艺流程如图6所示。

图6 甲苯氧化制苯甲酸反应-分离系统流程图Fig.6 Schematic diagram of the reaction-separation cycle for the oxidation of toluene to benzoic acid

流股1 为新鲜甲苯进料，流股1 与从分离器出来的循环流股3（主要为未反应的甲苯和未分离出的苯甲酸）混合后进入反应器，从反应器出来的混合流股2 进入分离器后，在塔底得到苯甲酸产品流股4，在塔顶将未反应完的甲苯和未采出的苯甲酸冷凝得到流股3 循环至反应器，以此为一个循环过程，不断重复下去。从图5 中可以看出甲苯氧化制苯甲酸是一个单回路循环过程，其中，反应器数学模型可以归结为反向模型，分离器数学模型可以归结为正向模型。

(1) 分离器数学模型 分离器数学模型的入口变量为反应器出来的混合流股2的苯甲酸摩尔分数x2、混合流股2 摩尔流量q2、产品流股4 摩尔流量q4，假设分离过程为等温过程，对分离器应用物料衡算，可以计算出循环流股3 的苯甲酸摩尔分率x3、产品流股4摩尔流量x4等出口变量，表示为

因此，分离器为正向模型。

(2) 反应器数学模型 反应器数学模型的入口变量为循环流股3的苯甲酸摩尔分率x3、循环流股3摩尔流量q3、新鲜甲苯进料流股1 摩尔流量q1，假设甲苯氧化制苯甲酸过程为等温反应，对反应器应用物料衡算，计算混合流股2 的苯甲酸摩尔分率x2等出口变量。由于反应器近似为全混流CSTR 反应器，在阿伦尼乌斯反应速率方程中反应器出口苯甲酸摩尔分率x2为非线性形式，由入口变量计算出口变量需要迭代计算，而由出口变量计算入口变量则可以直接计算，表示为

因此，反应器为反向模型。

根据图6 的甲苯氧化制苯甲酸的工艺流程，令X=x2,Y=x3，可得反应分离循环系统流程模拟的信息流如图7所示。

图7 反应-分离的信息流图Fig.7 Information flow diagram for reaction-separation

对图7中正向模型和反向模型之间进行流股断裂，X为迭代变量，Y为断裂变量，加入增益系数K后，得到循环物流系统迭代方程

只需要k值满足上述范围,则迭代方程式(10)必然收敛，进而循环物流系统在进行流程模拟时也会达到收敛。

如果按照通常的序贯模块法的计算流程，对式(9)，以x3为自变量，以x2为因变量，通过求解非线性方程进行由x3到x2的正向计算，同样在x3处进行断裂，分别采用直接迭代法、加权迭代法和牛顿法对x3进行迭代计算。

直接迭代法

式中，x(k-1)2与x(k-1)3之间存在隐函数关系x3=g(x2)，需要通过求解非线性方程由x3计算x2。

本文分别采用常规序贯模块法的3种迭代法和本文所述的正向/反向迭代法对甲苯氧化制苯甲酸过程进行流程模拟计算，模拟计算结果如图8所示。

由图8 可以看出，常规序贯模块法的3 种迭代法和本文所述的正向/反向迭代法的计算结果最后都能达到收敛，但常规序贯模块法需要进行求解非线性方程的内部迭代计算，会增加计算时间。此外，本文所述的正向/反向迭代法与直接迭代法相比，其优点是收敛过程稳定，不易出现振荡；与加权迭代法相比，其收敛系数K的稳定范围是可求的，只要K值在稳定范围里，迭代方程式必然收敛，而加权迭代法里的迭代权因子是未知的，需要进行试凑才可能使迭代式收敛；与牛顿迭代法相比，本文方法不需要在线求导，使用更方便，且通过选择合适的收敛系数，收敛速度不亚于牛顿迭代法。

图8 甲苯氧化反应-分离系统四种模拟计算结果对比Fig.8 Comparison simulation of four iterative methods for reaction-separation cycle of the toluene oxidation

4.2 催化裂化装置反应-再生系统

催化裂化装置[37-38]是一个高度非线性、强耦合、分布参数不确定的复杂系统，其反应-再生系统主要由提升管反应器、汽提段、再生器三部分组成，三者之间有着错综复杂、相互关联的流体力学、传热、传质、催化裂化反应动力学等关系。整个催化裂化装置反应-再生系统可以看成一个循环系统，其中催化剂含碳量和温度是贯穿整个系统的循环变量。

以某催化裂化装置为例进行研究，此催化裂化装置为高低并列式提升管催化裂化装置，处理量为600 kt/a，加工减压瓦斯油，单再生器单段再生，添加CO 助燃剂进行完全再生，实现焦炭的完全燃烧，再生效果好，无内外取热系统。其工艺流程如图9 所示，该催化裂化装置过程数据见表2。

表2 催化裂化装置过程数据Table 2 Process data of FCCU

图9 催化裂化装置结构示意图Fig.9 Schematic diagram of FCCU

本文建立单段普通密相床再生的催化裂化装置的数学模型，提升管反应器采用五集总反应动力学模型[39-40]，再生器采用分子筛催化剂烧碳、烧氢反应动力学模型[39]，再生器流化床采用两个虚拟CSTR 串联模型对再生器密相段进行模拟[41]。据此，对催化裂化装置反应-再生系统按提升管反应器、汽提段、再生器密相床分别建立稳态数学模型，对其进行流程模拟计算，该装置是由提升管反应器、汽提段、再生器三部分组成的一个催化剂循环系统。

提升管数学模型和汽提段数学模型可以归结为正向模型，再生器数学模型可以归结为反向模型。

(1)提升管稳态数学模型 取提升管反应器内微元高度分别进行物料和热量衡算，稳态模型数学表达式简化如下:

原料未转化率yA

柴油产率yD

式(24)和式(27)表明提升管反应器和汽提段的数学模型是正向模型，对其进行流程模拟迭代计算时，可以从入口到出口直接计算。

(3) 再生器稳态数学模型 对于再生器，用虚拟CSTR 串联模型进行模拟，假设待生催化剂上焦炭中碳氢组成为恒定，氢碳比为η；再生器催化剂藏量恒定，密相床藏量占总藏量的80%，虚拟CSTR 藏量均为密相床藏量的1/2。

第1 个CSTR 物料衡算和热量衡算的稳态数学模型表达式为：

第2 个CSTR 物料衡算和热量衡算的稳态数学模型表达式为：

含碳量Crg2

式(32)和式(36)说明再生器的两个虚拟CSTR 模型的数学表达式是反向模型，对其进行流程模拟迭代计算时，需要对其进行断裂处理后，才能进行迭代计算。由催化裂化反应-再生系统稳态模型向量方程式(24)、式(27)、式(32)、式(36)可知催化剂含碳量C和温度T是贯穿整个系统的循环变量，使提升管反应器、汽提段、再生器有机联系起来，在进行流程模拟计算时，需要对其进行整体考虑。

再生器的入口条件为：

对于催化裂化反应-再生系统，根据图9的催化裂化装置流程可得反应-再生系统稳态流程模拟如图10所示。

因为提升管数学模型和汽提段数学模型是正向模型，则对其进行流程模拟计算时，可以从入口到出口进行计算。

对于再生器，其数学模型为反向模型，在对其进行流程模拟迭代计算时不能直接计算，需要做断裂处理，因为再生器被分解成两个虚拟的CSTR 都是反向模型，所以在CSTR1 和CSTR2 之间加入一个单位正向模型1。对图10中正向模型和反向模型之间进行流股断裂并加入增益系数后如图11所示。

图10 催化裂化装置反应-再生系统的信息流Fig.10 Information flow diagram of the FCCU reaction-regeneration system

图11 断裂后催化裂化装置反应-再生系统的信息流Fig.11 Information flow diagram of the FCCU reaction-regeneration system after tearing stream

断裂后，对于再生器的两个反向模型的迭代方程有

式中，X2、X4为断裂变量；X3、X5为迭代变量；k为迭代步数；K1、K2为增益系数阵。

对催化裂化反应-再生系统中的迭代方程式(42)～式(45)中的非线性部分F、g1、g2进行线性化，得

即

所有根分布在单位圆内，离散线性系统式(48)稳定，迭代方程将收敛。

根据特征方程式(49)写出劳斯表，应用劳斯判据判断离散线性系统式(48)是否稳定，得到符合稳定性条件的增益系数阵K1,K2。

基于MATLAB 软件对催化裂化反应-再生系统进行流程模拟计算。仿真模拟结果如图12 和图13所示，横坐标代表迭代步数，纵坐标代表迭代变量。图12 表示提升管反应器和汽提段主要变量迭代收敛曲线，图13 表示再生器主要变量迭代收敛曲线。图12 和图13 的迭代曲线表明，在使用本文的方法对催化裂化反应-再生系统进行稳态流程模拟时，主要变量最终都能达到收敛至稳态点。本文模拟结果与生产数据误差如表3所示。

图12 提升管与汽提段主要变量迭代收敛曲线Fig.12 Iterative convergence curves for the main variables in riser and strip section

图13 再生器主要变量迭代收敛曲线Fig.13 Iterative convergence curves for the main variables in regenerator

表3 本文模拟计算的稳态值和生产数据对比Table 3 Comparison of simulated steady-state values and production data from this paper

由表3 可得，总体上本文模拟结果与生产数据基本一致，两种方法下最大的相对误差不超过5%，最小的相对误差为0.02%，证明本文提出的解决含循环物流的流程模拟方法是十分有效的。

5 结 论

基于催化裂化装置反应-再生系统的稳态数学模型，采用流股断裂并且加入增益系数的方法进行稳态流程模拟计算，从模拟结果可以得到以下结论：(1)针对复杂化工过程系统模型，可采用正向/反向模型对系统进行分解处理，变成由正向模型和反向模型组合的流程；(2)对分解后的过程系统模型，对正向模型进行直接计算，对于反向模型，为了解决其迭代时出现的代数环问题，对其与正向模型相邻位置进行流股断裂，得到切断变量和收敛变量，加入增益系数阵，得到迭代方程；(3)由迭代方程的收敛原则，联系到控制理论中的稳定性理论，依据稳定性判据，求解增益系数阵的稳定范围；(4)当增益系数位于稳定范围内时，迭代方程必然收敛，即流程模拟计算达到迭代收敛。

该方法弥补了传统方法难以处理含有多循环回路及多个嵌套循环回路的复杂大型过程系统的不足。对于复杂过程系统，都可以将其看成多个正向模型和反向模型的组合，然后再进行断裂，加入增益系数阵，解决了迭代时出现的代数环问题。然后用控制理论的稳定性理论来求解增益系数阵，找到增益系数的收敛范围并选择合适的增益系数，进而使流程模拟计算达到收敛。