CN108446480B - 一种适用于一维平推流反应器动态模型求解方法 - Google Patents

Abstract

本发明涉及一种适用于进料为气体的一维平推流反应器动态模型求解方法，从平推流反应器(理想管式反应器)的机理模型出发，采用有限差分法将所描述动态模型的偏微分方程组离散为代数方程组，通过气体的状态方程将模型方程中浓度转化为温度和气体摩尔流量的函数，再联立方程求解得到平推流反应器温度和组成随空间和时间的分布。本发明方法可以用于理想管式反应器进料组分为气相的计算方法，克服以往求解过程中必须将气体体积流量设为恒定参数的假设，使反应器动态模型求解的结果更加准确。

Description

一种适用于一维平推流反应器动态模型求解方法

技术领域

本发明涉及反应器模型技术领域，更具体地说，涉及一种适用于进料为气体的一维平推流反应器动态模型求解方法。

背景技术

按照变量是否和时间有关，可以将数学模型分为稳态模拟和动态模拟。稳态模拟不考虑过程的瞬态变化，因此无法描述装置开停车和异常工况等特性；而动态模拟能较好的体现化工过程的动态特性，可客观详实的描述化工实际过程。

动态模型需考虑时间项的影响，使得平推流反应器的动态数学模型用偏微分方程组描述。此外，气体物料广泛存在于化工反应及分离过程，而气体物料的体积流量容易受到环境参数(如温度、压力等)的变化而变化，该因素也将使得化工过程的动态模型的求解更难。

针对平推流反应器的建模已有相关文献报道，经严密推导得到一维平推流反应器动态模型方程为物料方程式(1)与热量方程式(2)，此模型方程中共有(m+1)个偏微分方程，但模型中包含(2m+1)个变量：分别为各组分浓度C_i、各组分摩尔流量F_i和温度T。由于模型中变量个数比方程个数多m个，故方程组无法求解，还需引入m个方程才能求解。

对于此动态模型的求解，常规方法将反应器中气体的体积流量Q设为恒定参数，其中摩尔流量F_i与浓度C_i的关系可用式(3)表示，此时物料方程式(1)可转化为式(4)，将式(2)和式(4)联立，此时偏微分方程组变量个数与方程个数相同，再用MOL法(也称线上法)将空间离散化，使偏微分方程组转换为常微分方程组，再结合相应的数值算法求解离散化后的常微分方程组即可求解浓度和温度等参数随空间和时间的分布。

F_i＝Q·C_i (3)

其中，C_i为组分i的浓度，mol/m³；F_i为组分i的摩尔流量，mol/s；t为反应时间，s；Q为反应器轴向一点处体积流量，m³/s；V为平推流反应器的体积，m³；r_i为组分i的反应速率，mol/(m³·s)；m为反应体系中物质组分数；U为传热系数，J/(m²·K·s)；a为单位体积的表面积，1/m；T为反应器轴向一点处温度，K；T_a为冷却剂的温度，K；C_pi为组分i的摩尔热容，J/(mol·K)；

为普遍反应速率，mol/(m³·s)；△H_Rx为反应热，J/mol。

但由于实际过程中气体体积流量Q会随温度和气体的摩尔流量、温度等参数变化而变化，此求解结果必然会产生一定误差。

发明内容

本发明的目的在于克服现有技术的不足，使用有限差分法将时间和空间均进行离散化，将整个偏微分方程组变成非线性方程组求解，在使用气体状态方程引入新的方程；由于此过程没有对气体体积流量Q进行假定，其求解结果相较于假设气体体积流量为恒定值的方法更加准确。

由式(1)和式(2)可知，一维平推流反应器的动态模型为一阶双曲型偏微分方程组，此方程组的初始条件和边界条件一般已知，由式(5)和式(6)表示:

t＝0,0<x<V；T＝T',F_i＝F_i',C_i＝C_i' i＝1,2,…m. (5)

V＝0,t＞0；T＝T₀,F_i＝F_i0,C_i＝C_i0 i＝1,2,…m. (6)

其中，x为反应器轴向一点到反应器入口处的体积，m³；T'为反应器初始温度，K；F_i'为组分i的初始摩尔流量，mol/s；C_i'为组分i的初始浓度，mol/m³；T₀为进料温度，K；F_i0为组分i的进料摩尔流量，mol/s；C_i0为组分i的进料浓度，mol/m³。

初始条件和边界条件求在空间和时间的求解域内温度和浓度的分布，先用有限差分法将描述动态模型的偏微分方程组离散为代数方程组，再使用理想气体状态方程将流率和温度转换为浓度的函数，通过求解代数方程求解，得到一维平推流反应器轴向各点不同时刻的温度、浓度、流率等数据。该方法的具体步骤如下：

