Modeling using 6 differential equation and a constraint?

Hi there

I've got 6 differential equations + some constraints for the answers.

I have tried using the code that you will find at the bottom, but it doesn't work. I don't know how to calculate the solutions following the constraint given by FT and v that are found at the end of the code, any idea?

clc
clear all
close

% Initial feed
FC3H8_0 = 74167.2; % Propane initial feed [mol/s]
FH2O_0 = FC3H8_0*1.59; % Water initial feed [mol/s]
FCH4_0 = 0;
FH2_0 = 0;
FPropylene_0 = 0;
FC2H4_0 = 0; 
FT0=FC3H8_0+FH2O_0+FCH4_0+FH2_0+FPropylene_0+FC2H4_0;

Y0=[FC3H8_0;FH2_0;FCH4_0;FPropylene_0;FC2H4_0;FH2O_0]; %Initial conditions
Vol = [100:100:100000];% Volume
[Vol,Y] = ode45(@f,Vol,Y0);
plot(Vol,Y(:,1),'-o',Vol, Y(:,2),'-+',Vol,Y(:,3),'-x',Vol,Y(:,4),'-s',Vol,Y(:,5),'-p' )
xlabel('Volum (L)')
ylabel('Cabal (mol/s)')
legend('','fH','fBe','fM','fBi')
fprintf('V=%f  %f   %f   %f   %f   %f %f \n', [Vol,Y]')


% PFR with 2 parallel reactions
function dydVol=f(Vol,Y)
R1=8.314;
R2 = 0.082;
T=1100;                      % [K]
P=5*10^5/101325;        % in atm, bar to Pa to atm

A1=4.69*10^10;                 % [L/mol·s]
Ea1=2.12*10^5;              % [J/mol]
k1=A1*exp(-Ea1/(R1*T));      % [L/mol·s]  

A2= 5.89*10^10;                % [L/mol·s]
Ea2= 2.14*10^5;                % [J/mol]
k2=A2*exp(-Ea2/(R1*T));      % [L/mol·s]    


FC3H8_0 = 74167.2; % Propane initial feed [mol/s]
FH20_0 = FC3H8_0*1.59; % Water initial feed [mol/s]
FCH4_0 = 0;
FH2_0 = 0;
FPropylene_0 = 0;
FC2H4_0 = 0; 


FC3H8=Y(1); FH2=Y(2); FCH4=Y(3); FPropylene=Y(4); Fethylene=Y(5); FH20=Y(6);
dydVol = zeros(6,1);
FT = [];
v = [];
    dydVol(1)=-k1*(FC3H8/v)-k2*(FC3H8/v); 
    dydVol(2)= k2*(FC3H8/v);
    dydVol(3)= k1*(FC3H8/v);
    dydVol(4)= k2*(FC3H8/v);
    dydVol(5)= k1*(FC3H8/v);
    dydVol(6)=0;
    FT = sum(dydVol);
    v=(FT*0.082*T/P);

end

Answer 1

[Sulaymon Eshkabilov]

The problem is with this [ ] assigning step in your code:

v = [ ];

What is the purpose of diving a number by an empty variable v?

dydVol(1)=-k1*(FC3H8/v)-k2*(FC3H8/v);

dydVol(2)= k2*(FC3H8/v);

dydVol(3)= k1*(FC3H8/v);

dydVol(4)= k2*(FC3H8/v);

dydVol(5)= k1*(FC3H8/v);

Answer 2

[Torsten]

v =
FT*0.082*T/P =
sum(dydVol)*0.082*T/P =
1/v*sum(dydVol*v)*0.082*T/P 
->
v^2 = sum(dydVol*v)*0.082*T/P

Thus define the dydVol without dividing by v, sum its elements, multiply by 0.082*T/P, take the positive or negative square root and then divide each element of dydVol by this v-value.