连续流动釜式反应器动态模拟分析

1.问题提出

连续流动釜式反应器（CSTR）在相同的操作条件下，可能会出现一个或多个不同的操作温度，这被称为反应器的多态现象，在数学上也称为多解。实际出现在哪一个操作温度取决于其起始状态。

现有一个夹套冷却的CSTR反应器，如图1-1所示，在其中进行一级放热反应：A→P + ΔH。需计算反应器达到稳定时的状态，并模拟开车时反应器的动态行为。

‍图1-1 非等温夹套冷却釜式反应器

已知：F = 10-8 m3/s，F为反应进料流量；CA0=5.0 kmol/m3，CA0为进料反应物初始浓度；T0= 300 K，T0为进料温度，T = 250 K，T为出料温度，TJ = 305K，TJ为冷却水温度，V = 2.0×106 m3，V 为反应釜的体积；ρ = 1000 kg/m3，ρ为冷却水的密度；cp = 4.187kJ/(kg·K)，cp为冷却水的比热容；UA = 5.68×106 kJ/(K·s)，U为反应器壁总传热系数，A为传热面积；H = -4.19×104 kJ/kmol，H为一级反应的热焓；k0 = 8.03×1012l/s，k0为一级反应的反应速率常数；E = 9.42×104 kJ/kmol，E为一级反应的活化能；R=8.314 kJ/(kmol·K)，R为气体常数。

2
计算模型

一级反应的动力学方程表达式如下，

对A组分做物料衡算，由积累量=进料量-出料量-反应量，可得到质量衡算方程：

定义空时τ为反应体积与进料体积流量的比值，故上式可写作：

对反应体系做热量衡算，可得到如下热量衡算方程

即 

式中，Q_g和Q_r分别为反应的放热量和移走的热量。

质量衡算和热量衡算两个常微分方程组成了CSTR的动态模型，当浓度C_A和温度T不再变化时即为反应器的稳定状态：

对于稳定状态的计算，实质是解非线性方程组，首先绘制出Q_g(Q_r)~T的关系图，估计出解的大致位置，再代入所估计的初值，由fsolve()函数计算非线性方程组的根。

对于反应器的动态模拟，实质是一个ODE方程组的初值问题，首先先定义ODE方程组，再用预估-校正方法求解在tspan=[0 5000]上浓度和温度随时间变化的动态数据。并用循环语句在一定范围内改变进料温度和浓度的初值，用以计算不同初值下CSTR的动态数据，由此可绘制出连续流动釜式反应器的相平面图。

3
流程框图和程序

反应器稳态计算的流程框图如图3-1所示，开车过程动态模拟的流程框图如图3-2所示。

‍图3-1 稳态计算流程框图            图3-2 动态模拟流程框图

Matlab程序如下：

functionCSTRReactor

global E R k0 FCA0 V T0 UA TJ rho Cp Hr tau

STOPTIME =5000;

F =1.0e-8;               % 液体体积流量, m3/s

CA0 =5.0;             % 进料浓度, kmol/m3

T0 =300;               % 进料温度, K  

TJ =305;               % 夹套温度, K

V =2.0e-6;              % 反应器体积,m3

rho =1000.0;             % 液体密度, kg/m3

Cp =4.187;             % 液体热容, kJ/kg K

UA =5.68e-6;           % 传热系数与传热面积的乘积, kJ/K s

Hr =-4.19e4;             % 反应热, kJ/kmol 

k0 =8.03E+12;          % 指前因子, 1/s

E =9.42e+4;            % 活化能, kJ/kmol

R =8.314;              % 气体常数, kJ/kmol K

tau = V/F;

% 绘制Qg (Qr)~T的关系图以观察稳态点

T =linspace(280 , 380);

k =RateConstant(T);

CA =CA0./(1+tau*k);

[Qg,Qr] =Heat(T,k,CA);

figure

plot(T,Qg,'b-',T,Qr,'k-.');

xlabel('T (K)');

ylabel('Q_g,Q_r (kJ/s)');

legend('Q_g','Q_r')

% 计算稳态温度

[x,f,h] =fsolve(@steadyfun,[5,300]);

CAS =x(1,1)

Ts =x(1,2)

% 动态模拟

tspan =[0  STOPTIME];

[t,y] =naeuler(@dynamicmodel,tspan,[0,T0],1);

CA =y(:,1);

T = y(:,2);

figure

plot(t,T,'b-')

xlabel('t'),ylabel('T(K)')

figure

plot(t,CA,'b-')

xlabel('t'),ylabel('C_A(mol/L)')

% 绘制反应器相平衡图

figure

hold on

TI =linspace(300,340,10);

