催化裂化反应系统过程模拟优化模型及其求解方法

摘要

本发明涉及催化裂化反应系统模拟优化模型及其求解方法，新模型基于Gupta等的模型，为便于其在过程模拟软件如AspenPlus

说明书

技术领域

本发明涉及石化过程模拟及优化技术，尤其涉及催化裂化反应系统模拟优化模型及其求解方法。

背景技术

作为我国目前处理量最大的原油二次加工手段的催化裂化在石油炼制工业中占有极其重要的地位，催化裂化装置设计及操作的好坏直接影响到炼厂的整体效益。在催化裂化装置的设计及操作阶段，借助催化裂化反应系统过程模拟优化模型，对其反应系统的设计及操作进行必要的模拟优化，能够提高整个装置的设计及操作水平，进而增加炼油厂的经济效益。

催化裂化反应系统过程模拟及优化的基础是需要一个准确实用的工艺计算模型。国内外已有工艺计算模型在催化裂化反应系统过程设计及操作中的应用，如中国石化集团洛阳石油化工工程公司的催化装置模拟优化软件FCCLK^TM，美国ASPENTECH公司的催化装置模拟优化软件Aspen FCC^TM等。这些工艺计算模型及软件包的成功应用对于催化裂化装置的设计及操作优化具有重要的推动意义，但是存在以下缺点：

1这些工艺计算模型是针对具体的催化裂化原料油(基础原料油)而开发的，对于与基础原料油性质差距较大的原料油，会出现计算模拟结果准确度下降的情况。

2由于具体催化裂化装置反应系统千差万别，因此为保证计算模拟的准确度，这些工艺计算模型在开发过程中均设置大量的装置因数，造成其具体应用的不便。

3对于计算模拟结果，这些工艺计算模型能够提供压力、温度、产品收率的变化等基本信息，但是不能提供空隙率、气固滑移因子等对过程设计及操作具有重要参考意义的信息。

4在催化裂化反应系统涉及物种的选择上，这些工艺计算模型均采用集总反应模型，即将原料及产品人为划分为多个集总组分，并将这些集总组分流量及组成的变化视为原料及产品流量及组成的变化。但是，目前过程模拟及优化软件，涉及物种一般为实组分或虚拟组分等。因此，在催化裂化反应系统模拟及优化模型中采用集总反应模型，会造成其难以与过程模拟及优化软件集成，进而难以进行全系统模拟优化的问题。

发明内容

本发明的技术方案如下：

一种催化裂化反应系统模拟优化模型，基于Gupta等[Gupta et al.A new generic approach for the modeling of fluid catalytic cracking(FCC)riser reactor.Chemical Engineering Science，2007，62：4510-4528.]提出的催化裂化反应系统工艺计算模型，为利于耦合过程模拟软件提供的物性计算子程序对模型计算所涉及的物性进行计算，针对物料衡算方程、能量衡算方程及气固流化状态模型中部分表达式提出新的形式，因此模型可作为过程模拟软件的用户单元模型，从而在过程模拟软件中完成催化裂化装置的模拟及优化工作。

以下对本发明作进一步的说明，具体如下：

1物料衡算方程数学表达形式，具体实现如下：

对裂化反应非参与组分i，

f_i，j＝f_i，j-1

对气相裂化反应参与组分i，

$f_{i,j}=f_{i,j-1}+M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{m,i,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_m\ge M_i+M_n,M_i\ge M_n}$

$+M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{m,n,i,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_m\ge M_n+M_i,M_n\ge M_i}-M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{i,m,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_i\ge M_m+M_n,M_m\ge M_n}$

对固相组分焦炭，

$f_{\mathrm{coke},j}=f_{\mathrm{coke},j-1}+{\left(\underset{i=1}{\overset{N}{\Sigma }}\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{i,m,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}(M_i-M_m-M_n)\right)}_{M_i\ge M_m+M_n,M_m\ge M_n}$

该方程原形式参看文献-Gupta et al.Chemical Engineering Science，2007，62：4510-4528.-中公式(9)及公式(11)：

对裂化反应非参与组分i，

f_i，j＝f_i，j-1

对气相裂化反应参与组分i，

对固相组分焦炭，

其中i、m及n为组分的代号；M为参与反应组分的相对分子质量，kg/kmol；j-1代表反应微元进口；j代表反应微元出口；cat代表裂化催化剂；coke代表焦炭；f为质量流量，kg/s；N为组分总数；为反应速率，kmol/(kg(cat)s)；τ代表裂化催化剂在反应微元内停留时间，s。

