all ODE function returns NaN

2 Ansichten (letzte 30 Tage)
Jonathan Hanandhito
Jonathan Hanandhito am 4 Mär. 2023
Hi, I have a project where I need to model a reaction in a pack bed. The codes are below. It was running well before I added the differential equation for temperature change inside the unit (shown as dydz(7)). When I ran the code, it kept on giving NaN as solutions to the differential equation. Does this mean my function is too stiff for the ODE solver, and I need to split the code in a different function? Thanks.
zspan=0.1:0.1:10;
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0';
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000;
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K
A= pi*(Dtube^2/4); %m^2
%%physical setting
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h
%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube);
[z,y]=ode15s(fun, zspan, IC);
Warning: Failure at t=1.000016e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.
figure
plot(z,y(:,1),z,y(:,2),z,y(:,3),z,y(:,4),z,y(:,5))
figure
plot(z,y(:,6));
figure
plot(z,y(:,7));
function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4;
Dax= 0;
t=0.002;%m
R=8.314;% kJ/kmol.K
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K
%%kinetics setting
ka_ref= 0.499;
kb_ref=6.62*10^-11;
kc_ref= 3453.38;
kd_ref= 1.07;
ke_ref= 1.22*10^10;
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;
%%equilibrium setting;
K1eq_A= 3066;
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029;
%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B);
K2eq= 10^((K2eq_A./y(7)) - K2eq_B);
%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7);
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number
%%thermal properties
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o));
Pr= mu.*cpg./kg; %prandtl number
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma));
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2);
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat;
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15);
Hr2= -39892 + 8.*(y(7)-498.15);
%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o);
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp);
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));
dydz=dydz';
end
  3 Kommentare
Sam Chak
Sam Chak am 5 Mär. 2023
First thing to check is to look for any division-by-zero terms that involve y(7), the newly added state variable. Multiples are found in code. The prime suspect is found.
Second, reduce the time step in the solver so that you see how y(7) evolves over time before the singularity occurs. Sulaymon has shown that y(7), the temperature drops so fast to nearly 0 K in a fraction of 0.1 sec.
Third, the rapid temperature drop must be caused by dT/dt or dydz(7). Some coefficients must have very large values. Are they logical for any physical processes on Earth?
Jonathan Hanandhito
Jonathan Hanandhito am 5 Mär. 2023
Hi Sam and Torsten, I've managed to fix it now based on what you all have said. Some expressions were left out, and the units for some parameters and equations were wrong. Thank you so much for pointing these out!

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov am 4 Mär. 2023
There are two ways to address this issue. (1) Let ode15s to select appropriate step size; (2) make the step size much smaller than it was specified. Also, note that the solutions are within the range of [0.1 ... 0.10002]. Here is the option 1:
zspan=[0.1, 0.10002]; % Alt. way: zspan=[0.1:1e-13:.1002];
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0';
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000;
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K
A= pi*(Dtube^2/4); %m^2
%%physical setting
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h
%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube);
[z,y]=ode15s(fun, zspan, IC);
Warning: Failure at t=1.000016e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.
figure
plot(z,y(:,1),z,y(:,2),z,y(:,3),z,y(:,4),z,y(:,5))
figure
plot(z,y(:,6));
figure
plot(z,y(:,7));
function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4;
Dax= 0;
t=0.002;%m
R=8.314;% kJ/kmol.K
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K
%%kinetics setting
ka_ref= 0.499;
kb_ref=6.62*10^-11;
kc_ref= 3453.38;
kd_ref= 1.07;
ke_ref= 1.22*10^10;
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;
%%equilibrium setting;
K1eq_A= 3066;
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029;
%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B);
K2eq= 10^((K2eq_A./y(7)) - K2eq_B);
%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7);
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number
%%thermal properties
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o));
Pr= mu.*cpg./kg; %prandtl number
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma));
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2);
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat;
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15);
Hr2= -39892 + 8.*(y(7)-498.15);
%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o);
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp);
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));
dydz=dydz';
end

