偏微分方程的数值解(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布
例 4 触煤反应装置内温度及转换率的分布以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应 系统之质量平衡及热平衡方程式如下：
https://img-blog.csdnimg.cn/20190430160312927.png
其中T 为温度(℃)， f 为反应率，L 为轴向距离，r 为径向距离。此系统的边界条件为
https://img-blog.csdnimg.cn/20190430160354234.png
此外，式中之相关数据及操作条件如下：
(i)反应速率式
https://img-blog.csdnimg.cn/20190430160411559.png
其中 P 表示分压(atm)，而速率参数为
https://img-blog.csdnimg.cn/20190430160454234.png
上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。R 为理想气体常数(1.987cal/mol·K)。
(ii)操作条件及物性数据
https://img-blog.csdnimg.cn/20190430160511294.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
https://img-blog.csdnimg.cn/20190430160530272.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
题意解析：
https://img-blog.csdnimg.cn/20190430160602531.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
将上式，连同反应速率式，带入平衡方程式中，配合边界条件，可利用 pdepe 求解。
MATLAB 程序设计 将原方程改写成如式（35）的标准式
https://img-blog.csdnimg.cn/20190430160648225.png
因此
https://img-blog.csdnimg.cn/2019043016112275.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzI5ODMxMTYz,size_16,color_FFFFFF,t_70
根据以上的分析，可编写 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=';
%**********************************
% 边界条件档
%**********************************
function =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=';
ql=';
pr=';
qr=';
ylabel('r')
zlabel('reaction rate')
%*************************************************
% PDE 函数
%*************************************************
function =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=';
f1='.*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)];
%**********************************
————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536
页:
[1]