Why is my function for a packed bed reactor not returning vectors for each column of F() when I run it like this?
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
Walter Roberson
2023-10-20
移动:Walter Roberson
2023-10-20
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 个评论
Walter Roberson
2023-10-20
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
2024-2-1
编辑:MOHD BELAL HAIDER
2024-2-1
check your rate expression. It is not following the balance. specifically the r3
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Interactive Control and Callbacks 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!