Sulaymon Eshkabilov
Sulaymon Eshkabilov am 4 Mär. 2023
Note that the found numerical solutions of y(:,1), .. y(:,5) are in different scales and therefore, it is appropriate to polt them in separate of using yyaxis option, e.g.:
zspan=[0.1,.1002];
%zspan=[0.1:1e-13:.1002];
n_co0= 39.757; %kmol/h
n_co20= 66.6; %kmol/h
n_h20= 200; %kmol/h
n_meoh0= 0; %kmol/h
n_h2o0= 0; %kmol/h
n0=[n_co0; n_co20; n_h20;n_meoh0;n_h2o0]; %kmol/h
nt0=sum(n0,'all'); %kmol/h
y0=n0./nt0;
y0=y0';
Pin= 50; %bar
Tin= 200+273.15; %K
IC= [y0 Pin Tin];
Ntube= 8000;
Dtube= 0.05; %m
R=8.314; %kJ/kmol.K
A= pi*(Dtube^2/4); %m^2
%%physical setting
mw_co = 28.01; %kg/kmol
mw_co2 = 44.01; %kg/kmol
mw_h2 = 2.016; %kg.kmol
mw_meoh = 32.04; % kg/kmol
mw_h2o = 18.02; % kg.kmol
mw= [mw_co; mw_co2; mw_h2; mw_meoh; mw_h2o];% kg.kmol
M= n0.*mw; %kg/h
Mt= sum(M, 'all'); %kg/h
M_flowt= Mt/Ntube; %kg/h
%kinetic expression
fun=@(z,y) odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube);
[z,y]=ode15s(fun, zspan, IC);
Warning: Failure at t=1.000016e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.
figure
subplot(311)
yyaxis left
plot(z,y(:,1))
ylabel('y(:,1)')
yyaxis right
plot(z,y(:,2))
ylabel('y(:,2)')
subplot(312)
yyaxis left
plot(z,y(:,3))
ylabel('y(:,3)')
yyaxis right
plot(z,y(:,4))
ylabel('y(:,4)')
subplot(313)
plot(z,y(:,5))
ylabel('y(:,5)')
figure
plot(z,y(:,6));
figure
plot(z,y(:,7));
function dydz=odefun_lovik_axial_pdrop_heattrf(z,y,M_flowt,mw_co,mw_co2,mw_h2,mw_meoh,mw_h2o,A,Dtube)
%%define variables
rho_cat= 1770; %kg/m^3
eps= 0.4;
Dax= 0;
t=0.002;%m
R=8.314;% kJ/kmol.K
dp=0.0042;%m
rp= dp/2;%m
rr=Dtube/2; %m
ro=rr+t; %m
Tc=260+273.15; %K
%%kinetics setting
ka_ref= 0.499;
kb_ref=6.62*10^-11;
kc_ref= 3453.38;
kd_ref= 1.07;
ke_ref= 1.22*10^10;
Ea=17197;
Eb=124119;
Ec=0;
Ed=3669;
Ee=-94765;
%%equilibrium setting;
K1eq_A= 3066;
K2eq_A= 2071;
K1eq_B= 10.592;
K2eq_B= 2.029;
%%kinetic expression
ka=ka_ref.*exp(Ea./(R.*y(7)));
kb=kb_ref.*exp(Eb./(R.*y(7)));
kc= kc_ref.*exp(Ec./(R.*y(7)));
kd= kd_ref.*exp(Ed./(R.*y(7)));
ke= ke_ref.*exp(Ee./(R.*y(7)));
K1eq= 10^((K1eq_A./y(7)) - K1eq_B);
K2eq= 10^((K2eq_A./y(7)) - K2eq_B);
%%to solve
%y(1)=yco
%y(2)=yco2
%y(3)=yh2
%y(4)=ymeoh
%y(5)=yh2o
%y(6)=P
%y(7)=Tr
%1==y(1)+y(2)+y(3)+y(4)+y(5);
%%partial pressures and gas properties
p_co= y(1).*y(6); %bar
p_co2= y(2).*y(6); %bar
p_h2= y(3).*y(6); %bar
p_meoh= y(4).*y(6); %bar
p_h2o= y(5).*y(6); %bar
mu=67.2.*10^-7 + 0.21875.*10^-7.*y(7);
c=y(6).*100./(R.*y(7)); %kmol/m^3
rho_g=c.*(y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2+y(4).*mw_meoh+y(5).*mw_h2o); %kg/m^3
ug=M_flowt./(rho_g.*eps*A.*3600); %m/s
G=M_flowt./(A.*3600); %kg/m^2.s
Re=G.*dp./mu; %Reynold's number
%%thermal properties
ks= 4.18e-3; %...thermal conductivity
kw= 20e-3; %...thermal conductivity
Rec= 1923; %cooling water Re
Prc= 0.843; %cooling water Pr
kg= 0.01234e-3 + (1.84375e-7.*y(7)); %gas thermal conductivity
cp_co= 30.87 + (-1.285e-2.*y(7)) + (2.789e-5.*y(7).^2) + (-1.272e-8.*y(7).^3); %j/mol.K
cp_co2= 19.8 + (7.344e-2.*y(7)) + (-5.602e-5.*y(7).^2) + (1.715e-8.*y(7).^3); %j/mol.K
cp_h2= 27.14 + (9.274e-3.*y(7)) + (-1.381e-5.*y(7).^2) + (7.645e-9.*y(7).^3); %j/mol.K
cp_meoh= 21.15 + (7.092e-2.*y(7)) + (2.587e-5.*y(7).^2) + (-2.852e-8.*y(7).^3); %j/mol.K
cp_h2o= 32.24 + (1.924e-3.*y(7)) + (1.055e-5.*y(7).^2) + (3.596e-9.*y(7).^3); %j/mol.K
cpg= (y(1).*cp_co + y(2).*cp_co2 + y(3).*cp_h2 + y(4).*cp_meoh + y(5).*cp_h2o)./((y(1).*mw_co + y(2).*mw_co2 + y(3).*mw_h2 + y(4).*mw_meoh + y(5).*mw_h2o));
Pr= mu.*cpg./kg; %prandtl number
gamma= ks./kg; %ratio ks:kg
kstat= kg.*((0.5 + 0.493.*((gamma-1)./(22+gamma))).*(gamma./(gamma.*eps + (1-eps))) + 1 - (0.5 + 0.493.*((gamma-1)./(22+gamma))).*(eps+(1-eps).*gamma));
kdyn= (kg.*Re.*Pr)./((7 + 135.8.*(rp/rr)).^2);
krg= kdyn + kstat; %radial thermal conductivity
krs= kstat;
hig= (1.23.*Re.*0.53.*kg)./(2.*rp); %gas inside heat transfer coefficient
his= (0.48 + 0.192.*((rr/rp)-1).^2.*krs)./rp; % inside surface heat transfer coefficient
hoc = (0.020.*Rec.^0.8.*Prc.^(1/3).*(ro/rr).^0.53.*kw)./(2.*rr); %outside cooling water heat transfer coefficient
Uwg= ((1./hig)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of gas
Uws= ((1./his)+((rr/ro)./hoc)+(rr.*log(rr./ro)./(2.*kw))).^-1; %overall heat transfer coefficient of surface
Ueff= ((1./Uwg).*(3.*rr./(8.*krg))).^-1 + ((1./Uws).*(3.*rr./(8.*krs))).^-1; %effective overall heat transfer coefficient of surface
Hr1= 5798 + 35.*(y(7)-498.15);
Hr2= -39892 + 8.*(y(7)-498.15);
%kinetic expression
r1=(kd.*p_co2.*p_h2.*(1-((1./K1eq).*(p_h2o.*p_meoh./(p_h2.^3.*p_co2)))))./(1+kc.*(p_h2o./p_h2)+ka.*p_h2.^0.5+kb.*p_h2o).^3;
r2=(ke.*p_co2.*((1-K2eq.*(p_h2o.*p_co./(p_co2.*p_h2)))))./(1+kc.*(p_h2o./p_h2)+ka.*(p_h2.^0.5)+kb.*p_h2o);
%%differential eq
dydz(1)= ((1-eps).*rho_cat.*(r2 - -2.*r1.*y(1)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(2)= ((1-eps).*rho_cat.*(-1.*r2 + -1.*r1 - -2.*r1.*y(2)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(3)= ((1-eps).*rho_cat.*(-1.*r2 + -3.*r1 - -2.*r1.*y(3)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(4)= ((1-eps).*rho_cat.*(r1 - -2.*r1.*y(4)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(5)= ((1-eps).*rho_cat.*(r2+r1 - -2.*r1.*y(5)))./((c.*ug.*eps).*(1-Dax.*c));
dydz(6)= -10^-5.*(1.75+(150.*(1-eps)./Re)).*((1-eps)./eps.^3).*(ug.^2.*rho_g./dp);
dydz(7)= (1./(G.*cpg)).*((2.*Ueff./rr).*(Tc-y(7))+((1-eps).*rho_cat.*(Hr1+Hr2)));
dydz=dydz';
end

Kategorien

Mehr zu Oil, Gas & Petrochemical finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by