2对于能量衡算方程，具体实现如下：

能量衡算方程的数学表达一般形式1，

$0=(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j-1}\right)h_{j-1}+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j-1}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})$

$-(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j}\right)h_j+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})-\mathrm{\alpha \pi d\Delta l}(T_{j-1}-T_{\mathrm{env}})$

能量衡算方程的数学表达的一般形式2，

$0=(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j-1}\right)h_{j-1}+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j-1}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})$

$-(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j}\right)h_j+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})-\mathrm{\alpha \pi d\Delta l}(T_j-T_{\mathrm{env}})$

如若假设反应系统为绝热反应系统，则上两式变为，

$0=(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j-1}\right)h_{j-1}+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j-1}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})-(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j}\right)h_j+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})$

其中，

h_j-1＝h_j-1(T_j-1，p_j-1，z_j-1)

h_j＝h_j(T_j，p_j，z_j)

该方程原形式参看文献-Gupta et al.Chemical Engineering Science，2007，62：4510-4528.-中公式(16)：

$(f_{\mathrm{cat}}{c}_{p,\mathrm{cat}}+f_{\mathrm{coke},j-1}{c}_{p,\mathrm{coke}}+\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j-1}{c}_{p,i})\left(T_{j-1}{T}_j\right)=\underset{i}{\Sigma }\underset{m}{\Sigma }\underset{n}{\Sigma }\left({\bar{r}}_{i,m,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}{\mathrm{\Delta H}}_{i,m,n,j}\right)$

α为反应器壁面总传热系数，kW/(m²K)；π为圆周率；d为反应系统直径，m；Δl为反应微元高度，m；T_env为环境温度，K；T_ref为组分热力学状态计算基准温度，K；h为比焓，kJ/kg；c_p为等压比热，kJ/(kg K)；T为裂化反应温度，K；ΔH为裂化反应热，kJ/kmol。h为局部温度、压力及组成的函数，该变量的计算调用过程模拟软件提供的物性计算子程序。

3对于气固流化状态计算模型，具体实现如下，

对于所涉及的气相速度计算方程、固相速度计算方程及气相粘度计算方程(Gupta et al 模型假设气相粘度为常数)依据如下公式：

$u_g=\frac{\underset{i=1}{\overset{N}{\Sigma }}{f}_i}{{\mathrm{A\rho }}_g{\delta }_g}$

ρ_g＝ρ_g(T，p，z)

$u_{c1}=\frac{f_{\mathrm{cat}}+f_{\mathrm{coke}}}{{\mathrm{A\rho }}_{c1}(1-{\delta }_g)}$

ρ_c1＝(ρ_cat(1-y_coke)+ρ_cokey_coke)(1-δ_g)

μ_g＝μ_g(T，p，z)

固相速度计算方程、气相速度计算方程原形式参看文献-Gupta et al.Chemical Engineering Science，2007，62：4510-4528.-中公式(26)及(38)：

$u_{c1}\approx \frac{f_{\mathrm{cat}}}{{\mathrm{A\rho }}_{\mathrm{cat}}(1-{\delta }_g)}$

$u_g=\frac{\underset{i=1}{\overset{N}{\Sigma }}{f}_i}{{\mathrm{A\rho }}_g{\delta }_g}$

${\rho }_g=\frac{P\underset{i=1}{\overset{N}{\Sigma }}\frac{f_i}{M_i}}{\mathrm{RT}}$

其中c1代表催化剂簇团；g代表气相；u为速度，m/s；ρ为密度，kg/m³；μ为粘度，Pa·s；A为提升管局部截面积，m²；δ为局部体积分数；z为局部气相摩尔组成；y为局部质量分数；R为通用气体常数，8.3145J/(mol K)；P为局部压力，Pa。ρ_g与μ_g为局部温度、压力及组成的函数，该变量计算需调用过程模拟软件提供的物性计算子程序。

4催化裂化反应系统模拟优化模型的求解方法，具体实现方法：将催化裂化反应系统划分为大量层状反应微元，根据反应微元进口的气相流股或者固相流股组成、温度及压力计算得到反应微元出口气相流股或者固相流股的组成、温度及压力，同时调用过程模拟软件提供的物性计算子程序计算模型求解所需的组分或者流股物性。

