Why does my graph of molar flowrates in a PFR look so weird

Hi, trying to plot a graph of molar flowrates against volume for the 3 reacting species in a PFR using the ode45 function.

However when it runs, it doesn't look at all like what you'd expect, Fb starts decreasing right away and Fa never reduces to 0.

Code used is:

%% CME3035
clc
Vspan = [0 1];
Y0 = [5 0 0 500.15];
[V, Y] = ode23(@Question1, Vspan, Y0);
subplot(2,1,1)
plot(V,Y(:,1),V,Y(:,2),V,Y(:,3));

legend('Fa','Fb','Fc');
ylabel('Molar flowrate, mol/min');
xlabel('Volume,dm3');
subplot(2,1,2)
plot(V,Y(:,4));
legend('Temperature (K)');
ylabel('Temperature (K)');
xlabel('Volume,dm3');

and in another file the actual function with the rates of reactions

function f = Question1(~,Y)
%series reactions Y is the concentrations of species and V PFR volume
Fa= Y(1);
Fb= Y(2);
Fc= Y(3);
T= Y(4);

%initial conditions
To=500.15; %K
Ta=523.15; %K
Cpa=20; % J/molK
Cpb=80; % J/molK
Cpc=100; % J/molK
Cpi=Cpa; % J/molK
DCp1=Cpb-(2*Cpa);
DCp2=Cpc-(2*Cpb)-Cpa;
dH1=-25000+(DCp1*(T - To));  
dH2=35000+(DCp2*(T - To));


Fio=1; % inert mol/min, remains constant

Ua=150; %J/dm3.min.K
CTo=0.3996; % mol/dm3
FT=Fa+Fb+Fc;
Ca=CTo*((Fa/FT)*(To/T));
Cb=CTo*((Fb/FT)*(To/T));
Cc=CTo*((Fc/FT)*(To/T));
Ci=CTo*((Fio/FT)*(To/T));

%defining rates
R=8.314462; %j / mol k
k1=50*exp((8000/R)*(1/315-1/T)); % dm3/mol.min
k2=400*exp((4000/R)*(1/310-1/T)); % dm3/mol.min
kc=k1/k2; % dm3/mol.min

ra=(k1*((Ca^2)-((1/kc)*Cb)))+(k2*Ca*(Cb^2));
rb=1/2*((k1*((Ca^2)-((1/kc)*Cb)))+(2*k2*(Cb^2)*Ca));
rc=k2*Ca*(Cb^2);

%differential equations 
dFadV=-ra;
dFbdV=-rb;
dFcdV=rc;
dTdV=(Ua*(Ta-T)-((ra-rc)*dH1)-(2*rc*(dH2)))/(Cpa*Fa+Cpb*Fb+Cpc*Fc+Cpi*Fio);
f=[dFadV; dFbdV; dFcdV; dTdV];
end

The above sets of code result in this graph

and this is the question i'm trying to answer

Any help would be much appreciated

Answer 1

[Kartik]

Hi,

From the code you provided, it seems that you are using the `ode23` function to solve a system of ordinary differential equations (ODEs) representing the molar flow rates and temperature in a plug flow reactor (PFR). However, there are a few issues in your code that might be causing unexpected results in the plotted graphs.

1. Initial conditions: The initial condition for species B (`Fb`) is set to 0, which might be causing the decrease in `Fb` right away. If you expect `Fb` to have a non-zero initial value, you should adjust the `Y0` vector accordingly.

2. Temperature calculation: In the `dTdV` equation, the term `2*rc*(dH2)` seems to be incorrect. It should be `rc*(dH2)`. Additionally, make sure that the signs and units of the terms in the equation are correct. Double-check the heat of reaction values (`dH1` and `dH2`) and the heat capacity values (`Cpa`, `Cpb`, `Cpc`, and `Cpi`) to ensure they are consistent.

3. Reaction rates: The calculation of the reaction rates `ra`, `rb`, and `rc` seems correct based on the given reaction equations. However, check that the stoichiometric coefficients and rate constants are accurate for the specific reaction you are modeling.

4. Volume range: The `Vspan` variable is set to `[0 1]`, representing the volume range for the PFR. Ensure that this range is appropriate for your system and that the simulation covers the desired volume range.

Make sure to review and verify these aspects of your code to ensure that the equations, initial conditions, and parameters are correctly implemented. Additionally, consider using the `ode45` function instead of `ode23` for better accuracy, especially if the solution involves steep changes or stiff ODEs.