Solving ODE system with a constraint/boundary condition
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.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. 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. ode, ode45, boundary condition MATLAB Answers — New Questions
