FCC六集总模型介绍

1六集总模型的化学说明2
1.1集总的划分 (2)
1.2模块算法流程 (4)
2相关数学运算说明 (5)
2.1计算反应器出口瞬时产率分布(四阶R UNGE-K UTTA法) (5)
2.2计算平均产率分布（G AUSS-L EGENDRE六点法） (6)
2.3估计动力学速率常数（改进的M ARQUARDT法） (7)
3实例说明 (9)
1六集总模型的化学说明
1.1集总的划分
将原料油和产品按照馏程分为6个集总，其中：
产物分为柴油（DS）、汽油（GL）、液化气（LPG）；油浆和原料油部基本为>350℃的馏分，因此将这两部分分为减压渣油（VR）和减压蜡油（VGO）两部分；由于干气和焦炭不是目标产品，因此视为一个集总考虑。

六个集总的馏程如下：
表1 集总馏程划分
集总间的反应网络如
图1所示。

图 1 集总间的化学反应图
假设提升管中气体流动状态假定为等温理想活塞流，质点内扩散以及催化剂的反混忽略不计，集总间的反应均视为一级反应。

由连续性方程和反应速度方程推导模型的基本方程式为：
y -------各集总的摩尔浓度, 集总mol/每小时进油量和进水量t ; X -------油气在提升管内的相对距离，无因次； P --------反应压力，Mpa ； φ(t ) --------催化剂失活函数；
WH S ---------真实重时空速； T------------反应温度, K ； k ------------反应速率常数；
其中催化剂失活函数φ(t )采用基于停留时间的Voorhies 失活函数 : φ(t ) = exp (- a × t )
j
j j
wh j y k y RT S t P X
y ∑-
=)
(d d ϕ
a -------失活速率常数，h-1；
t --------催化剂停留时间，h；
设Yinshu =P*φ(t)* /(
S*R*T)* j y)，
WH
则六个集总的反应动力学方程分别为：
dy(VR)/d(X)= –(k1+k2+k3+k4+k5) *y(VR)*Yinshu;
dy(VRO)/d(X)= (k1*y(VR) -(k6+k7+k8+k9)*y(VGO))*Yinshu;
dy(DS)/d(X)= (k2*y(VR) + k6*y(VGO) -(k10+k11+k12) *y(DS))*Yinshu;
dy(GL)/d(X)= (k2*y(VR) +k7 *y(VGO) +k10*y(DS) –(k13+k14)*y(GL))*Yinshu; dy(LPG)/d(X)= (k4*y(VR) + k8*y(VGO) +k11*y(DS)+k13*y(GL)-k15 *y(LPG)) *Yinshu;
dy(CODG)/d(X)= (k5*y(VR) + k9*y(VGO) +k12*y(DS)+k14*y(GL) +k15
*y(LPG)) *Yinshu;
1.2模块算法流程
要实现预测产物分布的目的，必须先求出个集总之间的反应速率常数k向量，因此，模块分为两个部分：一是根据原料油性质和反应条件求解反应速率常数k 向量，二是根据k向量针对不同性质的原料油计算理论产物分布。

模型第一部分的求解过程为：在一定反应温度下，首先给定反应网络15个速率常数k的初值，用四阶Runge-Kutta（四阶龙格库塔）法求出模型方程组的解，即反应器出口的瞬时产率分布计算值，再用Guass六点积分法计算出反应器出口的时间平均产率分布，然后用计算值与实验值建立目标函数S，根据目标函数采用改进的Marquardt法（非线性最小二乘法）反复迭代后计算反应速率常数k向量。

模型第二部分的求解过程为：输入原料油性质反应条件和第一部分计算出来的k向量，用四阶Runge-Kutta法求出反应器出口的瞬时产率分布计算值，再用Guass六点积分法计算出反应器出口的时间平均产率分布。

具体算法流程如图2和图3所示：
图 2 计算速率常数k 部分算法流程
图 3 产物预测部分算法流程
2 相关数学运算说明
2.1 计算反应器出口瞬时产率分布(四阶Runge-Kutta 法)
集总的反应动力学方程用四阶Runge-Kutta 法计算，其计算结果为反应器出口的瞬时产率分布。

