Method of lines to solve unsteady state PFR with multiple reactions

68 次查看（过去 30 天）
Sanjib Das Sharma
Sanjib Das Sharma 2022-9-28
编辑： Faidra 2024-11-12，13:21
Hi,
I am trying to solve an unsteady state PFR with simultaneous ODE for multiple reactions. The steady and unsteady state model with the same set of equations works fine. But when I am using method of lines to solve for an unsteady state PFR, the results are not comin. Actually, the code does not give any conversion. Also, I am not able to specify the inlet conditions correctly. The code compiles very well though.
  3 个评论
显示 1更早的评论
Sanjib Das Sharma
Sanjib Das Sharma 2022-9-28
Problem statement: Alkane to olefins in a plug flow reactor. The model contains 10 species involved in 15 reactions. The composite rate equation for each species are written as ODEs in the file ladh_optim_pfr_uns.m with the kinetic expressions. ladh_main_optim_uns.m calls the ODEs and gives the output.
I wrote two other codes with the same set of ODEs, one in transient (only time dimension) and the other in steady_state (only space dimension), solving only he composite rate expressions in both. Now I want to solve the following PDE:
dc/dt = -v0dc/dx - rate_expressions. For this I am using method of lines. The code is written based on an example.
The initial condition is no species in the reactor, so C=0; The boundary condition is C1 = 0.34, C2=0.33, C3=0.33. The rest of the species are predicted based on the model equations.
Let me know if you need additional info.
Thank you.
Sanjib
Faidra
Faidra 2024-11-12，13:21
编辑：Faidra 2024-11-12，13:21
Can you share the code in steady state? It will be very useful.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Davide Masiello
Davide Masiello 2022-9-28
Ran in:
I deleted many commented lines for clarity.
clear, clc
k0 = [0.1023 0.0011 0.0129 0.0000001 0.0008 0.000001 0.163 0.5 0.01 0.004 0.00001 0.00001 0.0002 0.0001 0.0001];
E = [20570 60200 25000 181700 135000 22570 135000 31500 31500 145420 34542 34570 256000 256000 256000];
C0 = [0.34 0.33 0.33 0 0 0 0 0 0 0]; % Initial Condition, e.g. total amount of moles
N = 100;
Check especially the following four lines, which is the right way to formulate your initial conditions.
init = zeros(N*10,1);
init(1) = 0.34;
init(N+1) = 0.33;
init(2*N+1) = 0.33;
volume = 20; % Total volume of the reactor
V = linspace(0,volume,N)';
time = [0,3000]; % seconds
[t,C] = ode15s(@(t,c)ladh_optim_pfr_uns(t,c,k0,E,V,N),time,init);
%% Plot of concentrations vs volume coordinate (V) at different times
time_idx = round(linspace(1,length(t),5));
species = {'Propane','Butane','Ethane','Propylene','Ethylene','Butylene','Methane','Hydrogen','Hexane','Coke'};
for idx = 1:10
figure
plot(V,C(time_idx,N*(idx-1)+1:idx*N))
xlabel('Volume');
ylabel('Concentration');
title(species{idx})
legend([num2str(t(time_idx)),[' s';' s';' s';' s';' s']],'Location','best')
end
function dcdt = ladh_optim_pfr_uns(~,c,k0,E,V,N)
c1 = c(1:N);
c2 = c(N+1:2*N);
c3 = c(2*N+1:3*N);
c4 = c(3*N+1:4*N);
c5 = c(4*N+1:5*N);
c6 = c(5*N+1:6*N);
c7 = c(6*N+1:7*N);
c8 = c(7*N+1:8*N);
c9 = c(8*N+1:9*N);
c10= c(9*N+1:10*N);
R = 8.314;
T = 580 + 273; % Reaction temperature, oK
A = 0.05; % Area of Reactor
V0 = 0.50; % Volumetric Flow Rate of Reacting Gas
Tref = 25 + 273; % Reference temperature, oK
% Initial guess Pre-exponential factors
k_0d = 0.1;
Ed = 15100;
% Defining the k values (reaction constants).
K1 = 8.49*10^8*exp(-118707/R/T); % Equilibrium constant for Propane Dehydrogenation
K2 = 8.49*10^8*exp(-118707/R/T); % Equilibrium constant for Butane Dehydrogenation
K3 = 8.49*10^8*exp(-118707/R/T); % Equilibrium constant for Ethane Dehydrogenation
% Initial rate constants
for i = 1
for j = 1:15
k(i,j) = k0(i,j)*exp(-E(i,j)/R/T);
end
end
% Deactivation Parameter - Exponential Model
kd = k_0d*exp(-Ed/R/T); % deactivation rate constant
a = exp(-kd/R/T); % deactivation rate
% Defining the net reaction rates % Net coke
dc1dt = zeros(N,1);
dc2dt = zeros(N,1);
dc3dt = zeros(N,1);
dc4dt = zeros(N,1);
dc5dt = zeros(N,1);
dc6dt = zeros(N,1);
dc7dt = zeros(N,1);
dc8dt = zeros(N,1);
dc9dt = zeros(N,1);
dc10dt = zeros(N,1);
dc1dt(2:N) = -V0*(c1(2:N)-c1(1:N-1))./(V(2:N)-V(1:N-1))-0.0005*c(2:N); %+(a*(-k(1,1).*(c1(2:N)- (c1(2:N)).*c8(2:N)./K1))-k(1,2).*c1(2:N)-k(1,3).*c1(2:N) - k(1,15).*c1(2:N)); % Propane Consumption
dc2dt(2:N) = -V0*(c2(2:N)-c2(1:N-1))./(V(2:N)-V(1:N-1))+(a*(-k(1,8).*c2(2:N)-k(1,9).*c2(2:N)-k(1,10).*c2(2:N))); % Butane Consumption
dc3dt(2:N) = -V0*(c3(2:N)-c3(1:N-1))./(V(2:N)-V(1:N-1))+(a*(-k(1,6).*c3(2:N) - k(1,7).*c3(2:N).*c8(2:N)+k(1,2).*c1(2:N).*c8(2:N)+k(1,10).*c2(2:N)+k(1,12).*c5(2:N).*c8(2:N)-k(1,14).*c3(2:N))); % Net Ethane Consumption
dc4dt(2:N) = -V0*(c4(2:N)-c4(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,1).*(c1(2:N)- (c1(2:N).*c8(2:N)./K1))- 2*k(1,4).*c4(2:N)- k(1,5).*c4(2:N).*c8(2:N)+k(1,9).*c2(2:N))); % Net Propylene equation
dc5dt(2:N) = -V0*(c5(2:N)-c5(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,3).*c1(2:N)+ k(1,5).*c4(2:N).*c8(2:N)+ k(1,11).*c6(2:N)+ k(1,6).*c3(2:N) - k(1,12).*c5(2:N))); % Net Ethylene equation
dc6dt(2:N) = -V0*(c6(2:N)-c6(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,8).*c2(2:N) - k(1,11).*c6(2:N))); % Net Butylene
dc7dt(2:N) = -V0*(c7(2:N)-c7(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,3).*c1(2:N) + k(1,2).*c1(2:N).*c8(2:N) + k(1,7).*c3(2:N).*c8(2:N)+ k(1,9).*c2(2:N) -k(1,13).*c7(2:N))); % Net Methane % Net Carbon
dc8dt(2:N) = -V0*(c8(2:N)-c8(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,1).*(c1(2:N)- (c1(2:N).*c8(2:N)/K1)) - k(1,2)*c1(2:N).*c8(2:N)+ k(1,6).*c3(2:N) +k(1,8).*c2(2:N)+k(1,13).*c7(2:N)-k(1,5).*c4(2:N).*c8(2:N)-k(1,7).*c2(2:N).*c8(2:N)-k(1,13).*c5(2:N).*c8(2:N)-k(1,4).*c4(2:N).^2.*c8(2:N))); % Net Hydrogen
dc9dt(2:N) = -V0*(c9(2:N)-c9(1:N-1))./(V(2:N)-V(1:N-1))+(a*(k(1,4).*c4(2:N).^2.*c8(2:N))); % Net Hexane
dc10dt(2:N)= -V0*(c10(2:N)-c10(1:N-1))./(V(2:N)-V(1:N-1))+(k(1,13).*c7(2:N)+ k(1,15).*c1(2:N) + k(1,14).*c3(2:N)); % Net coke
dcdt = [dc1dt;dc2dt;dc3dt;dc4dt;dc5dt;dc6dt;dc7dt;dc8dt;dc9dt;dc10dt];
end

类别

在 Help Center 和 File Exchange 中查找有关 Thermal Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by