필터 지우기
필터 지우기

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일)
Kayla Strathern
Kayla Strathern 2023년 10월 20일
편집: MOHD BELAL HAIDER 2024년 2월 1일
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

이 질문에 답변하려면 로그인하십시오.

답변 (2개)

Walter Roberson
Walter Roberson 2023년 10월 20일
이동: Walter Roberson 2023년 10월 20일
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
이전 댓글 1개 표시
Torsten
Torsten 2023년 10월 20일
편집: Torsten 2023년 10월 20일
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 2023년 10월 20일
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

댓글을 달려면 로그인하십시오.

MOHD BELAL HAIDER
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

카테고리

Help Center 및 File Exchange에서 Interactive Control and Callbacks에 대해 자세히 알아보기

태그

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by