偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

偏微分方程的数值解系列博文：

偏微分方程的数值解(一):定解问题 & 差分解法

偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

偏微分方程的数值解(五)： 二维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法

---

例 4 触煤反应装置内温度及转换率的分布

以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：

 其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为

此外，式中之相关数据及操作条件如下：

(i)反应速率式

其中 P 表示分压(atm)，而速率参数为

上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。

(ii)操作条件及物性数据

题意解析：

将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。

MATLAB 程序设计 将原方程改写成如式（35）的标准式

     因此

根据以上的分析，可编写 MATLAB 程序求解此 PDE 问题，其参考程序如下：

function ex60_3_1
%******************************
% 触媒反应器内温度及转化率的分布
%******************************
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
%******************************
% 给定数据
%******************************
Pt=1.25; %总压(atm)
rw=0.025; %管径(m)
Tw=100+273; %壁温(℃)
G=631; %质量流率(kg/m2hr)
M=30;
y0=0.0323;
Mav=4.47;
rho_B=1200;
Cp=1.74;
dHr=-49250;
h0=65.8;
T0=125+273;
Lw=1;
u=8.03;
R=1.987;
ke=0.65;
hw=112;
De=0.755;
%******************** 
m=1;
%********************
% 取点
%********************
r=linspace(0,rw,10);
L=linspace(0,Lw,10);
%***********************
% 利用 pdepe 求解
%***********************
sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L);
T=sol(:,:,1); %温度
f=sol(:,:,2); %反应率
%***********************
% 绘图输出
%***********************
figure(1)
surf(L,r,T'-273)
title('temp')
xlabel('L')
ylabel('r')
zlabel('temp (0C)')
%
figure(2)
surf(L,r,f')
title('reaction rate')
xlabel('L')
%初始条件函数
%**********************************
function u0=ex20_3_1ic(x)
u0=[125+273 0]';
%**********************************
% 边界条件档
%**********************************
function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L)
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
pl=[0 0]';
ql=[1 1]';
pr=[hw*(ur(1)-Tw) 0]';
qr=[G*Cp 1]'; 
ylabel('r')
zlabel('reaction rate')
%*************************************************
% PDE 函数
%*************************************************
function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr)
global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De
T=u1(1);
f=u1(2);
%
k=exp(-12100/(R*T)+32.3/R);
Kh=exp(15500/(R*T)-31.9/R);
Kb=exp(11200/(R*T)-23.1/R);
Kc=exp(8900/(R*T)-19.4/R);
%
a=1+M-3*f;
ph=Pt*(M-3*f)/a;
pb=Pt*(1-f)/a;
pc=Pt*f/a;
%
rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4;
%
c1=[1 1]';
f1=[ke/(G*Cp) De/u]'.*DuDr;
%s1=[ke/(G*Cp*r)*DuDr(1)-rA*rho_B*dHr/(G*Cp)-2*h0*(T-Tw)/(rw)
s1=[-rA*rho_B*dHr/(G*Cp);rA*rho_B*Mav/(G*y0)];
%********************************** 

 

---

偏微分方程的数值解系列博文：

偏微分方程的数值解(一):定解问题 & 差分解法

偏微分方程的数值解(二): 一维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布

偏微分方程的数值解(四): 化工应用————扩散系统之浓度分布

偏微分方程的数值解(五)： 二维状态空间的偏微分方程的 MATLAB 解法

偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法

---