步骤10，根据反应器大小和计算时长，将空间轴和时间轴进行等分，空间轴分为n等分，时间轴分为N等分。

步骤20，将初始条件和边界条件赋值给求解域上横坐标和纵坐标对应的坐标点，如下面两式：

F_i,j ⁰＝F_i' C_i,j ⁰＝C_i' T_j ⁰＝T' i＝1,2...m；j＝1,2...n (7)

F_i,0 ^k＝F_i C_i,0 ^k＝C_i T₀ ^k＝T i＝1,2...m；k＝1,2...N (8)

其中，F_i,j ⁰、C_i,j ⁰、T_j ⁰为在求解域横坐标上(j,0)点处各组分摩尔流量、浓度、温度；F_i,0 ^k、C_i,0 ^k、T₀ ^k为在求解域纵坐标上(0,k)点处各组分摩尔流量、浓度、温度；

步骤30，使用有限差分法在求解域(j,k)点处将原模型方程转换为代数方程；其中，j和k的初始值均为1；具体如下：

在求解域上(j,k)点处将原模型方程中浓度、摩尔流量和温度随时间或空间的偏导数用相邻两点间的一阶差分代替，模型方程经转换后的式子为：

其中，C_i,j ^k、F_i,j ^k、T_j ^k分别表示在求解域(j,k)点处组分i的浓度、摩尔流量、温度；C_i,j ^k-1、F_i,j ^k-1、T_j ^k-1分别表示在求解域(j,k-1)点处组分i的浓度、摩尔流量、温度；F_i,j-1 ^k、T_j-1 ^k分别表示在求解域(j-1,k)点处组分i的浓度、摩尔流量、温度；

表示在求解域(j,k-1)点处普遍反应速率；r_i,j ^k-1表示在求解域(j,k-1)点处组分i反应速率；△V表示空间步长；△t表示时间步长。

步骤40，将式(9)和式(10)中在(j,k)点处的浓度C_i,j ^k用摩尔流量F_i,j ^k和温度T_j ^k表示：

式中p为反应器总压力，Pa；R为理想气体常数。

步骤50，联立式(9)、式(10)以及式(11)，采用最小二乘法解代数方程组，得到在求解域(j,k)点处浓度C_i,j ^k、摩尔流量F_i,j ^k和温度T_j ^k的数据。

步骤60，判断j是否小于n，如果小于，令j＝j+1，重复步骤30～步骤60；否则，执行步骤70。

步骤70，判断k是否小于N，如果小于，令k＝k+1，重复步骤30～步骤70，否则，结束计算。

本发明的有益效果如下：

一维平推流反应器在求解域内空间轴和时间轴分割个数与实际计算精度有关，由用户自行设定，一般来说当差分步长减小时，计算误差会随之减小，最终结果也会更加精确；本发明一种适用于进料为气体的一维平推流反应器动态模型求解方法能适用于进料状态为气相，解决以往进料为气相时简单将流速设为定值的方法。

以下结合附图及实施例对本发明作进一步详细说明，但本发明的一种适用于进料为气体的一维平推流反应器动态模型求解方法不局限于实施例。

附图说明

图1是本发明的方法流程图；

图2是本发明的求解域分割图；

图3是本发明经典MOL法与本发明方法求解丙酮气相裂解过程动态模拟得到平推流反应器中点处各参数随时间的变化曲线。

图4是本发明经典MOL法与本发明方法求解丙酮气相裂解过程动态模拟得到平推流反应器出口处各参数随时间的变化曲线。

具体实施方式

以下结合附图及实施例对本发明进行进一步的详细说明。

参见图1所示，本发明一种适用于进料为气体的一维平推流反应器动态模型求解方法，包括：

步骤10，参见图2所示，根据反应器大小和计算时长，对求解域进行等分，将空间轴分为n等分，时间轴分为N等分；其中图2(a)为对求解域进行等分的总体分割图，图2(b)为(j,k)点处分割图；

步骤20，将初始条件和边界条件赋值给求解域上横坐标和纵坐标对应的坐标点；

步骤30，使用有限差分法在求解域(j,k)点处将原模型方程转换为代数方程；其中，j和k的初始值均为1；

步骤40，使用状态方程将各组分浓度转换为摩尔流量和温度的函数；

步骤50，联立所有代数方程求解(j,k)点处数据；

步骤60，判断j是否小于n，如果小于，令j＝j+1，重复步骤30～步骤60；否则，执行步骤70；

步骤70，判断k是否小于N，如果小于，令k＝k+1，重复步骤30～步骤70；否则，结束计算。

