Why is my function for a packed bed reactor not returning vectors for each column of F() when I run it like this?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Wvary = linspace(0,100,100);
initialguess = [50 20 0 10 0 1 0.5];
[W, F] = ode23s(@ODE, Wvary, initialguess);
FA = F(:,1);
FB = F(:,2);
FC = F(:,3);
FD = F(:,4);
FE = F(:,5);
P = F(:,6); %Pressure [Bar]
% T = F(:,8); %Temperature [K]
% Ta = F(:,7); % Coolant temperature [K]
X = F(:,7);
R = 8.314; %[J/mol.K]
plot(W,FA,W,FB,W,FC,W,FD,W,FE)
function d_dW = ODE(W,F)
FA = F(1); %Flowrate of CO
FB = F(2); %Flowrate of H2
FC = F(3); %Flowrate of methanol
FD = F(4); %Flowrate of CO2
FE = F(5); %FLowrate of water
P = F(6); %Pressure [Bar]
X = F(7);
R = 8.314; %[J/mol.K]
% Initial conditions:
P0 = 50; %[Bar]
T0 = 230+273.15; %[K]
FT0 = 200;
T = T0;
% Mole fractions:
FT = FA + FB + FC + FD +FE;
yA = FA./FT;
yB = FB./FT;
yC = FC./FT;
yD = FD./FT;
yE = FE./FT;
K1 = 10.^((5139./T)-12.621); %[bar]
K2 = 10.^((-2073./T)-2.029); %[bar]
K3 = 10.^((3066./T)-10.592); %[bar]
bCO = 2.16*10^-5.*exp(46800./(R.*T)); %[bar^-1]
bCO2 = 7.05*10^-7.*exp(61700/(R.*T)); %[bar^-1]
bH = 6.37*10^-9.*exp(84000/(R.*T)); %[bar^-0.5]
PA = FA./FT.*P;
PB = FB./FT.*P;
PC = FC./FT.*P;
PD = FD./FT.*P;
PE = FE./FT.*P;
k1 = 4.89*10^7.*exp(-11300./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k2 = 9.64*10^11.*exp(-152900./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k3 = 1.09*10^5.*exp(-87500./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
r1 = k1.*bCO.*((PA.*PB.^(2/3) - (PC./(PB^0.5.*K1)))./(1+bCO.*PA+bCO2.*PD).*(PB.^0.5+(bH.*PE)));
r2 = k2.*bCO2.*(PD.*PB - (PA.*PE./K2))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE));
r3 = k3.*bCO2.*((PA.*PB.^(3/2)-(PA.*PE./(PB.^(3/2).*K3)))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE)));
%Pressure effects:
rho_bulk = 2.192*yB(end)+43.474*yC(end)+63.512*yD(end)+24.651*yE(end); %[kg/m3]
viscosity = 1.3495e-05*yB(end)+1.772e-05*yC(end)+2.4572e-05*yD(end)+1.8283e-05*yE(end); %[N/m]
MW = (yA(end)*28.01 + yB(end)*2.016 + yC(end)*32.04 + yD(end)*44.01 + yE(end)*18.02)/1000; %[kg/mol]
dp = 50e-6;
rho_cat = 440; %[kg/m3]
phi = 0.51;
r = 0.5; %[m]
a = pi.*r.^2;
F_final = FA(end)+FB(end)+FC(end)+FD(end)+FE(end);
v0 = (FT0*MW)./rho_bulk;
v = v0.*(FT./FT0).*(P0./P);
u = (v)./ (pi*r^2*(phi));
%Rate expressions:
dFAdW = -r1 + r2;
dFBdW = -2.*r1 - r2 - 3.*r3;
dFCdW = r1 + r3;
dFDdW = -r2 + r3;
dFEdW = r2 + r3;
dPdW = -(1./(rho_cat.*(1-phi).*(a))).*(((150.*viscosity.*(1-phi).^2.*u)./(phi.^3.*dp.^2))+((1.75.*(1-phi).*rho_bulk.*u.^2)./(phi.^3.*dp)));
%Conversion expression:
FD0 = 50; %[mol/s]
dXdW = -dFDdW./FD0;
d_dW = [dFAdW; dFBdW; dFCdW; dFDdW; dFEdW; dPdW; dXdW];
end
0 commentaires
Réponses (2)
Walter Roberson
le 20 Oct 2023
Déplacé(e) : Walter Roberson
le 20 Oct 2023
You need to investigate to figure out why, but dPdW is suddenly going from relatively small to relatively large.
format long g
Wvary = [0 100]; %linspace(0,100,100);
initialguess = [50 20 0 10 0 1 0.5];
[W, F] = ode23s(@ODE, Wvary, initialguess);
FA = F(:,1);
FB = F(:,2);
FC = F(:,3);
FD = F(:,4);
FE = F(:,5);
P = F(:,6); %Pressure [Bar]
% T = F(:,8); %Temperature [K]
% Ta = F(:,7); % Coolant temperature [K]
X = F(:,7);
R = 8.314; %[J/mol.K]
tiledlayout('flow');
nexttile; plot(W,FA); title('FA');
nexttile; plot(W,FB); title('FB');
nexttile; plot(W,FC); title('FC');
nexttile; plot(W,FD); title('FD');
nexttile; plot(W,FE); title('FE');
nexttile; plot(W, F(:,6)); title(':6');
nexttile; plot(W, F(:,7)); title(':7');
F(end,:).'
ODE(W(end), F(end,:))
function d_dW = ODE(W,F)
FA = F(1); %Flowrate of CO
FB = F(2); %Flowrate of H2
FC = F(3); %Flowrate of methanol
FD = F(4); %Flowrate of CO2
FE = F(5); %FLowrate of water
P = F(6); %Pressure [Bar]
X = F(7);
R = 8.314; %[J/mol.K]
% Initial conditions:
P0 = 50; %[Bar]
T0 = 230+273.15; %[K]
FT0 = 200;
T = T0;
% Mole fractions:
FT = FA + FB + FC + FD +FE;
yA = FA./FT;
yB = FB./FT;
yC = FC./FT;
yD = FD./FT;
yE = FE./FT;
K1 = 10.^((5139./T)-12.621); %[bar]
K2 = 10.^((-2073./T)-2.029); %[bar]
K3 = 10.^((3066./T)-10.592); %[bar]
bCO = 2.16*10^-5.*exp(46800./(R.*T)); %[bar^-1]
bCO2 = 7.05*10^-7.*exp(61700/(R.*T)); %[bar^-1]
bH = 6.37*10^-9.*exp(84000/(R.*T)); %[bar^-0.5]
PA = FA./FT.*P;
PB = FB./FT.*P;
PC = FC./FT.*P;
PD = FD./FT.*P;
PE = FE./FT.*P;
k1 = 4.89*10^7.*exp(-11300./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k2 = 9.64*10^11.*exp(-152900./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k3 = 1.09*10^5.*exp(-87500./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
r1 = k1.*bCO.*((PA.*PB.^(2/3) - (PC./(PB^0.5.*K1)))./(1+bCO.*PA+bCO2.*PD).*(PB.^0.5+(bH.*PE)));
r2 = k2.*bCO2.*(PD.*PB - (PA.*PE./K2))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE));
r3 = k3.*bCO2.*((PA.*PB.^(3/2)-(PA.*PE./(PB.^(3/2).*K3)))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE)));
%Pressure effects:
rho_bulk = 2.192*yB(end)+43.474*yC(end)+63.512*yD(end)+24.651*yE(end); %[kg/m3]
viscosity = 1.3495e-05*yB(end)+1.772e-05*yC(end)+2.4572e-05*yD(end)+1.8283e-05*yE(end); %[N/m]
MW = (yA(end)*28.01 + yB(end)*2.016 + yC(end)*32.04 + yD(end)*44.01 + yE(end)*18.02)/1000; %[kg/mol]
dp = 50e-6;
rho_cat = 440; %[kg/m3]
phi = 0.51;
r = 0.5; %[m]
a = pi.*r.^2;
F_final = FA(end)+FB(end)+FC(end)+FD(end)+FE(end);
v0 = (FT0*MW)./rho_bulk;
v = v0.*(FT./FT0).*(P0./P);
u = (v)./ (pi*r^2*(phi));
%Rate expressions:
dFAdW = -r1 + r2;
dFBdW = -2.*r1 - r2 - 3.*r3;
dFCdW = r1 + r3;
dFDdW = -r2 + r3;
dFEdW = r2 + r3;
dPdW = -(1./(rho_cat.*(1-phi).*(a))).*(((150.*viscosity.*(1-phi).^2.*u)./(phi.^3.*dp.^2))+((1.75.*(1-phi).*rho_bulk.*u.^2)./(phi.^3.*dp)));
%Conversion expression:
FD0 = 50; %[mol/s]
dXdW = -dFDdW./FD0;
d_dW = [dFAdW; dFBdW; dFCdW; dFDdW; dFEdW; dPdW; dXdW];
end
3 commentaires
Walter Roberson
le 20 Oct 2023
If you switch to a non-stiff solver then you can see that the values go complex. If you zoom on the plots you can see that they start to get a non-zero imaginary component at pretty much the place that ode23s stops. And that the difficulty appears to correspond closely to the place where the 6th parameter goes negative.
format long g
Wvary = [0 100]; %linspace(0,100,100);
initialguess = [50 20 0 10 0 1 0.5];
[W, F] = ode45(@ODE, Wvary, initialguess);
FA = F(:,1);
FB = F(:,2);
FC = F(:,3);
FD = F(:,4);
FE = F(:,5);
P = F(:,6); %Pressure [Bar]
% T = F(:,8); %Temperature [K]
% Ta = F(:,7); % Coolant temperature [K]
X = F(:,7);
R = 8.314; %[J/mol.K]
S = @(V) [real(V), imag(V)];
tiledlayout('flow');
nexttile; plot(W,S(FA)); title('FA');
nexttile; plot(W,S(FB)); title('FB');
nexttile; plot(W,S(FC)); title('FC');
nexttile; plot(W,S(FD)); title('FD');
nexttile; plot(W,S(FE)); title('FE');
nexttile; plot(W, S(F(:,6))); title(':6');
nexttile; plot(W, S(F(:,7))); title(':7');
F(end,:).'
ODE(W(end), F(end,:))
figure();
tiledlayout('flow');
L = 8e-8;
nexttile; plot(W,S(FA)); title('FA'); xlim([0 L]);
nexttile; plot(W,S(FB)); title('FB'); xlim([0 L]);
nexttile; plot(W,S(FC)); title('FC'); xlim([0 L]);
nexttile; plot(W,S(FD)); title('FD'); xlim([0 L]);
nexttile; plot(W,S(FE)); title('FE'); xlim([0 L]);
nexttile; plot(W, S(F(:,6))); title(':6'); xlim([0 L]);
nexttile; plot(W, S(F(:,7))); title(':7'); xlim([0 L]);
function d_dW = ODE(W,F)
FA = F(1); %Flowrate of CO
FB = F(2); %Flowrate of H2
FC = F(3); %Flowrate of methanol
FD = F(4); %Flowrate of CO2
FE = F(5); %FLowrate of water
P = F(6); %Pressure [Bar]
X = F(7);
R = 8.314; %[J/mol.K]
% Initial conditions:
P0 = 50; %[Bar]
T0 = 230+273.15; %[K]
FT0 = 200;
T = T0;
% Mole fractions:
FT = FA + FB + FC + FD +FE;
yA = FA./FT;
yB = FB./FT;
yC = FC./FT;
yD = FD./FT;
yE = FE./FT;
K1 = 10.^((5139./T)-12.621); %[bar]
K2 = 10.^((-2073./T)-2.029); %[bar]
K3 = 10.^((3066./T)-10.592); %[bar]
bCO = 2.16*10^-5.*exp(46800./(R.*T)); %[bar^-1]
bCO2 = 7.05*10^-7.*exp(61700/(R.*T)); %[bar^-1]
bH = 6.37*10^-9.*exp(84000/(R.*T)); %[bar^-0.5]
PA = FA./FT.*P;
PB = FB./FT.*P;
PC = FC./FT.*P;
PD = FD./FT.*P;
PE = FE./FT.*P;
k1 = 4.89*10^7.*exp(-11300./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k2 = 9.64*10^11.*exp(-152900./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
k3 = 1.09*10^5.*exp(-87500./(R.*T)); %[mol.s^-1.bar^-1.kgcat^-1]
r1 = k1.*bCO.*((PA.*PB.^(2/3) - (PC./(PB^0.5.*K1)))./(1+bCO.*PA+bCO2.*PD).*(PB.^0.5+(bH.*PE)));
r2 = k2.*bCO2.*(PD.*PB - (PA.*PE./K2))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE));
r3 = k3.*bCO2.*((PA.*PB.^(3/2)-(PA.*PE./(PB.^(3/2).*K3)))./((1+bCO.*PA+bCO2.*PD).*((PB.^0.5)+bH.*PE)));
%Pressure effects:
rho_bulk = 2.192*yB(end)+43.474*yC(end)+63.512*yD(end)+24.651*yE(end); %[kg/m3]
viscosity = 1.3495e-05*yB(end)+1.772e-05*yC(end)+2.4572e-05*yD(end)+1.8283e-05*yE(end); %[N/m]
MW = (yA(end)*28.01 + yB(end)*2.016 + yC(end)*32.04 + yD(end)*44.01 + yE(end)*18.02)/1000; %[kg/mol]
dp = 50e-6;
rho_cat = 440; %[kg/m3]
phi = 0.51;
r = 0.5; %[m]
a = pi.*r.^2;
F_final = FA(end)+FB(end)+FC(end)+FD(end)+FE(end);
v0 = (FT0*MW)./rho_bulk;
v = v0.*(FT./FT0).*(P0./P);
u = (v)./ (pi*r^2*(phi));
%Rate expressions:
dFAdW = -r1 + r2;
dFBdW = -2.*r1 - r2 - 3.*r3;
dFCdW = r1 + r3;
dFDdW = -r2 + r3;
dFEdW = r2 + r3;
dPdW = -(1./(rho_cat.*(1-phi).*(a))).*(((150.*viscosity.*(1-phi).^2.*u)./(phi.^3.*dp.^2))+((1.75.*(1-phi).*rho_bulk.*u.^2)./(phi.^3.*dp)));
%Conversion expression:
FD0 = 50; %[mol/s]
dXdW = -dFDdW./FD0;
d_dW = [dFAdW; dFBdW; dFCdW; dFDdW; dFEdW; dPdW; dXdW];
end
MOHD BELAL HAIDER
le 1 Fév 2024
Modifié(e) : MOHD BELAL HAIDER
le 1 Fév 2024
check your rate expression. It is not following the balance. specifically the r3
0 commentaires
Voir également
Catégories
En savoir plus sur Interactive Control and Callbacks dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!