Solving ODE system with a constraint/boundary condition

5 views (last 30 days)
Hello,
I am trying to simulate a chemical reaction in a tubular reactor. The reaction is in this form: A -> B + 6 C.
The reaction is a gas phase reaction and the reaction rate is a function of the partial pressures of the three components. I solved the mass balances for steady state (and some other assumtions) and came up with these equations:
I implemented them as an ODE system and want to solve this system via ode45. The code runs without an error but the solution does not fulfill the physical boundary condition, that the sum of all partial pressure has to be the system pressure at all times, or at all positions.
Is there any way to implement a sort of boundary condition, which can be solved by ode45?
See my code below. Please consider that the values for the parameters are arbitrary.
%% Definition of constant parameters
d = 0.028; % tube diameter [m]
length = 0.84; % tube length [m]
A = pi()/4*d^2; % Area [m^2]
V_reactor = A*length; % Reaction volume [m^3]
m_kat = 0.3; % catalyst mass [kg]
rho_packing = m_kat/(V_reactor)*1000; % catalyst density [kg m^-3]
T_shell = 330; % Reaction temperature [°C]
T_shell = T_shell+273.15; % Reaction temperature [K]
p = 4; % reaction pressure [bara]
V_in = 5/60; %inlet stream [m^3 s^-1]
v_z = V_in/A; % velocity [m s^-1]
R = 8.314; % J mol^-1 K^-1
%Starting partial pressures
p_B_in = 1e-4;
p_C_in = 1e-4;
p_A_in = p-p_B_in-p_C_in;
T_in = T_shell;
%Reaction parameters
ReactionRate.k0 = 150;
ReactionRate.EA = 100000;
ReactionRate.a = -0.5;
ReactionRate.m = 1;
ReactionRate.n = -0.5;
%% Calculation of equilibrium constant
eq = @(p,T)(exp((0.06206+p.*0.00176).*((773.15-T)-(443.46876-99.02783.*log(p+5.73037)))))./(1+exp((0.06206+p*0.00176).*((773.15-T)-(443.46876-99.02783.*log(p+5.73037)))));
K_eq_function = @(p,T)((1-eq(p,T))/(eq(p,T)));
K_eq = K_eq_function(p,T_shell);
% Solving ODE system
dz = 0.01; %step size
[z,c] = ode45(@(z, c) ODE_System(z, c, v_z, rho_packing, ReactionRate, R, T_shell, K_eq),0:dz:1,[p_B_in, p_A_in, p_C_in]);
%%Definition of ODE system
function [ddz] = ODE_System(~,c,v_z,rho_packing,ReactionRate,R,T_shell,K_eq)
ddz(1) = ((R*T_shell)/v_z*rho_packing*(((c(3))^ReactionRate.a*ReactionRate.k0*exp(-ReactionRate.EA/ (R*T_shell))*(c(2))^ReactionRate.m*(c(1))^ReactionRate.n*(1-(c(1)/(c(2)*K_eq)))))); % mass balance for B
ddz(2) = (-(R*T_shell)/v_z*rho_packing*(((c(3))^ReactionRate.a*ReactionRate.k0*exp(-ReactionRate.EA/(R*T_shell))*(c(2))^ReactionRate.m*(c(1))^ReactionRate.n*(1-(c(1)/(c(2)*K_eq)))))); % mass balance for A
ddz(3) = 6*((R*T_shell)/v_z*rho_packing*(((c(3))^ReactionRate.a*ReactionRate.k0*exp(-ReactionRate.EA/(R*T_shell))*(c(2))^ReactionRate.m*(c(1))^ReactionRate.n*(1-(c(1)/(c(2)*K_eq)))))); % mass balance for C
ddz = ddz';
end
The condition which has to be fulfilled at any time would be:
Thank you in advance.

Answers (2)

Harald
Harald on 25 Apr 2024
Hi,
if I look at it correctly, your function basically returns
ddz(1) = x;
ddz(2) = -x;
ddz(3) = 6*x;
with some very complicated x.
When summed up, the first two terms cancel each other, and the third one remains. The only way the sum would not change is if x was 0. I suspect that there is some mistake in your ODE function and once that is fixed, the problem will solve itself (at least, up to numerical accuracy).
Best wishes,
Harald
  3 Comments
Harald
Harald on 29 Apr 2024
Hi Luca,
I can't spot an error but that may well be due to my limited knowledge around chemistry.
However if you require to be constant, this implies that . This again would only be true if in my above notation, x was 0 - which is not the case.
Best wishes,
Harald
Harald
Harald on 29 Apr 2024
Reading your comment to Torsten's reply and thinking about it a bit more, I would interpret x as a "reaction rate", i.e. the speed at which A --> B + 6 C happens. Since I suppose the moles will have different masses, the above ODEs then seem to be for number of moles rather than masses?

Sign in to comment.


Torsten
Torsten on 25 Apr 2024
Edited: Torsten on 25 Apr 2024
p_A + p_B should remain constant since dp_A/dz + dp_B/dz = 0 according to your equations.
But this is not the case for p_A + p_B + p_C.
What you write will only hold true for a mass balance, not for a molar balance - there is no such thing as mole conservation.
  2 Comments
Luca
Luca on 29 Apr 2024
Hello Torsten,
I know that there is no molar conservation, in this reaction i make 7 moles out of one mole (A -> B + 6C). What i wrote is the mass balance, just in the partial pressure form. See the derivation below.
So I start with the mass balance in molar form for an ideal PFTR in stationarity. (https://en.wikipedia.org/wiki/Plug_flow_reactor_model). Here is the concentration of component i, u is the flow velocity, and the the source term, which would be the sum of the reaction rate r times the stochiometric factor over all reactions j.
Since there is just one equation to consider, the sum term becomes:
Now I introduce the ideal gas law. . With the definition of the concentration i rewrite the concentration.
No i insert this into the mass balance.
Since the factor RT is constant in this case, i can pull it outside the differential part and put it on the right hand side.
Now since the reaction rate r is the same for all components the only factor which is different in the three equations is the stochiometric factor. This leads to the given equations.
But now there has to be the boundary condition that the sum of all partial pressures has to be equal the system pressure, which i provide. Is there a way to implement that with ode45 or do i have to use other solver. Unfortunately i just worked with ode45.
Torsten
Torsten on 29 Apr 2024
Edited: Torsten on 29 Apr 2024
What i wrote is the mass balance, just in the partial pressure form.
No, these are molar balances written in partial pressure form. If they were mass balances, molar masses for the substances would appear somewhere in your equations.
I don't think that system pressure can be constant if V and T remain constant and the sum of molar amounts of the substances changes. This simply follows from the identity
P = sum_i Ni *RT/V
and the fact that in your case, sum_i Ni is not constant.

Sign in to comment.

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!