查看: 2815|回复: 0

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

[复制链接]

字体大小: 正常放大

浅夏110

542 主题	15 听众	1万积分

TA的每日心情

	开心 2020-11-14 17:15

签到天数: 74 天

[LV.6]常住居民II

群组: 2019美赛冲刺课程

群组: 站长地区赛培训

群组: 2019考研数学桃子老师

群组: 2018教师培训（呼伦贝

群组: 2019考研数学站长系列

电梯直达

1^#

发表于 2020-6-10 10:27 |只看该作者 |倒序浏览

|招呼Ta 关注Ta |邮箱已经成功绑定

例 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)];
%**********************************

————————————————
版权声明：本文为CSDN博主「wamg潇潇」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_29831163/article/details/89711536

zan