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
Best Answer