Filter löschen
Filter löschen

Modeling using 6 differential equation and a constraint?

6 Ansichten (letzte 30 Tage)
Albert Escudé Alcover
Albert Escudé Alcover am 27 Feb. 2023
Bearbeitet: Torsten am 27 Feb. 2023
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

Antworten (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov am 27 Feb. 2023
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);

Torsten
Torsten am 27 Feb. 2023
Bearbeitet: Torsten am 27 Feb. 2023
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.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by