Azzera filtri
Azzera filtri

Why is my function for a packed bed reactor not returning vectors for each column of F() when I run it like this?

1 visualizzazione (ultimi 30 giorni)
Kayla Strathern
Kayla Strathern il 20 Ott 2023
Modificato: MOHD BELAL HAIDER il 1 Feb 2024
Ran in:
Wvary = linspace(0,100,100);
initialguess = [50 20 0 10 0 1 0.5];
[W, F] = ode23s(@ODE, Wvary, initialguess);
Warning: Failure at t=6.622339e-08. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.352727e-22) at time t.
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

Accedi per rispondere a questa domanda.

Risposte (2)

Walter Roberson
Walter Roberson il 20 Ott 2023
Spostato: Walter Roberson il 20 Ott 2023
Ran in:
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);
Warning: Failure at t=6.622339e-08. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.352727e-22) at time t.
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,:).'
ans = 7×1
49.9903109892204 19.9806219784384 0.00968901078038299 10.0000000000002 8.26565235832224e-13 1.27093572771167e-05 0.499999999999996
ODE(W(end), F(end,:))
ans = 7×1
1.0e+00 * 2.0141469481385 4.02829389627659 -2.01414694813836 -1.17585044769789e-12 1.33361356930335e-13 -3.06779853443709e+16 2.35170089539577e-14
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 Commenti
Mostra 1 commento meno recente
Torsten
Torsten il 20 Ott 2023
Modificato: Torsten il 20 Ott 2023
According to the error message, the solver already stops integration at W = 6.622e-8 s because it has difficulties with your equations. So plotting for Wvary = [0 100] is impossible.
Is W the reactor length ?
Walter Roberson
Walter Roberson il 20 Ott 2023
Ran in:
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,:).'
ans =
-1.878517206251e+27 + 3.86125190294145e+24i -3.75703441250233e+27 + 7.72250380588349e+24i 1.87851720625111e+27 - 3.86125190294165e+24i 108076250255856 - 195486055907.456i 108076250255847 - 195486055907.456i 5.75744702126953e+16 - 85214938520674.1i -2161525005116.43 + 3909721118.14912i
ODE(W(end), F(end,:))
ans =
-1.17420982323434e+26 + 2.02756847268434e+23i -2.34841964646887e+26 + 4.05513694536895e+23i 1.1742098232344e+26 - 2.02756847268443e+23i 5944877735638.84 - 8798974484.64832i 5944877735638.84 - 8798974484.64832i 2.59197022656932e+15 - 2982814166215.26i -118897554712.777 + 175979489.692966i
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

Accedi per commentare.

MOHD BELAL HAIDER
MOHD BELAL HAIDER il 1 Feb 2024
Modificato: MOHD BELAL HAIDER il 1 Feb 2024
check your rate expression. It is not following the balance. specifically the r3

Categorie

Scopri di più su Interactive Control and Callbacks in Help Center e File Exchange

Tag

Prodotti


Release

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by