# Modeling using 6 differential equation and a constraint?

6 Ansichten (letzte 30 Tage)
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
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (2)

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);
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

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.
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### 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