利用计算机高级语言FORTRAN编写AspenPlus^TM用户模型子程序USER2，然后利用AspenPlus^TM软件提供的编译链接环境AspenPlusSimulationEngine^TM对其进行编译链接，生成.DLL文件，然后在AspenPlus^TM中调用该.DLL文件实现本模型与AspenPlus^TM软件的完全集成，从而可以利用AspenPlus^TM软件虚拟组分表征系统、组分物性参数数据库、物性计算系统以及物性计算模型库；对于气相比焓、气相密度等组分或者流股物性的计算通过AspenPlus^TM软件提供的物性计算子程序的调用来实现。

本发明的有益效果是：新模型采用了实组分及虚拟组分反应动力学，便于调用过程模拟软件提供的组分或者流股物性计算子程序，因此本发明催化裂化反应系统模拟优化模型及求解方法有效避免了上文背景技术中所提到的其他模型及求解程序的一般缺点。

附图说明

图1为示出催化裂化反应系统模拟优化模型涉及的实组分或虚拟组分裂化反应网络的描述例示意图。

图2为示出求解催化裂化反应系统模拟及优化模型采用的一般有限元法的示意图。

图3为示出催化裂化反应系统出口反应终产物实组分及虚拟组分组成计算结果的描述例示意图。

图4为示出反应系统装置温度变化趋势计算结果的描述例示意图。

图5为示出反应系统装置压力变化趋势计算结果的描述例示意图。

图6为示出反应系统装置空隙率变化趋势计算结果的描述例示意图。

图7为示出反应系统装置气固速率变化趋势计算结果的描述例示意图。

具体实施方式

结合附图对本发明做进一步的详细说明：

本发明催化裂化反应系统模拟优化模型的求解环境为AspenPlus^TM过程模拟软件，若其他软件提供所需用户模型编程接口亦可以采用。

图1为实组分或虚拟组分裂化反应网络，假设反应系统发生的反应均为裂化反应，且一单位的组分i裂化反应生成一单位的组分m及一单位的组分n。求解过程中涉及的裂化反应动力学模型采用Gupta等[Gupta et al.A new generic approach for the modeling of fluid catalytic cracking(FCC)riser reactor.Chemical Engineering Science，2007，62：4510-4528.]提供的关联模型。

对于物料衡算方程，依据如下公式计算：

对裂化反应非参与组分i，

f_i，j＝f_i，j-1

对气相裂化反应参与组分i，

$f_{i,j}=f_{i,j-1}+M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{m,i,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_m\ge M_i+M_n,M_i\ge M_n}$

$+M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{m,n,i,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_m\ge M_n+M_i,M_n\ge M_i}-M_i{\left(\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{i,m,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}\right)}_{M_i\ge M_m+M_n,M_m\ge M_n}$

对固相组分焦炭，

$f_{\mathrm{coke},j}=f_{\mathrm{coke},j-1}+{\left(\underset{i=1}{\overset{N}{\Sigma }}\underset{m=1}{\overset{N}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{\bar{r}}_{i,m,n,j}{f}_{\mathrm{cat},j}{\tau }_{\mathrm{cat},j}(M_i-M_m-M_n)\right)}_{M_i\ge M_m+M_n,M_m\ge M_n}$

其中i、m及n为组分的代号；M为参与反应组分的相对分子质量，kg/kmol；j-1代表反应微元进口；j代表反应微元出口；cat代表裂化催化剂；coke代表焦炭；f为质量流量，kg/s；N为组分总数；为反应速率，kmol/(kg(cat)s)；τ代表裂化催化剂在反应微元内停留时间，s。

对于能量衡算方程，因被模拟装置工程设计中考虑了有效的保温，故假设反应器绝热(即α＝0)，采用了如下公式计算：

$0=(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j-1}\right)h_{j-1}+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j-1}\underset{T_{\mathrm{ref}}}{\overset{T_{j-1}}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})-(\left(\underset{i=1}{\overset{N}{\Sigma }}{f}_{i,j}\right)h_j+f_{\mathrm{cat}}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{cat}}\mathrm{dT}+f_{\mathrm{coke},j}\underset{T_{\mathrm{ref}}}{\overset{T_j}{\int }}{c}_{p,\mathrm{coke}}\mathrm{dT})$

h_j-1＝h_j-1(T_j-1，p_j-1，z_j-1)

h_j＝h_j(T_j，p_j，z_j)

其中h为比焓，kJ/kg；c_p为等压比热，kJ/(kg K)；T为裂化反应温度，K；T_ref为组分热力学状态计算基准温度，K。h为局部温度、压力及组成的函数，该变量的计算调用AspenPlus^TM过程模拟软件提供的物性计算子程序。