利用已知条件，可构成下面关于一个初值问题的常微分方程，即：
i C WH i KC C RT
S MW
P dX dC )(ϕ=
00y y X ==
其中y 0为实验用原料油的集总浓度分布向量。

由于油气馏分的停留时间比催化剂停留时间要短得多，故可以认为在X ∈[0，1]积分区间内，t c 是个常数。

采用龙格—库塔法，将积分区间X ∈[0，1]分成若干个小区间，从y 0出发逐次确定各节点处所对应得瞬时产率分布，最后得出反应器出口处（X=1）所对应得瞬时产率分布。

设X i =iH ，y(X i )=y(iH)≈y i ，其中H 为计算时所用的步长，则有递推公式：
)222(643211Z Z Z Z H y y i i
++++=+
其中： ),(1y X f Z i =
)2
1,2/(12Z H y H X f Z i i ++=
)2
1,2/(23Z H y H X f Z i i ++=
),(34Z H y H X f Z i i ++=
该方法实质是从初始点出发的一个递推过程，所以总是把初始点作为起点，把要计算函数值的点作为终点。

因此，步长H 的大小直接影响求解的精度。

为了满足一定精度的要求同时又不增加太多的计算次数，可采用变步长计算。

2.2 计算平均产率分布（Gauss-Legendre 六点法）
利用Runge-Kutta 法由模型方程得出的是反应器出口瞬时产率，而实验测定的值以及目标函数中观察变量采用的是时间平均产率。

因此，我们用高斯六点积分公式由瞬时产率计算时间平均产率。

高斯求积方法的基本思想是：在节点数一定时，通过适当选择节点位置，从而尽可能地提高求积的精确度。

Gauss-Legendre 求积公式为：
求积分⎰-1
1)(dx x f ，权ρ(x)=1，x ∈[-1,1]
∑⎰
=-≈n
i i i x f A dx x f 0
1
1
)()(
其中x i （i=0,1,…,n ）是n+1阶Legendre 多项式P n+1(x)的零点，求积系数为：
⎰
-++'-=1
111)()()
(dx x P x x x P A i
n i n i （i=0,1,…,n ）
实际上，x i 和A i 可以计算并由表查出；对于六点积分公式，n=5，其积分节点和求积系数为：
对一般区间[a,b]上得积分⎰b
a dx x f )(通过变量代换
x=
2)(2)(b a t a b ++-，dx=dt a
b 2
- 将x ∈[a,b]变换为t ∈[-1,1]，再用Gauss-Legendre 求积公式计算积分。

对于集总模型来说，x ∈[0,1]，x=2
1
21+t ，其中t c 为催化剂停留时间。

则反应
器出口时间平均产率分布为：
∑=→→
='6
1)(21i i i X y A y
2.3 估计动力学速率常数（改进的Marquardt 法）
改进的Marquardt 法即非线性最小二乘法，
目标函数为：∑=→→→→→'-'-=M j j j j T
j j a y Y Q y Y K S 1
)()()(
（1）
其中： a K
——待求模型参数向量。

M ——实验的组数；
j Y →
——第j 次实验所测得的时间平均产率分布向量；
j y →'——模型观察变量的计算值向量。

j y →'=f(t c ，→y )，→
y 为反应器出口瞬时产
率分布的计算值向量。

→
j Q 为加权矩阵，对f(t c ，→
y )在初值点)0(a K
附近进行一阶泰勒展开：
a K a K K c j K K y y f y t f y f y a a a
∆∂∂∙∂∂+=='→
→
)0()0()0()/()/(),()，tc ( （2）
将（2）式代入（1）式得：
（)()(a a K S K S
∆=）=
a K a K K c j j M
j T
a K a K K c j K K y y f y t f Y Q K K y y f y t f Y a a a a a a ∆∙∂∂∙∂∂--∆∙∂∂∙∂∂--∑=)
0()0()0(1
)0()0()0()()(),([])()(),([[（3） 从式（3）可以看出，目标函数转化为a K ∆的函数，问题变为对a K
∆求极小值问题：令
])()(),([)()(20)()(1)0()0()0()0(=∆∂∂∂∂--∂∂∂∂-⇒=∆∂∆∂∑=M j a K a K c j j T K T K a a a K K y y f y t f Y Q y f K y K K S a a a a
上式整理得：
a
K a
M
j M j K j T K T K a K c j j T K T K a K K y y f Q y f K y y t f Y Q y f K y a a a a a a a
∆∂∂∂∂∂∂∂∂=-∂∂∂∂∑∑==])()()()(2[]),([)()(2)0(11)0()0()0()0()0()0(
（4）
在（4）中，Y 为实验观测值，)
0(),(a K c y t f
的值可以根据已知条件，利用高斯六点法计算。

又因为1])(1[0=∂∂=∂∂⎰y dt t y t y f
c t c
，故只要)0()(a K a K y ∂∂能求出，a K ∆即可用通常解线性方程组的方法求出，使迭代搜索模型参数过程得以继续进行下去。

令)0()0()(a K a
K y
∂∂=Λ，)0(Λ 称为灵敏性系数矩阵。

为求出)0(Λ 矩阵中的各元素，对模型方程),,(a K y t h dX y d
=关于K a 取偏导并变换积分次序，得：
)0()0()0()0()()()()()(a a
a a a K a
K a K K a K a K y dX d K y y h K h dx y d K
∂∂=∂∂∂∂+∂∂=∂∂
即Λ∂∂+∂∂=Λ
)0()0()()(a a K K a
y h K h dX d
解上式关于Λ 的常微分方程，即可得Λ。

令式（4）中：
)0(1)0()0()0()()()()(2a a a a K a M
j K j T K T K a
K y y f Q y f K y
∂∂∂∂∂∂∂∂∑==)]0([a K T
∑=-∂∂∂∂M
j K c j j T K T K a
a a a y t f Y Q y f K y 1)0()0()0(]),([)()(2
=)]0([(a K g
则：)]0([a K T a K ∆=)]0([(a K g
即：a K ∆=)]0([1a K T -∙)]0([(a K g
由)()()1(n K n K n K a a a ∆+=+反复迭代搜索)(a K S
的最小值。

如果函数),(y t f c
与线性函数相差不大，则上述迭代过程容易收敛到参数的最小二乘解。

当),(y t f c
与线性函数差别相当大，且)0(a K 又远离极小点时，若仍忽略泰勒展开式的高阶项，则迭代过程往往发散。

因此，选取一个较好的初值非常重要，但对于集总模型这样的多参数估计来说，要做到这点相当困难。

另外，当线性方程组的系数矩阵呈病态时，矩阵求逆困难或无法实现，因此，可以使用改进的高斯—牛顿法或阻尼最小二乘法，即在线性方程组系数矩阵的对角元上加上一个阻尼因子d ，其迭代公式如下：
)]([))]([()()1(1n K g I d n K T n K n K a a a a
∙++=+-
其中d 在整个迭代过程中不断变化，其原则是保证矩阵的正定性。

在收敛的
情况下，d 取较小的值或0；仅当目标函数)]([)]1([n K S n K S a a
>+时，才选取较大的d 以保证矩阵正定，使目标函数值下降。

3 程序说明
程序是用Matlab 上完成的，改为C#时只需注意rk.m 中使用的ode45（）函数需要按照四阶龙格库塔法重新编译。

具体方法参见2.1。

计算部分程序文件分为4部分：main.m ；guass.m ；rk.m ；fangcheng.m 。

fangcheng.m ：描述集总间的反应动力学方程式；
rk.m：调用fangcheng.m，用四阶龙格库塔法计算反应器出口瞬时产率；
guass.m：调用rk.m，将瞬时产率转化为平均产率；
main.m：调用guass.m，用改进的最小二乘法加速k的收敛效率。

预测部分程序文件与算部分程序文件相同，main中的算法略有不同。

即在
计算部分程序的基础上把k的迭代计算改为一次计算，直接输入k值后计算出产物分布。

需要输入的量：
本程序使用2007年标定数据为例进行试运算，
具体数据如下：
表 2 原料油馏程性质
从原料油馏程中可知，原料油中98%左右的组分馏程为350℃以上，其中500℃馏出率为32%，因此可将原料油分为VR（质量分数68%）和VGO（质量分数32%）两部分，进一步计算VR和VGO的进料摩尔浓度为：
zu(1)=0.6200*q(1)/M1/(q(1)+q(2));
zu(2)=0.3800*q(1)/M2/(q(1)+q(2));
zu(1) ------VR的摩尔浓度，mol/g；
zu(2) ------VGO的摩尔浓度，mol/g；
M1-------- VR的摩尔质量，g / mol；
M1-------- VR的摩尔质量，g / mol；
q(1)-------原料油进料速率，t/h
q(1)-------水蒸气进提升管总速率，t/h
则进料的集总向量为[zu(1), zu(2),0,0,0,0]。

输入反应速率常数k的初值。

jisuanzhi =
0.0058 0.0694 0.1460 0.4061 0.2449 0.1277
shiyanzhi =
0.0038 0.0718 0.1467 0.4108 0.2344 0.1294
wucha =
53.3503 3.3158 0.4699 1.1410 4.4854 1.2922。