为验证本发明的有效性，本发明实施例采用丙酮气相裂解生成乙烯酮和甲烷绝热条件下进行分析。其中，反应器压力p为162kPa，进料丙酮浓度C_a0为18.8mol/m³，摩尔流量为38.3mol/s，温度T₀为1035K，原料以体积流量Q₀为2.037m³/s从反应器入口连续流入；初始时刻可假设反应器轴向各点处丙酮浓度C_a'、乙烯酮浓度C_b'、甲烷浓度C_c'均为0mol/m³以及轴向各点处初始温度T'为1035K。该反应的动力学模型如下式所示：

一维平推流反应器动态模型求解方法如下：

步骤a，根据反应器大小为1m³和计算时间为1s，将空间轴和时间轴进行等分，空间轴分为40等分，时间轴分为320等分。

步骤b，将初始条件和边界条件赋值给求解域上横坐标和纵坐标对应的坐标点：

F_a,j ⁰＝F_b,j ⁰＝F_c,j ⁰＝0mol/s j＝1,2...40 C_a,j ⁰＝C_b,j ⁰＝C_c,j ⁰＝0mol/m³ j＝1,2...40

T_j ⁰＝1035K j＝1,2...40

F_a,0 ^k＝38.3mol/s F_b,0 ^k＝F_c,0 ^k＝0mol/s k＝1,2...320

C_a,0 ^k＝18.8mol/m³ C_b,0 ^k＝C_c,0 ^k＝0mol/s k＝1,2...320

T₀ ^k＝1035K k＝1,2...320

步骤c，初始化j＝1,k＝1，求解在求解域上(1,1)点处将原模型方程中浓度、摩尔流量和温度随时间或空间的偏导数用相邻两点间的一阶差分代替，模型方程经转换后的式子为：

步骤d，将求解域在(1,1)点处的浓度C_i,1 ¹用摩尔流量F_i,1 ¹和温度T₁ ¹表示：

步骤e，联立所有代数方程组，使用最小二乘法解代数方程组，得到在求解域(1,1)点处浓度C_i,1 ¹、摩尔流量F_i,1 ¹和温度T₁ ¹的数据。

步骤f，令j＝j+1，重复步骤c至步骤e，即可求当前时刻下一个空间点的浓度、摩尔流量和温度的数据。

步骤g，令k＝k+1，重复步骤c至步骤f，即可求下一时刻所有空间点的浓度、摩尔流量和温度的数据。

为验证本发明方法的有效性，分别给出了将体积流量设为恒定参数的MOL法和本发明方法求解得到平推流反应器中点(V＝0.5m³)和出口(V＝1m³)处各参数的稳态值与文献中的理论值的比较结果如表1和2所示，其中表1表示反应器中点处稳态值与理论值比较，表2表示反应器出口处稳态值与理论值比较。

1)两种方法得到的温度稳态值和理论值相对误差较小(均低于0.2％)。

2)但本发明方法得到浓度稳态值较MOL法得到的浓度稳态值更精确，反应器中点处MOL法求解得到的浓度稳态值与理论值相比的相对误差为6.065％，而本文方法得到模拟结果的相对误差为2.030％；同样反应器出口处MOL法的得到结果的相对误差为8.341％，而本文方法的相对误差为1.073％。

3)由于MOL法将气体体积流量设为恒定参数，在反应器中点和出口处体积流量的稳态值与理论值差距较大，而本文方法的体积流量的稳态值和理论值基本一致。

表1

表2

参见图3和图4所示，分别为经典MOL法与本发明方法求解丙酮气相裂解过程动态模拟得到平推流反应器中点和出口处各参数随时间的变化曲线。其中，图3(a)为反应器中点处温度随时间分布图，图3(b)为反应器中乙烯酮浓度随时间分布图；图4(a)为反应器出口处温度随时间分布图，图4(b)为反应器出口处乙烯酮浓度随时间分布图。由图3及图4可以清晰的发现两种方法动态模拟得到温度和浓度随时间变化趋势基本一致。由于初始时刻反应器轴向各点处丙酮浓度假设为零，只有当原料丙酮从入口流到中点和出口处才能开始反应，为此在中点和出口处反应发生前(如图3中约0.1s前和图4中约0.2s前)温度和浓度均为初始值。反应开始后各项参数均随时间发生变化：由于该体系发生的反应为吸热反应，温度会随时间下降；产物乙烯酮浓度曲线会随时间上升。经过一段时间中点和出口处各参数会达到稳态值(图3中约0.4s后和图4中约0.7s后)。

虽然两种方法求解得到各参数随时间变化的曲线大致相同，但图3和图4显示出不同方法得到结果的动态特性并不相同，本文方法得到轴向各点处的各参数开始变化的时间点早于MOL法的结果，越靠近反应器出口这种现象越明显。此结果就是由气体体积流量变化所引起，根据表1和表2可清晰的发现反应器中点处及出口处气体体积流量的理论值均大于MOL法的假设值，即反应器实际流速比MOL法的假设值要快，故实际过程物料在反应器中某点处的动态响应点应较MOL法所得到的结果要早。而发明方法中所求体积流量与理论值一致，因此本发明所得动态响应点比MOL法的动态响应点更准确。