对于流化状态计算模型所涉及气相速度、固相速度及气相粘度的计算，依据如下公式：

$u_g=\frac{\underset{i=1}{\overset{N}{\Sigma }}{f}_i}{{\mathrm{A\rho }}_g{\delta }_g}$

ρ_g＝ρ_g(T，p，z)

$u_{c1}=\frac{f_{\mathrm{cat}}+f_{\mathrm{coke}}}{{\mathrm{A\rho }}_{c1}(1-{\delta }_g)}$

ρ_c1＝(ρ_cat(1-y_coke)+ρ_cokey_coke)(1-δ_g)

μ_g＝μ_g(T，p，z)

其中c1代表催化剂簇团；g代表气相；u为速度，m/s；ρ为密度，kg/m³；μ为粘度，Pa·s；A为提升管局部截面积，m²；δ为局部体积分数；z为局部气相摩尔组成；y为局部质量分数。ρ_g与μ_g为局部温度、压力及组成的函数，该变量计算需调用AspenPlus^TM过程模拟软件提供的物性计算子程序。

图2为本模型求解采用的一般有限元法(见附图2)，其基本原理为将催化裂化反应系统划分为大量层状反应微元，根据反应微元进口的气相流股或者固相流股组成、温度及压力计算得到反应微元出口气相流股或者固相流股的组成、温度及压力，同时调用AspenPlus^TM软件提供的物性计算子程序计算组分或者流股物性。

利用计算机高级语言FORTRAN编写AspenPlus^TM用户模型子程序USER2，然后利用AspenPlus^TM软件提供的编译链接环境AspenPlusSimulationEngine^TM对其进行编译链接(具体方法请参考软件帮助说明)，生成.DLL文件，然后在AspenPlus^TM中调用该.DLL文件实现本模型与AspenPlus^TM软件的完全集成，从而可以利用AspenPlus^TM软件虚拟组分表征系统、组分物性参数数据库、物性计算系统以及物性计算模型库。对于气相比焓、气相密度等组分或者流股物性的计算通过AspenPlus^TM软件提供的物性计算子程序的调用来实现，具体调用程序及接口规定参考该软件帮助说明。

表1为催化裂化反应系统模拟优化模型输入原料的性质，表2为模拟优化模型的输入操作条件。利用表1及表2的输入数据，采用AspenPlus^TM软件对模型进行求解，图3-图7为本催化裂化反应系统模拟优化模型在表1及表2条件下的模拟计算结果。

图3为出反应系统终产物组分累计液相体积，由于本模型的求解在AspenPlus^TM软件完成，可针对工业数据使用AspenPlus^TM软件的参数估计功能对模型参数进行估计，保证其模拟计算结果与工业数据一致，且能够根据估算得到的模型参数对其进行适度外延，显示本模型在预测产品收率上的灵活性。图4为出反应系统装置的流股温度变化趋势计算结果，可以看出，由于吸热裂化反应的存在，反应系统温度迅速下降并维持在500℃左右，其计算结果与工业观察数据(即反应系统出口气固相温度500℃)一致，显示了本模型对于反应系统轴向温度预测的准确度。图5为反应系统轴向压力降变化示意图，可以看出，压力降大约为32kPa，该结果与工业观察数据30.5kPa基本一致。图6为反应系统轴向空隙率的变化趋势图，可以看出，反应系统轴向空隙率迅速上升，并维持在0.95m³/m³左右，该值与流态化文献结论一般值0.8～0.98m³/m³一致。图7为反应系统轴向气固相速率的变化趋势图，可以看出，气固两相速率均持续上升，气相速率最终维持在8m/s，固相速率维持在2m/s，气固滑移因子维持在4左右，这与流态化相关文献给出的气固滑移因子大约为4.3的结论是吻合的，进一步证明了本模型对于气固相速率计算的准确度。

表1催化裂化反应系统模拟优化模型的输入原料性质

表2催化裂化反应系统模拟优化模型基础输入条件

本发明公开催化裂化反应系统模拟优化模型及求解方法，本领域技术人员可通过借鉴本文内容，适当改变模型构成、求解方法等环节实现。本发明已通过具体实施实例进行了描述，相关技术人员明显能在不脱离本发明内容、精神和范围内对本文所述的系统进行改动或适当变更与组合，来实现本发明技术。所有相类似的替换和改动对本领域技术人员来说是显而易见的，他们都被视为包括在本发明精神、范围和内容中。