Solving ODE system with a constraint/boundary condition
5 views (last 30 days)
Show older comments
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.
0 Comments
Answers (2)
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
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
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?
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
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.
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!