Error using daeic12 (line 76) This DAE appears to be of index greater than 1. Error in ode15s (line 310) [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,... Error in PSA_CYCLE_FOR_AIR_SEPERATION (line 255) [t,dudt]=ode15
1 次查看（过去 30 天）
显示 更早的评论
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
%-------------------------------------------------------------
0 个评论
回答（0 个）
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
您也可以从以下列表中选择网站：
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)