Method of lines to solve unsteady state PFR with multiple reactions

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.

[Davide Masiello]

Can you please also share a descriptive summary of the problem, inclusive of the governing equations, initial conditions and boundary conditions?

[Sanjib Das Sharma]

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]

Can you share the code in steady state? It will be very useful.

Answer 1

[Davide Masiello]

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