上述实施例仅是用来说明本发明，而并非用作对本发明的限定，只要是依据本发明的技术实质，对上述实施例进行变化、变型等都将落在本发明的权利要求的范围内。

Claims (1)

1.一种适用于进料为气体的一维平推流反应器动态模型求解方法，其特征在于，包括：

步骤10，根据反应器大小和计算时长，对求解域进行等分，将空间轴分为n等分，时间轴分为N等分；

步骤20，将初始条件和边界条件赋值给求解域上横坐标和纵坐标对应的坐标点；

步骤30，使用有限差分法在求解域(j,k)点处将原模型方程转换为代数方程；其中，j和k的初始值均为1；

步骤40，使用状态方程将各组分浓度转换为摩尔流量和温度的函数；

步骤50，联立所有代数方程求解(j,k)点处数据；

步骤60，判断j是否小于n，如果小于，令j＝j+1，重复步骤30～步骤60；否则，执行步骤70；

步骤70，判断k是否小于N，如果小于，令k＝k+1，重复步骤30～步骤70；否则，结束计算；

所述动态模型通过物料方程式和热量方程式表示，如下：

其中，C_i为组分i的浓度，mol/m³；F_i为组分i的摩尔流量，mol/s；t为反应时间，s；V为平推流反应器的体积，m³；m为反应体系中物质组分数；U为传热系数，J/(m²·K·s)；a为单位体积的表面积，1/m；T为反应器轴向一点处温度，K；T_a为冷却剂的温度，K；C_pi为组分i的摩尔热容，J/(mol·K)；

为普遍反应速率，mol/(m³·s)；△H_Rx为反应热，J/mol；

所述将初始条件和边界条件赋值给求解域上横坐标和纵坐标对应的坐标点，如下表示：

F_i,j ⁰＝F_i' C_i,j ⁰＝C_i' T_j ⁰＝T' i＝1,2...m；j＝1,2...n (7)

F_i,0 ^k＝F_i C_i,0 ^k＝C_i T₀ ^k＝T i＝1,2...m；k＝1,2...N (8)

其中，F_i,j ⁰、C_i,j ⁰、T_j ⁰分别表示在求解域横坐标上(j,0)点处各组分的摩尔流量、浓度、温度；F_i,0 ^k、C_i,0 ^k、T₀ ^k分别表示在求解域纵坐标上(0,k)点处各组分的摩尔流量、浓度、温度；

所述初始条件和边界条件通过如下方式表示：

t＝0,0<x<V；T＝T',F_i＝F_i',C_i＝C_i' i＝1,2,…m (5)

V＝0,t＞0；T＝T₀,F_i＝F_i0,C_i＝C_i0 i＝1,2,…m (6)

其中，x为反应器轴向一点到反应器入口处的体积，m³；T'为反应器初始温度，K；F_i'为组分i的初始摩尔流量，mol/s；C_i'为组分i的初始浓度，mol/m³；T₀为进料温度，K；F_i0为组分i的进料摩尔流量，mol/s；C_i0为组分i的进料浓度，mol/m³；

所述使用有限差分法在求解域(j,k)点处将原模型方程转换为代数方程，包括：

在求解域上(j,k)点处将原模型方程中浓度、摩尔流量和温度随时间或空间的偏导数用相邻两点间的一阶差分代替，模型方程经转换后的式子为：

其中，C_i,j ^k、F_i,j ^k、T_j ^k分别表示在求解域(j,k)点处组分i的浓度、摩尔流量、温度；C_i,j ^k-1、F_i,j ^k-1、T_j ^k-1分别表示在求解域(j,k-1)点处组分i的浓度、摩尔流量、温度；F_i,j-1 ^k、T_j-1 ^k分别表示在求解域(j-1,k)点处组分i的浓度、摩尔流量、温度；

表示在求解域(j,k-1)点处普遍反应速率；r_i,j ^k-1表示在求解域(j,k-1)点处组分i反应速率；△V表示空间步长；△t表示时间步长；

所述使用状态方程将各组分浓度转换为摩尔流量和温度的函数，包括：

将式(9)和式(10)中在(j,k)点处的浓度C_i,j ^k用摩尔流量F_i,j ^k和温度T_j ^k表示：

其中，p为反应器总压力，Pa；R为理想气体常数；

所述联立所有代数方程求解(j,k)点处数据，具体为：

联立式(9)、式(10)以及式(11)，采用最小二乘法解代数方程组，得到在求解域(j,k)点处浓度C_i,j ^k、摩尔流量F_i,j ^k和温度T_j ^k的数据。

Images