CAI =linspace(0,CA0,2);

fori=1:length(TI)

    for j=1:length(CAI)

        [t,y] = ode45(@dynamicmodel,tspan,[CAI(j),TI(i)]);

        CA = y(:,1);

        T = y(:,2);

        plot(T,CA)

    end

end

xlabel('T (K)');

ylabel('C_A(mol/L)');

hold off

% 速率常数计算函数---------------------------------------------

function k =RateConstant(T)

global k0 E R

k = k0 *exp(-E/R./T);

% 移热量与放热量计算函数---------------------------------------

function [Qg,Qr] =Heat(T,k,CA)

global F T0 V UATJ rho Cp Hr

Qg =k.*CA*V*(-Hr);                   % 产生的热量,kJ/s

Qr =F*rho*Cp*(T-T0) + UA*(T-TJ);    % 移走的热量,kJ/s

% 计算稳态状态下的温度函数-------------------------------------

functionf=steadyfun(x)

global CA0 tauk0 F rho Cp T0 UA TJ Hr E R V

k =RateConstant(x(2));

f(1)=CA0-x(1) - k.*x(1)*tau;

f(2)=F*rho*Cp*(x(2)-T0)+UA*(x(2)-TJ)-k.*x(1)*V*(-Hr);

%ODE function-------------------------------------------------

function f =dynamicmodel(t,y)

global CA0 tauk0 F rho Cp T0 UA TJ Hr E R V

CA = y(1);T = y(2);

k =RateConstant(T);

f(1) =(CA0-CA)/tau - k*CA;

[Qg,Qr] =Heat(T,k,CA);

f(2) = (Qg- Qr)/(V*rho*Cp);

f = f(:);

% 预估-校正模型解常微分方程组-----------------------------------

function [x,y] =naeuler(dyfun,xspan,y0,h)

x=xspan(1):h:xspan(2);

y=zeros(length(y0),length(x));

y(:,1)=y0(:);

for n =1:length(x)-1

    K1 = dyfun(x(n),y(:,n));

    y(:,n+1) = y(:,n) + h * K1;

    K2 = dyfun(x(n+1),y(:,n+1));

    y(:,n+1) = y(:,n) + h*(K1+K2)/2;

end

x=x';y=y';

4
图像与计算结果

图4-1所示为Q_g(Q_r)~T的关系图，由图可看出所讨论的反应器其Q_g与Q_r曲线只有一个交点，即不存在多解，其唯一解出现在300 K附近，以300 K为解非线性方程组的初值进行迭代计算。经计算稳定浓度C_AS=4.4638mol/L，稳定温度T_s=305.1373K。

‍图4-1 Q_g(Q_r)~T关系图

图4-2所示为反应器开车过程温度随时间变化图。由图可看出，开车过程中温度T在快速上升后，在1500 s左右稳定于T_s，即305.1373 K。

‍图4-2反应器开车过程温度随时间变化图

图4-3所示为反应器开车过程A的浓度随时间变化图。由图可看出，开车过程中浓度C_A在快速上升后又小幅下降，在1500 s左右稳定于C_As，即4.4638 mol/L。

图4-3反应器开车过程A组分浓度随时间变化图

改变进料温度T₀和浓度C_A0的初值，计算不同初值下CSTR开车的动态数据，由此绘制出连续流动釜式反应器的相平面图如图4-4所示，从图中可看出，该反应器在不同初始浓度和进料温度下开车，最后都会达到唯一的稳定状态，各曲线汇聚的点即为该平衡稳定点（T_s, C_As）。

‍图4-4 连续流动釜式反应器的相平面图

5
结论

由模拟结果可知，该反应器的反应转化率较低，但稳定性较好，进一步计算可对其进行改造并确定最佳操作参数。在求解ODE-IVP问题时，本文采用了预估-校正模型，也可如同在计算相平衡图时那样，直接调用Matlab的ode45()函数，由四五阶Runge-Kutta法求解，计算稳定可靠。

非等温反应器中常出现多定态现象，出现一个或者三个稳定点。根据系统和过程中的动力学、热力学规律以及传递现象建立数学模型，并用数值方法求解，可对多定态体系进行模拟计算和过程分析，进而根据模拟结果指导反应器的开停车，为其安全操作提供重要指导。

参考文献：

[1]刘世超,侯言超,刘宏超,朱建华.水合法制丙二醇的CSTR开车模拟及反应器稳定性分析[J].化学反应工程与工艺,2012,28(05):469-474. 

（注：关注微信公众号，并在后台回复“CSTR动态模拟”即可获得MATALB源文件及本文的参考文献~）

本期编辑：邵青楠

本期排版：伍青林