Error using daeic12 (line 76) This DAE appears to be of index greater than 1. Error in ode15s (line 310) [y,yp,f0,d​fdy,nFE,nP​D,Jfac] = daeic12(od​eFcn,odeAr​gs,t,ICtyp​e,Mt,y,yp0​,f0,... Error in PSA_CYCLE_​FOR_AIR_SE​PERATION (line 255) [t,dudt]=ode15

1 view (last 30 days)
ANDREAS CHARALAMBOUS
ANDREAS CHARALAMBOUS on 19 Aug 2018
Edited: Walter Roberson on 20 Aug 2018
Hello everyone,i am trying to solve a dae system of equations which resulting of a space dicr. with central differences. I have 6 odes and 1 algebraic equation but the space dicr. leads to 427 equations. I create a mass matrix but i get the error that my dae is index greater than 1. I have to reduce index but i don't how to do it with so much number of equations and variables. I attach my equations above and my code it as follow
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=(P/R*Tf); %mol/m^3;
cff=1e5*cf;
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 cff].';
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
%------------------------------------------------------------
%Create the mass matrix in order to solve the coupled differential and
%algebraic equation resulting from the space dicretization
for i=1:7*n
for j=1:7*n
if i<=n && j<=n && i==j
M(i,j)=1;
elseif i>n && j>n && i<=2*n && j<=2*n &&i==j
M(i,j)=0;
elseif i>2*n && j>2*n && i<=3*n && j<=3*n && i==j
M(i,j)=ep;
elseif i>3*n && j>3*n && i<=4*n && j<=4*n && i==j
M(i,j)=ep;
elseif i>4*n && j>4*n && i<=5*n && j<=5*n && i==j
M(i,j)=1;
elseif i>5*n && j>5*n && i<=6*n && j<=6*n && i==j
M(i,j)=1;
elseif i>6*n && j>6*n && i<=7*n && j<=7*n && i==j
M(i,j)=1;
end
end
end
%------------------------------------------------------------
u=zeros(1,7*n); %matrix of solutions
for i=1:n
c1(i)=0;
velo(i)=0;
cp1(i)=0;
cp2(i)=0;
q1(i)=0;
q2(i)=0;
T(i)=0;
end
u(1:n)=c1(1:n);
u(n+1:2*n)=velo(1:n);
u(2*n+1:3*n)=cp1(1:n);
u(3*n+1:4*n)=cp2(1:n);
u(4*n+1:5*n)=q1(1:n);
u(5*n+1:6*n)=q2(1:n);
u(6*n+1:7*n)=T(1:n);
tspan=linspace(t0,tf,m);
options=odeset('Mass',M);
[t,dudt]=ode15s(@(t,u)ADSORPTION_fun(u,n,dx,param ),tspan,u0,options);
function [ dudt ] = ADSORPTION_fun(u,n,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);
cff=param(46,1);
const1=((1-eb)/eb)*tr;
const2=1/(Peth*(1+Le));
const3=(1-(1/(1+Le)));
dPdt=0;
const4=((tr/cff)*rho_particle*qe1f);
const5=((tr/cff)*rho_particle*qe20);
const6=((1-eb)/(eb*(1+Le)))*(1/(rho_gas*cp_gas*Tf));
%-------------------------------------------------------------
dudt=zeros(7*n,1);
dc1dt=zeros(n,1);
dvelodt=zeros(n,1);
dcp1dt=zeros(n,1);
dcp2dt=zeros(n,1);
dq1dt=zeros(n,1);
dq2dt=zeros(n,1);
dTdt=zeros(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
dc1dt(i)=(1/Pe)*((c1(i+1)-2*c1(i)+c1f/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));
elseif i==n
dc1dt(i)=(1/Pe)*((-2*c1(i)+2*c1(i-1))/(dx^2))...
+const1*k1*(c1(i)-cp1(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));
end
end
for i=1:n
if i==1
dvelodt(i)=const2*(T(i+1)-2*T(i)+Tf/Tf)/(dx^2)...
+(2/(Pe*T(i)))*(T(i+1)+Tf/Tf)/(2*dx)...
-(1/Pe)*(T(i+1)-2*T(i)+Tf/Tf)/(dx^2) ...
+velo(i)*const3*(T(i+1)-Tf/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-(1/T(i)))-cp2(i));
elseif i==n
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-(1/T(i)))-cp2(i));
elseif i==2:n-1
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-(1/T(i)))-cp2(i));
end
end
for i=1:n
if i==1
dcp1dt(i)=(k1*tr)*(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));
elseif i==n
dcp1dt(i)=(k1*tr)*(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));
elseif i==2:n-1
dcp1dt(i)=(k1*tr)*(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));
end
end
for i=1:n
if i==1
dcp2dt(i)=(k2*tr)*((1-(1/T(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));
elseif i==n
dcp2dt(i)= (k2*tr)*((1-(1/T(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));
elseif i==2:n-1
dcp2dt(i)=(k2*tr)*((1-(1/T(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));
end
end
for i=1:n
if i==1
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));
elseif i==n
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));
elseif i==2:n-1
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));
end
end
for i=1:n
if i==1
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));
elseif i==n
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));
elseif i==2:n-1
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));
end
end
for i=1:n
if i==1
dTdt(i)=const2*(T(i+1)-2*T(i)+Tf/Tf)/(dx^2)...
-(velo(i)/(1+Le))*(T(i+1)-Tf/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==n
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));
elseif i==2:n-1
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));
end
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
%-------------------------------------------------------------

Sign in to answer this question.

Answers (0)

Categories

Find more on General Physics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!