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.

Answer 1

[Harald]

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

[Luca]

Hello Harald,

That's right, since the balances are derived from just one chemical reaction, the resulting equations look similar. I will show quickly my derivation so you could mark errors. But I am pretty sure the equations are setup right.

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.

[Harald]

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]

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?

Answer 2

[Torsten]

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.

[Luca]

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]

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.