Warning: Failure at t=1.130644e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.710505e-20) at time t.

2 visualizaciones (últimos 30 días)
Hello everyone,i am trying to simulated a PSA process with matlab.I am a really new in matlab.I start from trying to simluated the adsorption step which is my first step,i get thiw error and i dont really interstand why.Any help will be much appreciated.I attach the model which trying to simulating
if true
%PROCESS FOR AIR SEPARATION.THE MODEL IS IN DIMENSIONLESS FORM.
% NUMBER OF STEPS:
% ADSORPTION:1
% CO-CURRENT DEPRESSURIZATION:2
% COUNTE CURRENT BLOWDOWN:3
% PURGE:4
% PRESSUREIZATION:5
%for nstep=1:5 %if nstep==1 %ADSORPTION STEP xl=0; xr=1; n=61; dx=(xr-xl)/(n-1); velo_f=0.23; %m/s; P=1.;%bar; Tf=298;%K R=8.314;%J/mol*K P=1.;%BAR
T0=Tf; Tw=Tf;
th0=T0/Tf; thw= Tw/Tf;
%Mixture composition at inlet and at start %inlet x=0 yf1=0.78; yf2=1-yf1;
%at start t=0 y10=0; y20=1-y10; P10=P*y10; P20=P*y20;
%Air:78% N2(1),22% O2(2) cf=1.e5*(P/R*Tf); %mol/m^3;
c1f=cf*yf1/cf; c2f=cf*yf2/cf;
c10=cf*y10/cf; c20=cf*y20/cf;
% Mixture: 1=N2,2=O2 MW1=0.028; %kg/mol MW2=0.032; %kg/mol
%Mixture transport properties in the bulk at P,T conditions DM0=1.5e-5; %m2/s Pref=1.5; %ref.pressure DM=DM0*Pref/P; %m2/s,scales with 1/P mu=1.846e-5; % Pa*s,assumed constant with P
%Fixed bed properties eb=0.4 ; L=2.5; %m Rbed=0.5; % m
%Pellet structural properties ep=0.39; Rp=0.475e-3; % m rho_particle =720/(1-eb); %kg/m3 so that rho_bulk =(1-eb)*rho_particle=720 kg/m3
cp_gas=995.3406; %heat capacity of air J/kgK cp_solid=1171.575; %heat capacity of the solid J/kgK lamda=0.418; %solid thermal conductivity
area_bed=2/Rbed; hw=0;
%---------------------------------------------------------------------------------------- % Equilibrium adsorption isotherm data for air
bp1a=(1.05e-3)/1.e5; % attention bp1a must be in 1/Pa in matlab code
bp1b=2012.92;
bp2a=(3.51e-3)/1.e5;
bp2b=1544.43;
qm1a=1.4796;
qm1b=-0.00167;
qm2a=0.3182;
qm2b=-0.000145;
% Enthalpies of adsorption
DH1=-18033.89029; %J/mol not cal/mol DH2=-13222.06341; %J/mol not cal/mol
% Determine Langmuir parameters at T=T0 bp10=bp1a*exp(bp1b/T0) ; bp20=bp2a*exp(bp2b/T0) ;
qm10=qm1a+qm1b*T0; qm20=qm2a+qm2b*T0;
Ke10=qm10*bp10 ; Ke20=qm20*bp20 ;
% Determine Langmuir parameters at T=Tf bp1f=bp1a*exp(bp1b/Tf) ; bp2f=bp2a*exp(bp2b/Tf) ;
qm1f=qm1a+qm1b*Tf; qm2f=qm2a+qm2b*Tf; Ke1f=qm1f*bp1f ; Ke2f=qm2f*bp2f ;
AMW=yf1*MW1+yf2*MW2; %kg/mol rho_gas=cf*AMW; tr=L/velo_f;
qe1f = Ke1f*(P*10^5)*yf1 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 ); qe2f = Ke2f*(P*10^5)*yf2 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 );
qe10 = Ke10*(P*10^5)*y10 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 ); qe20 = Ke20*(P*10^5)*y20 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 ); %---------------------------------------------------------------------------- De1=2.22e-6; De2=De1*sqrt(MW1/MW2);
Re=(rho_gas*(eb*velo_f)*2*Rp)/mu; Sc=mu/((rho_gas*DM)); Sh=2+1.1*(Sc^(1/3))*(Re^(0.6)); kf=(Sh*DM)/(2*Rp); tdif1=1/(15*ep*De1/Rp^2); tdif2=1/(15*ep*De2/Rp^2); tfilm=Rp/(3*kf); k1inv=tdif1+tfilm; k1=1/k1inv; k2inv=tdif2+tfilm; k2=1/k2inv;
%k1=1000; %k2=1000;
%"LINE DRIVING FORCE RATE APPROXIMATION GAS SEPARATION BY ADSORPTION %PROCESSES",RALPH T.YANG,pg.296,eq.8.47 kc1=15*De1/(Rp^2); kc2=15*De2/(Rp^2); %---------------------------------------------------------------------------- %kc1 = 1000; % N2 mass transfer in the micropores (zeolite crystals) %kc2 = 1000; % O2 mass transfer in the micropores (zeolite crystals)
% Check these correlations
Dz = velo_f/667; lamz= Dz*cp_gas; lame=lamz+lamda/eb;
%tdisp = (Dz/uf^2)*(1-eb)/eb
% Here we non-dimensionalize wrt tr=L/u0 Pe=velo_f*L/Dz; Peth=rho_gas*cp_gas*velo_f*L/lame;
% I put these values to comply with Rege and Fortran file Pe=150; Peth=2; %-------------------------------------
Bi=hw*area_bed*tr/eb*rho_gas*cp_gas; Le=((1-eb)/eb)*rho_particle*cp_solid/(rho_gas*cp_gas); F=((1-eb)/eb)*(tr/c1f)*(3/Rp)*(R/P); A1=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe1f; A2=Bi/(1+Le); A3=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe20; dPdt=0;
param =[Pe Peth Dz DM A1 A2 F rho_particle rho_gas k1 k2 kc1 kc2 ep tr... DH1 DH2 eb Le n thw L cp_gas Rp c1f Bi dPdt bp1a bp2a bp1b bp2b... Tf qm1a qm2a qm1b qm2b qe1f qe2f qe10 qe20 A3 velo_f P cf R]';
for i=1:n c10(i)=0.78; velo_0(i)=1; cp10(i)=c10(i); cp20(i)=1-c10(i); q10(i) =qe10/qe1f; q20(i) =qe20/qe20; th0(i) =T0/Tf; end
u0=[c10 velo_0 cp10 cp20 q10 q20 th0].';
for i=1:n if i==1 x(i) = 0; else x(i) = x(i-1)+dx ; end end t0=0; %lower bound for time tf=10; %upper bound for time m=101; %time points tspan=linspace(t0,tf,m); [t,dudt]=ode15s(@(t,u)ADSORPTION_fun(u,dx,param ),tspan,u0); function [ dudt ] = ADSORPTION_fun(u,dx,param ) %UNTITLED7 Summary of this function goes here %----------------------------------------------------------- %PARAMETERS Pe= param(1,1); Peth=param(2,1); Dz =param(3,1); DM =param(4,1); A1 =param(5,1); A2 =param(6,1); F =param(7,1); rho_particle =param(8,1); rho_gas=param(9,1); k1=param(10,1); k2=param(11,1); kc1=param(12,1); kc2=param(13,1); ep=param(14,1); tr=param(15,1); DH1=param(16,1); DH2=param(17,1); eb=param(18,1); Le=param(19,1); n=param(20,1); thw=param(21,1); L=param(22,1); cp_gas=param(23,1); Rp=param(24,1); c1f=param(25,1); Bi=param(26,1); dPdt=param(27,1); bp1a=param(28,1); bp2a=param(29,1); bp1b=param(30,1); bp2b=param(31,1); Tf=param(32,1); qm1a=param(33,1); qm2a=param(34,1); qm1b=param(35,1); qm2b=param(36,1); qe1f=param(37,1); qe2f=param(38,1); qe10=param(39,1); qe20=param(40,1); A3=param(41,1); velo_f=param(42,1); P=param(43,1); cf=param(44,1); R=param(45,1);
const1=((1-eb)/eb)*tr; const2=1/(Peth*(1+Le)); const3=(1-(1/(1+Le))); dPdt=0; const4=(tr/c1f)*rho_particle*qe1f; const5=(tr/(1-c1f))*rho_particle*qe20; const6=((1-eb)/(eb*(1+Le)))*(1/(rho_gas*cp_gas*Tf)); %------------------------------------------------------------- dudt=zeros(7*n,1); dc1dt=zeros(7*n,1); dvelodt=zeros(7*n,1); dcp1dt=zeros(7*n,1); dcp2dt=zeros(7*n,1); dq1dt=zeros(7*n,1); dq2dt=zeros(7*n,1); dTdt=zeros(7*n,1); %--------------------------------------------------------------- c1(1:n)=u(1:n); velo(1:n)=u(n+1:2*n); cp1(1:n)=u(2*n+1:3*n); cp2(1:n)=u(3*n+1:4*n); q1(1:n)=u(4*n+1:5*n); q2(1:n)=u(5*n+1:6*n); T(1:n)=u(6*n+1:7*n); %------------------------------------------------------- %qe1 = (qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) ); %qe2 = (qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) ); %------------------------------------------------------- for i=1:n if i==1 %BOUNDARY CONDITIONS AT X=0 %--------------------------------------------------------- dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1f)/(dx^2))... -velo(i)*((c1(i+1)-c1f)/2*dx)... -c1(i)*((velo(i+1)-velo_f)/2*dx)... +const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+Tf)/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+Tf)/(dx^2) ...
+velo(i)*const3*(T(i+1)-Tf)/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo_f)/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-Tf)/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
(1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i)+...
(bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))
+A2*(thw-T(i));
%---------------------------------------------------------
elseif i==2:n-1
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1(i-1))/(dx^2))...
-velo(i)*((c1(i+1)-c1(i-1))/2*dx)...
-c1(i)*((velo(i+1)-velo(i-1))/2*dx)...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+T(i-1))/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+T(i-1))/(dx^2) ...
+velo(i)*const3*(T(i+1)-T(i-1))/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo(i-1))/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(Tw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i))/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-T(i-1))/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i)...
/( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe1f)-q1(i))
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe20)-q2(i))
+A2*(thw-T(i));
elseif i==n
%BOUNDARY CONDITIONS AT X=L
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((-2*c1(i)+2*c1(i-1))/(dx^2))...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-(1/Pe)*(-2*T(i)+T(i-1))/(dx^2) ...
-(T(i)/P)*dPdt...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)= ((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*...
(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))...
+A2*(thw-T(i));
%--------------------------------------------------------
end
dudt(i)=dc1dt(i);
dudt(n+i)=dvelodt(i);
dudt(2*n+i)=dcp1dt(i);
dudt(3*n+i)=dcp2dt(i);
dudt(4*n+i)=dq1dt(i);
dudt(5*n+i)=dq2dt(i);
dudt(6*n+i)=dTdt(i);
end %-------------------------------------------------------------
end
end
  1 comentario
Jan
Jan el 14 de Ag. de 2018
Please use the "{} Code" button to format the code. Post the complete error message, such that the readers can know in which line of which function the problem occurs.
Note: 10^5 is a very expensive power operation, while 1e5 is a cheap constant.

Iniciar sesión para comentar.

Respuesta aceptada

ANDREAS CHARALAMBOUS
ANDREAS CHARALAMBOUS el 14 de Ag. de 2018
if true
xl=0;
xr=1;
n=61;
dx=(xr-xl)/(n-1);
velo_f=0.23; %m/s;
Tf=298;%K
R=8.314;%J/mol*K
P=1.;%BAR
T0=Tf;
Tw=Tf;
th0=T0/Tf;
thw= Tw/Tf;
%Mixture composition at inlet and at start
%inlet x=0
yf1=0.78;
yf2=1-yf1;
%at start t=0
y10=0;
y20=1-y10;
P10=P*y10;
P20=P*y20;
%Air:78% N2(1),22% O2(2)
cf=1.e5*(P/R*Tf); %mol/m^3;
c1f=cf*yf1/cf;
c2f=cf*yf2/cf;
c10=cf*y10/cf;
c20=cf*y20/cf;
% Mixture: 1=N2,2=O2
MW1=0.028; %kg/mol
MW2=0.032; %kg/mol
%Mixture transport properties in the bulk at P,T conditions
DM0=1.5e-5; %m2/s
Pref=1.5; %ref.pressure
DM=DM0*Pref/P; %m2/s,scales with 1/P
mu=1.846e-5; % Pa*s,assumed constant with P
%Fixed bed properties
eb=0.4 ;
L=2.5; %m
Rbed=0.5; % m
%Pellet structural properties
ep=0.39;
Rp=0.475e-3; % m
rho_particle =720/(1-eb); %kg/m3 so that rho_bulk =(1-eb)*rho_particle=720 kg/m3
cp_gas=995.3406; %heat capacity of air J/kgK
cp_solid=1171.575; %heat capacity of the solid J/kgK
lamda=0.418; %solid thermal conductivity
area_bed=2/Rbed;
hw=0;
%----------------------------------------------------------------------------------------
% Equilibrium adsorption isotherm data for air
bp1a=(1.05e-3)/1.e5; % attention bp1a must be in 1/Pa in matlab code
bp1b=2012.92;
bp2a=(3.51e-3)/1.e5;
bp2b=1544.43;
qm1a=1.4796;
qm1b=-0.00167;
qm2a=0.3182;
qm2b=-0.000145;
% Enthalpies of adsorption
DH1=-18033.89029; %J/mol not cal/mol
DH2=-13222.06341; %J/mol not cal/mol
% Determine Langmuir parameters at T=T0
bp10=bp1a*exp(bp1b/T0) ;
bp20=bp2a*exp(bp2b/T0) ;
qm10=qm1a+qm1b*T0;
qm20=qm2a+qm2b*T0;
Ke10=qm10*bp10 ;
Ke20=qm20*bp20 ;
% Determine Langmuir parameters at T=Tf
bp1f=bp1a*exp(bp1b/Tf) ;
bp2f=bp2a*exp(bp2b/Tf) ;
qm1f=qm1a+qm1b*Tf;
qm2f=qm2a+qm2b*Tf;
Ke1f=qm1f*bp1f ;
Ke2f=qm2f*bp2f ;
AMW=yf1*MW1+yf2*MW2; %kg/mol
rho_gas=cf*AMW;
tr=L/velo_f;
qe1f = Ke1f*(P*10^5)*yf1 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 );
qe2f = Ke2f*(P*10^5)*yf2 /( 1 + bp1f*(P*10^5)*yf1 + bp2f*(P*10^5)*yf2 );
qe10 = Ke10*(P*10^5)*y10 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 );
qe20 = Ke20*(P*10^5)*y20 /( 1 + bp10*(P*10^5)*y10 + bp20*(P*10^5)*y20 );
%----------------------------------------------------------------------------
De1=2.22e-6;
De2=De1*sqrt(MW1/MW2);
Re=(rho_gas*(eb*velo_f)*2*Rp)/mu;
Sc=mu/((rho_gas*DM));
Sh=2+1.1*(Sc^(1/3))*(Re^(0.6));
kf=(Sh*DM)/(2*Rp);
tdif1=1/(15*ep*De1/Rp^2);
tdif2=1/(15*ep*De2/Rp^2);
tfilm=Rp/(3*kf);
k1inv=tdif1+tfilm;
k1=1/k1inv;
k2inv=tdif2+tfilm;
k2=1/k2inv;
%k1=1000;
%k2=1000;
%"LINE DRIVING FORCE RATE APPROXIMATION GAS SEPARATION BY ADSORPTION
%PROCESSES",RALPH T.YANG,pg.296,eq.8.47
kc1=15*De1/(Rp^2);
kc2=15*De2/(Rp^2);
%----------------------------------------------------------------------------
%kc1 = 1000; % N2 mass transfer in the micropores (zeolite crystals)
%kc2 = 1000; % O2 mass transfer in the micropores (zeolite crystals)
% Check these correlations
Dz = velo_f/667;
lamz= Dz*cp_gas;
lame=lamz+lamda/eb;
%tdisp = (Dz/uf^2)*(1-eb)/eb
% Here we non-dimensionalize wrt tr=L/u0
Pe=velo_f*L/Dz;
Peth=rho_gas*cp_gas*velo_f*L/lame;
Pe=150;
Peth=2;
%-------------------------------------
Bi=hw*area_bed*tr/eb*rho_gas*cp_gas;
Le=((1-eb)/eb)*rho_particle*cp_solid/(rho_gas*cp_gas);
F=((1-eb)/eb)*(tr/c1f)*(3/Rp)*(R/P);
A1=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe1f;
A2=Bi/(1+Le);
A3=(1/(1+Le))+((1-eb)/eb)*(1/(rho_gas*cp_gas*Tf))*qe20;
dPdt=0;
param =[Pe Peth Dz DM A1 A2 F rho_particle rho_gas k1 k2 kc1 kc2 ep tr...
DH1 DH2 eb Le n thw L cp_gas Rp c1f Bi dPdt bp1a bp2a bp1b bp2b...
Tf qm1a qm2a qm1b qm2b qe1f qe2f qe10 qe20 A3 velo_f P cf R]';
for i=1:n
c10(i)=0;
velo_0(i)=1;
cp10(i)=c10(i);
cp20(i)=1-c10(i);
q10(i) =qe10/qe1f;
q20(i) =qe20/qe20;
th0(i) =T0/Tf;
end
u0=[c10 velo_0 cp10 cp20 q10 q20 th0].';
for i=1:n
if i==1
x(i) = 0;
else
x(i) = x(i-1)+dx ;
end
end
t0=0; %lower bound for time
tf=10; %upper bound for time
m=101; %time points
tspan=linspace(t0,tf,m);
[t,dudt]=ode15s(@(t,u)ADSORPTION_fun(u,dx,param ),tspan,u0);
Pe= param(1,1);
Peth=param(2,1);
Dz =param(3,1);
DM =param(4,1);
A1 =param(5,1);
A2 =param(6,1);
F =param(7,1);
rho_particle =param(8,1);
rho_gas=param(9,1);
k1=param(10,1);
k2=param(11,1);
kc1=param(12,1);
kc2=param(13,1);
ep=param(14,1);
tr=param(15,1);
DH1=param(16,1);
DH2=param(17,1);
eb=param(18,1);
Le=param(19,1);
n=param(20,1);
thw=param(21,1);
L=param(22,1);
cp_gas=param(23,1);
Rp=param(24,1);
c1f=param(25,1);
Bi=param(26,1);
dPdt=param(27,1);
bp1a=param(28,1);
bp2a=param(29,1);
bp1b=param(30,1);
bp2b=param(31,1);
Tf=param(32,1);
qm1a=param(33,1);
qm2a=param(34,1);
qm1b=param(35,1);
qm2b=param(36,1);
qe1f=param(37,1);
qe2f=param(38,1);
qe10=param(39,1);
qe20=param(40,1);
A3=param(41,1);
velo_f=param(42,1);
P=param(43,1);
cf=param(44,1);
R=param(45,1);
const1=((1-eb)/eb)*tr;
const2=1/(Peth*(1+Le));
const3=(1-(1/(1+Le)));
dPdt=0;
const4=(tr/c1f)*rho_particle*qe1f;
const5=(tr/(1-c1f))*rho_particle*qe20;
const6=((1-eb)/(eb*(1+Le)))*(1/(rho_gas*cp_gas*Tf));
%-------------------------------------------------------------
dudt=zeros(7*n,1);
dc1dt=zeros(7*n,1);
dvelodt=zeros(7*n,1);
dcp1dt=zeros(7*n,1);
dcp2dt=zeros(7*n,1);
dq1dt=zeros(7*n,1);
dq2dt=zeros(7*n,1);
dTdt=zeros(7*n,1);
%---------------------------------------------------------------
c1(1:n)=u(1:n);
velo(1:n)=u(n+1:2*n);
cp1(1:n)=u(2*n+1:3*n);
cp2(1:n)=u(3*n+1:4*n);
q1(1:n)=u(4*n+1:5*n);
q2(1:n)=u(5*n+1:6*n);
T(1:n)=u(6*n+1:7*n);
%-------------------------------------------------------
%qe1 = (qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) );
%qe2 = (qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) );
%-------------------------------------------------------
for i=1:n
if i==1
%BOUNDARY CONDITIONS AT X=0
%---------------------------------------------------------
dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1f)/(dx^2))...
-velo(i)*((c1(i+1)-c1f)/2*dx)...
-c1(i)*((velo(i+1)-velo_f)/2*dx)...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+Tf)/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+Tf)/(dx^2) ...
+velo(i)*const3*(T(i+1)-Tf)/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo_f)/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+Tf)/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-Tf)/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
(1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp1(i)+...
(bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))
+A2*(thw-T(i));
%---------------------------------------------------------
elseif i==2:n-1
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1(i-1))/(dx^2))...
-velo(i)*((c1(i+1)-c1(i-1))/2*dx)...
-c1(i)*((velo(i+1)-velo(i-1))/2*dx)...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+T(i-1))/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+T(i-1))/(dx^2) ...
+velo(i)*const3*(T(i+1)-T(i-1))/(2*dx)...
-(T(i)/P)*dPdt...
-T(i)*(velo(i+1)-velo(i-1))/(2*dx)...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(Tw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)=((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i))/qe20)-q2(i));
dq1dt(i)=kc1*tr(i)*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(T(i+1)-2*T(i)+T(i-1))/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-T(i-1))/(2*dx)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i)...
/( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe1f)-q1(i))
-const6*kc2*tr*qe20*DH2*(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i))/qe20)-q2(i))
+A2*(thw-T(i));
elseif i==n
%BOUNDARY CONDITIONS AT X=L
%--------------------------------------------------------
dc1dt(i)=(1/Pe)*((-2*c1(i)+2*c1(i-1))/(dx^2))...
+const1*k1*(c1(i)-cp1(i));
dvelodt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-(1/Pe)*(-2*T(i)+T(i-1))/(dx^2) ...
-(T(i)/P)*dPdt...
-A1*dq1dt(i)...
-A3*dq2dt(i)...
+A2*(thw-T(i))...
+F*(T(i)^2)*k1*tr*(c1(i)-cp1(i))...
+F*(T(i)^2)*k2*tr*((1-c1(i))-cp2(i));
dcp1dt(i)=((k1*tr)/ep)*(c1(i)-cp1(i))
-kc1*const4*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dcp2dt(i)= ((k2*tr)/ep)*((1-c1(i))-cp1(i))
-kc2*const5*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dq1dt(i)=kc1*tr*(((qm1a+qm1b*(Tf*T(i)))*bp1a*exp(bp1b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe1f)-q1(i));
dq2dt(i)=kc2*tr*(((qm2a+qm2b*(Tf*T(i)))*bp2a*exp(bp2b/(Tf*T(i)))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) /( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*...
(cf*R*Tf*T(i)*10^5)*cp2(i) )/qe20)-q2(i));
dTdt(i)=const2*(-2*T(i)+2*T(i-1))/(dx^2)...
-const6*kc1*tr*qe1f*DH1*(((qm1a+qm1b*(Tf*T(i)))*...
bp1a*exp(bp1b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp1(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe1f)-q1(i))...
-const6*kc2*tr*qe20*DH2*...
(((qm2a+qm2b*(Tf*T(i)))*...
bp2a*exp(bp2b/(Tf*T(i)))*(cf*R*Tf*T(i)*10^5)*cp2(i) /...
( 1 + (bp1a*exp(bp1b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp1(i) + (bp2a*exp(bp2b/(Tf*T(i))))*(cf*R*Tf*T(i)*10^5)*...
cp2(i) )/qe20)-q2(i))...
+A2*(thw-T(i));
%--------------------------------------------------------
end
dudt(i)=dc1dt(i);
dudt(n+i)=dvelodt(i);
dudt(2*n+i)=dcp1dt(i);
dudt(3*n+i)=dcp2dt(i);
dudt(4*n+i)=dq1dt(i);
dudt(5*n+i)=dq2dt(i);
dudt(6*n+i)=dTdt(i);
end %-------------------------------------------------------------
end
  2 comentarios
ANDREAS CHARALAMBOUS
ANDREAS CHARALAMBOUS el 14 de Ag. de 2018
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. > In ode15s (line 589) In PSA_CYCLE_FOR_AIR_SEPERATION (line 208) Warning: Failure at t=9.295078e-06. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.710505e-20) at time t. > In ode15s (line 668) In PSA_CYCLE_FOR_AIR_SEPERATION (line 208) >>

Iniciar sesión para comentar.

Más respuestas (1)

Jan
Jan el 14 de Ag. de 2018
Although the code is hard to read yet, the error "Unable to meet integration tolerances" means most likely that the function to be integrated has a pole or a discontinuity at this point. Simply try this by integrating until 1.1306e-05 and plot the results. What do you observe?

Categorías

Más información sobre General Physics en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by