Spatio-temporal plug flow reactor simulation with method of lines

I'm trying out and tinkering a few Matlab codes and techniques in my spare time and encountered this particular problem. Note that I'm fairly new to Matlab.

The code and problem statement that I started from is from the below link, which is transient (time-dependent) simulation of a plug flow reactor with a single component.

http://matlab.cheme.cmu.edu/2011/11/17/modeling-a-transient-plug-flow-reactor/

What I want is that another component, say B, is specified as product of the reaction. Hence, it has two components involved with this reaction scheme:

A --> B

with additional rate of B formation as -vo*diff(C)./diff(V)-k*CA(2:end).^2

The problem is that the vectorization allows only representation of C(i) as concentration of A along the volume/length of the reactor, instead of type of reacting species. As result, the type of reacting species cannot be specified in the code as matrix of C. The code I wrote so far and the error message is given as follows.

function PFR_t_L_sim_v2_1

% function method_of_lines
clear all; clc; close all


    function dCdt = reaction(t,C)
        r_1 = k*C(1,2:end).^2;
        % we use vectorized operations to define the odes at each node
        % point.
        dCdt_1 = -vo*diff(C)./diff(V)- r_1;
        dCdt_2 = -vo*diff(C)./diff(V)+ r_1;
        dCdt = [0; dCdt_1; dCdt_2];
       
        dCdt = dCdt(:);
    end

Ca0 = 2; % Entering concentration
Cb0 = 1; % Entering concentration
vo = 2; % volumetric flow rate
volume = 20;  % total volume of reactor, spacetime = 10
k = 1;  % reaction rate constant

N = 100; % number of points to discretize the reactor volume on

init = zeros(N,1); % Concentration in reactor at t = 0
init(1) = Ca0; % concentration component 1 at entrance
init(2) = Cb0; % concentration component 1 at entrance

V = linspace(0,volume,N)'; % discretized volume elements, in column form

tspan = [0 20];
[t,C] = ode15s(@reaction,tspan,init);

plot(t,C(:,end))
xlabel('time')
ylabel('Concentration at exit')
legend('C_A', 'C_B','Location','N')

end

Error message

Error using odearguments (line 95)
PFR_T_L_SIM_V2_1/REACTION returns a vector of length 1, but the length of initial conditions vector is 100. The vector returned by PFR_T_L_SIM_V2_1/REACTION and
the initial conditions vector must have the same number of elements.

Could anyone help me to sort this out?

Thanks in advance.

Answer 1

[Torsten]

function PFR_t_L_sim_v2_1
% function method_of_lines
clear all; clc; close all
    function dCdt = reaction(t,C)
        C1 = C(1:N);
        C2 = C(N+1:2*N);
        dC1dt = zeros(N,1);
        dC2dt = zeros(N,1);
        dC1dt(2:N) = -vo*(C1(2:N)-C1(1:N-1))./(V(2:N)-V(1:N-1)) - k*C1(2:N).^2;
        dC2dt(2:N) = -vo*(C2(2:N)-C2(1:N-1))./(V(2:N)-V(1:N-1)) + k*C1(2:N).^2;
        dCdt = [dC1dt;dC2dt];
        
        %r_1 = k*C(1,2:end).^2;
        % we use vectorized operations to define the odes at each node
        % point.
        %dCdt_1 = -vo*diff(C)./diff(V)- r_1;
        %dCdt_2 = -vo*diff(C)./diff(V)+ r_1;
        %dCdt = [0; dCdt_1; dCdt_2];
        
        %dCdt = dCdt(:);
    end
Ca0 = 2; % Entering concentration
Cb0 = 1; % Entering concentration
vo = 2; % volumetric flow rate
volume = 20;  % total volume of reactor, spacetime = 10
k = 1;  % reaction rate constant
N = 100; % number of points to discretize the reactor volume on
init = zeros(2*N,1); % Concentration in reactor at t = 0
init(1) = Ca0; % concentration component 1 at entrance
init(N+1) = Cb0; % concentration component 1 at entrance
%init = zeros(N,1); % Concentration in reactor at t = 0
%init(1) = Ca0; % concentration component 1 at entrance
%init(2) = Cb0; % concentration component 1 at entrance
V = linspace(0,volume,N)'; % discretized volume elements, in column form
tspan = [0 20];
[t,C] = ode15s(@reaction,tspan,init);
plot(t,[C(:,N),C(:,2*N)])
%plot(t,C(:,end))
xlabel('time')
ylabel('Concentration at exit')
legend('C_A', 'C_B','Location','N')
end

[Anthony Hamzah]

Thanks a lot! I think keeping C_1, C_2, ... C(i) in one column matrix is a nice workaround.

By the way, does anyone have idea to visualize the data as shown below (x and y axes as time and type of species, respectively)? I guess plot3 could be enough but I'm not quite sure how to express the Y-axis length in the "meshgrid".