Filter löschen
Filter löschen

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

1 Ansicht (letzte 30 Tage)
Kayla Strathern
Kayla Strathern am 20 Okt. 2023
Bearbeitet: MOHD BELAL HAIDER am 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

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (2)

Walter Roberson
Walter Roberson am 20 Okt. 2023
Verschoben: Walter Roberson am 20 Okt. 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 Kommentare
1 älteren Kommentar anzeigen
Torsten
Torsten am 20 Okt. 2023
Bearbeitet: Torsten am 20 Okt. 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 am 20 Okt. 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

Melden Sie sich an, um zu kommentieren.

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

Kategorien

Mehr zu Interactive Control and Callbacks finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by