使用s函数时提示call must be a real vector
提示：State derivatives returned by S-function 'sfunmc' in 'mcFC/S-Function' during flag=1 call must be a real vector of length 160.检查后发现是在else中的log内出现了负数，可是为什么会出现负数呢，如何更改呢。程序比较长不知道能不能显示完全，非常感谢您的关注，我的qq是5471282，如果哪位大虾能帮助解决不胜感激，可以请你吃饭或其他有偿的形式表达我的感谢，再次感谢。sizes = simsizes;
sizes.NumContStates= 16*N;
sizes.NumDiscStates= 0;
sizes.NumOutputs = 16*N;
sizes.NumInputs = 16;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1; % at least one sample time is needed
sys = simsizes(sizes);
x0= [];
function sys=mdlDerivatives(t,x,u,N)
L=0.4;%为电池长度
W=0.1;%为电池宽度
h=L/N;%N为离散的点数
for i=1:N
F=96487;%法拉第常数
R=8.314;
hf=0.0008;%燃料流道的高度
Af=W*hf;%通道截面积
f=hf/W;
f_w=24*(1-1.3553*f+1.9467*f^2-1.0712*f^3+0.9564*f^4-0.2537*f^5);%计算f(W,h)
dhf=(2*W*hf)/(W+hf); % 阳极气流通道的hydraulic diameter
if i==1
roua(i)=u(14)*0.032+u(15)*0.028+u(16)*0.044;
rouf(i)=u(9)*0.016+u(10)*0.002+u(11)*0.028+u(12)*0.044+u(13)*0.018;% the unit of rou
nd_ch4=-9.9989+529.37*(u(5)/1000)-543.82*(u(5)/1000).^2+548.11*(u(5)/1000).^3-367.06*(u(5)/1000).^4+140.48*(u(5)/1000).^5-22.92*(u(5)/1000).^6;
nd_co2=-20.434+680.07*(u(5)/1000)-432.49*(u(5)/1000).^2+244.22*(u(5)/1000).^3-85.929*(u(5)/1000).^4+14.45*(u(5)/1000).^5-0.4564*(u(5)/1000).^6;
nd_h2o=-6.7541+244.93*(u(5)/1000)+419.5*(u(5)/1000).^2-522.38*(u(5)/1000).^3+348.12*(u(5)/1000).^4-126.96*(u(5)/1000).^5+19.591*(u(5)/1000).^6;
nd_co=-4.9137+793.65*(u(5)/1000)+875.9*(u(5)/1000).^2+883.75*(u(5)/1000).^3-572.14*(u(5)/1000).^4+208.42*(u(5)/1000).^5-32.298*(u(5)/1000).^6;
nd_h2=15.553+299.78*(u(5)/1000)-244.34*(u(5)/1000).^2+249.41*(u(5)/1000).^3-167.51*(u(5)/1000).^4+62.966*(u(5)/1000).^5-9.9892*(u(5)/1000).^6;
mol_ch4=u(9)/(u(9)+u(10)+u(11)+u(12)+u(13));
mol_h2=u(10)/(u(9)+u(10)+u(11)+u(12)+u(13));
mol_co=u(11)/(u(9)+u(10)+u(11)+u(12)+u(13));
mol_co2=u(12)/(u(9)+u(10)+u(11)+u(12)+u(13));
mol_h2o=u(13)/(u(9)+u(10)+u(11)+u(12)+u(13));
nd_fuel= (nd_co2*mol_co2*0.044^(0.5)+nd_h2o*mol_h2o*0.018^(0.5)+nd_ch4*mol_ch4*0.016^(0.5)+nd_co*mol_co*0.028^(0.5)+nd_h2*mol_h2*0.002^(0.5))*10^(-7)/(mol_co2*0.044^(0.5)+mol_h2o*0.018^(0.5)+mol_ch4*0.016^(0.5)+mol_co*0.028^(0.5)+mol_h2*0.002^(0.5));
cp_co2=4.3669+204.6*(u(5)/1000)-471.33*(u(5)/1000).^2+657.88*(u(5)/1000).^3-519.9*(u(5)/1000).^4+214.58*(u(5)/1000).^5-35.992*(u(5)/1000).^6;
cp_h2o=37.737-41.205*(u(5)/1000)+146.01*(u(5)/1000).^2-217.08*(u(5)/1000).^3+181.54*(u(5)/1000).^4-79.409*(u(5)/1000).^5+14.105*(u(5)/1000).^6;
cp_co=30.429-8.1781*(u(5)/1000)+5.2062*(u(5)/1000).^2+41.974*(u(5)/1000).^3-66.346*(u(5)/1000).^4+37.456*(u(5)/1000).^5-7.6538*(u(5)/1000).^6;
cp_ch4=47.964-178.59*(u(5)/1000)+712.55*(u(5)/1000).^2-1068.7*(u(5)/1000).^3+856.93*(u(5)/1000).^4-358.75*(u(5)/1000).^5+61.321*(u(5)/1000).^6;
cp_h2=21.157+56.036*(u(5)/1000)-150.55*(u(5)/1000).^2+199.29*(u(5)/1000).^3-136.15*(u(5)/1000).^4+46.903*(u(5)/1000).^5-6.4725*(u(5)/1000).^6;
cp_mol=mol_ch4*cp_ch4+mol_h2*cp_h2+mol_co*cp_co+mol_co2*cp_co2+mol_h2o*cp_h2o;
cpf=1000*cp_mol/(mol_ch4*16+ mol_h2*2+mol_co*28+mol_co2*44+mol_h2o*18);
Mf=(16*mol_ch4+2*mol_h2+28*mol_co+44*mol_co2+18*mol_h2o)/1000
KfPEN=4.25;% J/s*m2*K
C11=KfPEN/(rouf(i)*cpf*hf);
C12=C11;
nd_co2_ca=-20.434+680.07*(u(6)/1000)-432.49*(u(6)/1000).^2+244.22*(u(6)/1000).^3-85.929*(u(6)/1000).^4+14.45*(u(6)/1000).^5-0.4564*(u(6)/1000).^6;
nd_o2=-1.6918+889.75*(u(6)/1000)-892.79*(u(6)/1000).^2+905.98*(u(6)/1000).^3-598.36*(u(6)/1000).^4+221.64*(u(6)/1000).^5-34.754*(u(6)/1000).^6;
nd_n2=1.2719+771.45*(u(6)/1000)-809.2*(u(6)/1000).^2+832.47*(u(6)/1000).^3-553.93*(u(6)/1000).^4+206.15*(u(6)/1000).^5-32.43*(u(6)/1000).^6;
mol_co2_ca=u(16)/(u(14)+u(15)+u(16));
mol_o2=u(14)/(u(14)+u(15)+u(16));
mol_n2=u(15)/(u(14)+u(15)+u(16));
nd_oxid= (nd_co2_ca*mol_co2_ca*0.044^(0.5)+nd_o2*mol_o2*0.032^(0.5)+nd_n2*mol_n2*0.028^(0.5))*10^(-7)/(mol_co2_ca*0.044^(0.5)+mol_o2*0.032^(0.5)+mol_n2*0.028^(0.5));
cp_co2=4.3669+204.6*(u(6)/1000)-471.33*(u(6)/1000).^2+657.88*(u(6)/1000).^3-519.9*(u(6)/1000).^4+214.58*(u(6)/1000).^5-35.992*(u(6)/1000).^6;
cp_o2=34.85-57.975*(u(6)/1000)+203.68*(u(6)/1000).^2-300.37*(u(6)/1000).^3+231.72*(u(6)/1000).^4-91.821*(u(6)/1000).^5+14.776*(u(6)/1000).^6;
cp_n2=29.027+4.8987*(u(6)/1000)-38.04*(u(6)/1000).^2+105.17*(u(6)/1000).^3-113.56*(u(6)/1000).^4+55.554*(u(6)/1000).^5-10.35*(u(6)/1000).^6;
cp_mol_ox=mol_o2*cp_o2+mol_n2*cp_n2+mol_co2_ca*cp_co2;
cpa=1000*cp_mol_ox/(mol_co2_ca*44+ mol_o2*32+mol_n2*28);
Ma=(44*mol_co2_ca+32*mol_o2+28*mol_n2)/1000
ha=0.0008
KaPEN=10.1;%w/m2*k
C21=KaPEN/(roua(i)*cpa*ha);
C22=C21;
lmdPEN=22;
rouPEN=2000;
cpPEN=800;
C31=lmdPEN/(rouPEN*cpPEN);
hPEN=0.001;%pen 厚度
C32=KfPEN/(rouPEN*cpPEN*hPEN);
C33=KaPEN/(rouPEN*cpPEN*hPEN);
dh1= (206*10^3 - 16379.71 + 8.314*(7.951*u(5) - 4.354*10^(-3)*u(5)^2 + 0.7213*10^(-6)*u(5)^3 - 0.097*10^5/u(5)));
dh2= (-41.1*10^3 - 7657.295 + 8.314*(1.86*u(5) - 0.27*10^(-3)*u(5)^2 + 1.164*10^5/u(5)));
dh3= (-241*10^3 - 4098.471 + 8.314*(-1.5985*u(7) + 0.3875*10^(-3)*u(7)^2 - 0.1515*10^5/u(7)))
Ca=2.04*10^(-3);
C1=3.28*10^(-9);
C2=3.39*10^(-6);
Cir=1.12*10^(-2);
detaHc1=132000;
detaHc2=67100;
detaHa=23700;
detaHir=23000;
vopt=0.7;
detaG=-242000+45.8*u(7);
E=-detaG/(2*F);
Ph2_an=u(3)*mol_h2*10^(-5);%
Po2_ca=u(4)*mol_o2*10^(-5);
Pco2_ca=u(4)*mol_co2_ca*10^(-5);
Ph2o_an=u(3)*mol_h2o*10^(-5);
Pco2_an=u(3)*mol_co2*10^(-5);
Nerst=R*u(7)*log((Ph2_an*Pco2_ca*Po2_ca^0.5)/(Ph2o_an*Pco2_an))/(2*F);
Ran=Ca*(Ph2_an)^(-0.5)*exp(detaHa/(R*u(7)));
Rca=C1*(Po2_ca)^(-0.75)*(Pco2_ca)^(0.5)*exp(detaHc1/(R*u(7)))+C2*exp(detaHc2/(R*u(7)))*mol_co2_ca^(-1);
Rir=Cir*exp(detaHir/(R*u(7)));% change the pressure unit from pa to atm
Rtol=Ran+Rca+Rir;
j=10000*(E-vopt-Nerst)/Rtol;%current density unit is A/m2
k0=4274*10^(-5);
Ea=82000;
kwgsr=2.35;%
Keq=exp(4276/u(5)-3.961);
P_ch4=u(3)*mol_ch4;%
P_co=u(3)*mol_co;%
rate1 = k0*P_ch4*exp(- Ea/(R*u(5)));
rate2 = kwgsr*P_co*(1-(mol_co2*mol_h2)/(mol_co*mol_h2o*Keq))/10^5;
rate3 = j/(2*F);
C34=(1/(rouPEN*cpPEN*hPEN))*[-dh3*rate3-j*vopt];
sigma=5.6704*10^(-8);
epsI=0.55
epsPEN=0.55
C35=(1/(rouPEN*cpPEN*hPEN))*(sigma/(1/epsI+1/epsPEN-1));
lmdI=22; %
rouI=8100;%
cpI=462;%
C41=lmdI/(rouI*cpI);
hI=0.0008;
KfI=KfPEN;
C42=KfI/(rouI*cpI*hI);
KaI=KaPEN;
C43=KaI/(rouI*cpI*hI);
C44=(1/(rouI*cpI*hI))*(sigma/(1/epsI+1/epsPEN-1));
sys(i)=-Af*(x(2*N+i)-u(3))/h-2*f_w*nd_fuel*x(i)/(rouf(i)*dhf^2);
sys(N+i)=-Af*(x(3*N+i)-u(4))/h-2*f_w*nd_oxid*x(N+i)/(roua(i)*dhf^2);
sys(2*N+i)=R*x(4*N+i)*((-(x(i)-u(1))/(Af*h)))/Mf;
%-(-0.5*0.016*rate3-0.044*rate3)/hf
sys(3*N+i)=R*x(5*N+i)*((-(x(N+i)-u(2))/(Af*h)))/Ma;
sys(4*N+i)=(-x(i)*x(4*N+i))/(rouf(i)*h*Af)+u(1)*u(5)/(rouf(i)*h*Af)+C11*(x(6*N+i)-x(4*N+i))+C12*(x(7*N+i)-x(4*N+i))+C11*(-dh1*rate1-rate2*dh2)/KfPEN;
sys(5*N+i)=(-x(N+i)*x(5*N+i))/(roua(i)*h*Af)+u(2)*u(6)/(roua(i)*h*Af)+C21*(x(6*N+i)-x(5*N+i))+C22*(x(7*N+i)-x(5*N+i));
sys(6*N+i)= C31*(x(6*N+i)-u(7))/h/h-C32*(x(6*N+i)-x(4*N+i))-C33*(x(6*N+i)-x(5*N+i))+C35*((x(7*N+i))^4-(x(6*N+i))^4)+C34;
sys(7*N+i)= C41*(x(7*N+i)-u(8))/h/h-C42*(x(7*N+i)-x(4*N+i))-C43*(x(7*N+i)-x(5*N+i))-C44*((x(7*N+i))^4-(x(6*N+i))^4);
sys(8*N+i)=(-x(i)*x(8*N+i)+u(1)*u(9))/(h*Af*rouf(i))-rate1/hf;
sys(9*N+i)=(-x(i)*x(9*N+i)+u(1)*u(10))/(h*Af*rouf(i))+(3*rate1+rate2-rate3)/hf;
sys(10*N+i)=(-x(i)*x(10*N+i)+u(1)*u(11))/(h*Af*rouf(i))-rate2/hf;
sys(11*N+i)=(-x(i)*x(11*N+i)+u(1)*u(12))/(h*Af*rouf(i))+(rate2+rate3)/hf;
sys(12*N+i)=(-x(i)*x(12*N+i)+u(1)*u(13))/(h*Af*rouf(i))-(rate1+rate2-rate3)/hf;
sys(13*N+i)=(-x(N+i)*x(13*N+i)+u(2)*u(14))/(h*Af*roua(i))-rate3/(2*hf);
sys(14*N+i)=(-x(N+i)*x(14*N+i)+u(2)*u(15))/(h*Af*roua(i));
sys(15*N+i)=(-x(N+i)*x(15*N+i)+u(2)*u(16))/(h*Af*roua(i))-rate3/(hf);
else if i==2
roua(i)=x(13*N+i-1)*0.032+x(14*N+i-1)*0.028+x(15*N+i-1)*0.044;
rouf(i)=x(8*N+i-1)*0.016+x(9*N+i-1)*0.002+x(10*N+i-1)*0.028+x(11*N+i-1)*0.044+x(12*N+i-1)*0.018;%
nd_ch4=-9.9989+529.37*(x(4*N+i-1)/1000)-543.82*(x(4*N+i-1)/1000).^2+548.11*(x(4*N+i-1)/1000).^3-367.06*(x(4*N+i-1)/1000).^4+140.48*(x(4*N+i-1)/1000).^5-22.92*(x(4*N+i-1)/1000).^6;
nd_co2=-20.434+680.07*(x(4*N+i-1)/1000)-432.49*(x(4*N+i-1)/1000).^2+244.22*(x(4*N+i-1)/1000).^3-85.929*(x(4*N+i-1)/1000).^4+14.45*(x(4*N+i-1)/1000).^5-0.4564*(x(4*N+i-1)/1000).^6;
nd_h2o=-6.7541+244.93*(x(4*N+i-1)/1000)+419.5*(x(4*N+i-1)/1000).^2-522.38*(x(4*N+i-1)/1000).^3+348.12*(x(4*N+i-1)/1000).^4-126.96*(x(4*N+i-1)/1000).^5+19.591*(x(4*N+i-1)/1000).^6;
nd_co=-4.9137+793.65*(x(4*N+i-1)/1000)+875.9*(x(4*N+i-1)/1000).^2+883.75*(x(4*N+i-1)/1000).^3-572.14*(x(4*N+i-1)/1000).^4+208.42*(x(4*N+i-1)/1000).^5-32.298*(x(4*N+i-1)/1000).^6;
nd_h2=15.553+299.78*(x(4*N+i-1)/1000)-244.34*(x(4*N+i-1)/1000).^2+249.41*(x(4*N+i-1)/1000).^3-167.51*(x(4*N+i-1)/1000).^4+62.966*(x(4*N+i-1)/1000).^5-9.9892*(x(4*N+i-1)/1000).^6;
mol_ch4=x(8*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_h2=x(9*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_co=x(10*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_co2=x(11*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_h2o=x(12*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
nd_fuel= (nd_co2*mol_co2*44^(0.5)+nd_h2o*mol_h2o*18^(0.5)+nd_ch4*mol_ch4*16^(0.5)+nd_co*mol_co*28^(0.5)+nd_h2*mol_h2*2^(0.5))*10^(-7)/(mol_co2*44^(0.5)+mol_h2o*18^(0.5)+mol_ch4*16^(0.5)+mol_co*28^(0.5)+mol_h2*2^(0.5));%
cp_co2=4.3669+204.6*(x(4*N+i-1)/1000)-471.33*(x(4*N+i-1)/1000).^2+657.88*(x(4*N+i-1)/1000).^3-519.9*(x(4*N+i-1)/1000).^4+214.58*(x(4*N+i-1)/1000).^5-35.992*(x(4*N+i-1)/1000).^6;
cp_h2o=37.737-41.205*(x(4*N+i-1)/1000)+146.01*(x(4*N+i-1)/1000).^2-217.08*(x(4*N+i-1)/1000).^3+181.54*(x(4*N+i-1)/1000).^4-79.409*(x(4*N+i-1)/1000).^5+14.105*(x(4*N+i-1)/1000).^6;
cp_co=30.429-8.1781*(x(4*N+i-1)/1000)+5.2062*(x(4*N+i-1)/1000).^2+41.974*(x(4*N+i-1)/1000).^3-66.346*(x(4*N+i-1)/1000).^4+37.456*(x(4*N+i-1)/1000).^5-7.6538*(x(4*N+i-1)/1000).^6;
cp_ch4=47.964-178.59*(x(4*N+i-1)/1000)+712.55*(x(4*N+i-1)/1000).^2-1068.7*(x(4*N+i-1)/1000).^3+856.93*(x(4*N+i-1)/1000).^4-358.75*(x(4*N+i-1)/1000).^5+61.321*(x(4*N+i-1)/1000).^6;
cp_h2=21.157+56.036*(x(4*N+i-1)/1000)-150.55*(x(4*N+i-1)/1000).^2+199.29*(x(4*N+i-1)/1000).^3-136.15*(x(4*N+i-1)/1000).^4+46.903*(x(4*N+i-1)/1000).^5-6.4725*(x(4*N+i-1)/1000).^6;
cp_mol=mol_ch4*cp_ch4+mol_h2*cp_h2+mol_co*cp_co+mol_co2*cp_co2+mol_h2o*cp_h2o;
cpf=1000*cp_mol/(mol_ch4*16+ mol_h2*2+mol_co*28+mol_co2*44+mol_h2o*18);
Mf=(16*mol_ch4+2*mol_h2+28*mol_co+44*mol_co2+18*mol_h2o)/1000;%fuel mol mass kg/mol
KfPEN=4.25;
C11=KfPEN/(rouf(i)*cpf*hf);
C12=C11;
nd_co2_ca=-20.434+680.07*(x(5*N+i-1)/1000)-432.49*(x(5*N+i-1)/1000).^2+244.22*(x(5*N+i-1)/1000).^3-85.929*(x(5*N+i-1)/1000).^4+14.45*(x(5*N+i-1)/1000).^5-0.4564*(x(5*N+i-1)/1000).^6;
nd_o2=-1.6918+889.75*(x(5*N+i-1)/1000)-892.79*(x(5*N+i-1)/1000).^2+905.98*(x(5*N+i-1)/1000).^3-598.36*(x(5*N+i-1)/1000).^4+221.64*(x(5*N+i-1)/1000).^5-34.754*(x(5*N+i-1)/1000).^6;
nd_n2=1.2719+771.45*(x(5*N+i-1)/1000)-809.2*(x(5*N+i-1)/1000).^2+832.47*(x(5*N+i-1)/1000).^3-553.93*(x(5*N+i-1)/1000).^4+206.15*(x(5*N+i-1)/1000).^5-32.43*(x(5*N+i-1)/1000).^6;
mol_co2_ca=x(15*N+i-1)/(x(13*N+i-1)+x(14*N+i-1)+x(15*N+i-1));
mol_o2=x(13*N+i-1)/(x(13*N+i-1)+x(14*N+i-1)+x(15*N+i-1));
mol_n2=x(14*N+i-1)/(x(13*N+i-1)+x(14*N+i-1)+x(15*N+i-1));
nd_oxid= (nd_co2_ca*mol_co2_ca*44^(0.5)+nd_o2*mol_o2*32^(0.5)+nd_n2*mol_n2*28^(0.5))*10^(-7)/(mol_co2_ca*44^(0.5)+mol_o2*32^(0.5)+mol_n2*28^(0.5));%
cp_co2=4.3669+204.6*(x(5*N+i-1)/1000)-471.33*(x(5*N+i-1)/1000).^2+657.88*(x(5*N+i-1)/1000).^3-519.9*(x(5*N+i-1)/1000).^4+214.58*(x(5*N+i-1)/1000).^5-35.992*(x(5*N+i-1)/1000).^6;
cp_o2=34.85-57.975*(x(5*N+i-1)/1000)+203.68*(x(5*N+i-1)/1000).^2-300.37*(x(5*N+i-1)/1000).^3+231.72*(x(5*N+i-1)/1000).^4-91.821*(x(5*N+i-1)/1000).^5+14.776*(x(5*N+i-1)/1000).^6;
cp_n2=29.027+4.8987*(x(5*N+i-1)/1000)-38.04*(x(5*N+i-1)/1000).^2+105.17*(x(5*N+i-1)/1000).^3-113.56*(x(5*N+i-1)/1000).^4+55.554*(x(5*N+i-1)/1000).^5-10.35*(x(5*N+i-1)/1000).^6;
cp_mol_ox=mol_o2*cp_o2+mol_n2*cp_n2+mol_co2_ca*cp_co2;
cpa=1000*cp_mol_ox/(mol_co2_ca*44+ mol_o2*32+mol_n2*28);
Ma=(44*mol_co2_ca+32*mol_o2+28*mol_n2)/1000;% oxidant mol mass
ha=0.0008;%氧化剂流道的高度
KaPEN=10.1;
C21=KaPEN/(roua(i)*cpa*ha);
C22=C21;
lmdPEN=22;
rouPEN=2000
cpPEN=800;
C31=lmdPEN/(rouPEN*cpPEN);
hPEN=0.001;%pen 厚度
C32=KfPEN/(rouPEN*cpPEN*hPEN);
C33=KaPEN/(rouPEN*cpPEN*hPEN);
dh1= (206e3 - 16379.71 + 8.314*(7.951*x(4*N+i-1) - 4.354*10^(-3)*x(4*N+i-1)^2 + 0.7213*10^(-6)*x(4*N+i-1)^3 - 0.097*10^5/x(4*N+i-1)));
dh2= (-41.1*10^3 - 7657.295 + 8.314*(1.86*x(4*N+i-1) - 0.27*10^(-3)*x(4*N+i-1)^2 + 1.164*10^5/x(4*N+i-1)));
dh3= (-241*10^3 - 4098.471 + 8.314*(-1.5985*x(6*N+i-1) + 0.3875*10^(-3)*x(6*N+i-1)^2 - 0.1515*10^5/x(6*N+i-1)));%
Ca=2.04*10^(-3);
C1=3.28*10^(-9);
C2=3.39*10^(-6);
Cir=1.12*10^(-2);
detaHc1=132000;
detaHc2=67100;
detaHa=23700;
detaHir=23000;
vopt=0.7;%
detaG=-242000+45.8*x(6*N+i-1);
E=-detaG/(2*F);
Ph2_an=x(2*N+i-1)*mol_h2*10^(-5);
Po2_ca=x(3*N+i-1)*mol_o2*10^(-5);
Pco2_ca=x(3*N+i-1)*mol_co2_ca*10^(-5);
Ph2o_an=x(2*N+i-1)*mol_h2o*10^(-5);
Pco2_an=x(2*N+i-1)*mol_co2*10^(-5);
Nerst=R*x(6*N+i-1)*log((Ph2_an*Pco2_ca*Po2_ca^0.5)/(Ph2o_an*Pco2_an))/(2*F);
%(Ph2_an*Pco2_ca*Po2_ca^0.5)/(Ph2o_an*Pco2_an)
Ran=Ca*(Ph2_an)^(-0.5)*exp(detaHa/(R*x(6*N+i-1)));
Rca=C1*(Po2_ca)^(-0.75)*(Pco2_ca)^(0.5)*exp(detaHc1/(R*x(6*N+i-1)))+C2*exp(detaHc2/(R*x(6*N+i-1)))*mol_co2_ca^(-1);
Rir=Cir*exp(detaHir/(R*x(6*N+i-1)));% change the pressure unit from pa to atm
Rtol=Ran+Rca+Rir;
j=10000*(E-vopt-Nerst)/Rtol;%current density unit is A/m2
k0=4274*10^5;%常数，unit is mol/s*m2*pa
Ea=82000;%j/mol
kwgsr=2.35;% water gas shift reaction
Keq=exp(4276/x(4*N+i-1)-3.961);
P_ch4=x(2*N+i-1)*mol_ch4;% ch4的分压力， pa
P_co=x(2*N+i-1)*mol_co;%co分压力，pa
rate1 = k0*P_ch4*exp(- Ea/(R*x(4*N+i-1))); rate2 = kwgsr*P_co*(1-(mol_co2*mol_h2)/(mol_co*mol_h2o*Keq))/10^5;
rate3 = j/(2*F);
C34=(1/(rouPEN*cpPEN*hPEN))*[-dh3*rate3*10^3-j*vopt];
sigma=5.6704*10^(-8); epsI=0.55
epsPEN=0.55;%
C35=(1/(rouPEN*cpPEN*hPEN))*(sigma/(1/epsI+1/epsPEN-1));
lmdI=22; %
rouI=8100;% cpI=462;% C41=lmdI/(rouI*cpI);
hI=0.0008; % KfI=KfPEN;
C42=KfI/(rouI*cpI*hI);
KaI=KaPEN;
C43=KaI/(rouI*cpI*hI);
C44=(1/(rouI*cpI*hI))*(sigma/(1/epsI+1/epsPEN-1));
sys(i)=-Af*(x(2*N+i)-x(2*N+i-1))/h-2*f_w*nd_fuel*x(i)/(rouf(i)*dhf^2);
sys(N+i)=-Af*(x(3*N+i)-x(3*N+i-1))/h-2*f_w*nd_oxid*x(N+i)/(roua(i)*dhf^2);
sys(2*N+i)=R*x(4*N+i)*(-(x(i)-x(i-1))/(Af*h))/Mf;
sys(3*N+i)=R*x(5*N+i)*(-(x(N+i)-x(N+i-1))/(Af*h))/Ma;
sys(4*N+i)=(-x(i)*x(4*N+i))/(rouf(i)*h*Af)+x(i-1)*x(4*N+i-1)/(rouf(i)*h*Af)+C11*(x(6*N+i)-x(4*N+i))+C12*(x(7*N+i)-x(4*N+i))+C11*(-dh1*rate1-dh2*rate2)/KfPEN;
sys(5*N+i)=(-x(N+i)*x(5*N+i))/(roua(i)*h*Af)+x(N+i-1)*x(5*N+i-1)/(roua(i)*h*Af)+C21*(x(6*N+i)-x(5*N+i))+C22*(x(7*N+i)-x(5*N+i));
sys(6*N+i)=C31*(x(6*N+i)-2*x(6*N+i-1)+u(7))/h/h-C32*(x(6*N+i)-x(4*N+i))-C33*(x(6*N+i)-x(5*N+i))+C35*((x(7*N+i))^4-(x(6*N+i))^4)+C34;
sys(7*N+i)=C41*(x(7*N+i)-2*x(7*N+i-1)+u(8))/h/h-C42*(x(7*N+i)-x(4*N+i))-C43*(x(7*N+i)-x(5*N+i))+C44*((x(7*N+i))^4-(x(6*N+i))^4);
sys(8*N+i)=(-x(i)*x(8*N+i)+x(i-1)*x(8*N+i-1))/(h*Af*rouf(i))-rate1/hf;
sys(9*N+i)=(-x(i)*x(9*N+i)+x(i-1)*x(9*N+i-1))/(h*Af*rouf(i))+(3*rate1+rate2-rate3)/hf;
sys(10*N+i)=(-x(i)*x(10*N+i)+x(i-1)*x(10*N+i-1))/(h*Af*rouf(i))-rate2/hf;
sys(11*N+i)=(-x(i)*x(11*N+i)+x(i-1)*x(11*N+i-1))/(h*Af*rouf(i))+(rate2+rate3)/hf;
sys(12*N+i)=(-x(i)*x(12*N+i)+x(i-1)*x(12*N+i-1))/(h*Af*rouf(i))-(rate1+rate2-rate3)/hf;
sys(13*N+i)=(-x(N+i)*x(13*N+i)+x(N+i-1)*x(13*N+i-1))/(h*Af*roua(i))-rate3/(2*hf);
sys(14*N+i)=(-x(N+i)*x(14*N+i)+x(N+i-1)*x(14*N+i-1))/(h*Af*roua(i));
sys(15*N+i)=(-x(N+i)*x(15*N+i)+x(N+i-1)*x(15*N+i-1))/(h*Af*roua(i))-rate3/(hf);
else
roua(i)=x(13*N+i-1)*0.032+x(14*N+i-1)*0.028+x(15*N+i-1)*0.044;
rouf(i)=x(8*N+i-1)*0.016+x(9*N+i-1)*0.002+x(10*N+i-1)*0.028+x(11*N+i-1)*0.044+x(12*N+i-1)*0.018;% nd_ch4=-9.9989+529.37*(x(4*N+i-1)/1000)-543.82*(x(4*N+i-1)/1000).^2+548.11*(x(4*N+i-1)/1000).^3-367.06*(x(4*N+i-1)/1000).^4+140.48*(x(4*N+i-1)/1000).^5-22.92*(x(4*N+i-1)/1000).^6;
nd_co2=-20.434+680.07*(x(4*N+i-1)/1000)-432.49*(x(4*N+i-1)/1000).^2+244.22*(x(4*N+i-1)/1000).^3-85.929*(x(4*N+i-1)/1000).^4+14.45*(x(4*N+i-1)/1000).^5-0.4564*(x(4*N+i-1)/1000).^6;
nd_h2o=-6.7541+244.93*(x(4*N+i-1)/1000)+419.5*(x(4*N+i-1)/1000).^2-522.38*(x(4*N+i-1)/1000).^3+348.12*(x(4*N+i-1)/1000).^4-126.96*(x(4*N+i-1)/1000).^5+19.591*(x(4*N+i-1)/1000).^6;
nd_co=-4.9137+793.65*(x(4*N+i-1)/1000)+875.9*(x(4*N+i-1)/1000).^2+883.75*(x(4*N+i-1)/1000).^3-572.14*(x(4*N+i-1)/1000).^4+208.42*(x(4*N+i-1)/1000).^5-32.298*(x(4*N+i-1)/1000).^6;
nd_h2=15.553+299.78*(x(4*N+i-1)/1000)-244.34*(x(4*N+i-1)/1000).^2+249.41*(x(4*N+i-1)/1000).^3-167.51*(x(4*N+i-1)/1000).^4+62.966*(x(4*N+i-1)/1000).^5-9.9892*(x(4*N+i-1)/1000).^6;
mol_ch4=x(8*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_h2=x(9*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_co=x(10*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_co2=x(11*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
mol_h2o=x(12*N+i-1)/(x(8*N+i-1)+x(9*N+i-1)+x(10*N+i-1)+x(11*N+i-1)+x(12*N+i-1));
nd_fuel= (nd_co2*mol_co2*44^(0.5)+nd_h2o*mol_h2o*18^(0.5)+nd_ch4*mol_ch4*16^(0.5)+nd_co*mol_co*28^(0.5)+nd_h2*mol_h2*2^(0.5))*10^(-7)/(mol_co2*44^(0.5)+mol_h2o*18^(0.5)+mol_ch4*16^(0.5)+mol_co*28^(0.5)+mol_h2*2^(0.5));%
cp_co2=4.3669+204.6*(x(4*N+i-1)/1000)-471.33*(x(4*N+i-1)/1000).^2+657.88*(x(4*N+i-1)/1000).^3-519.9*(x(4*N+i-1)/1000).^4+214.58*(x(4*N+i-1)/1000).^5-35.992*(x(4*N+i-1)/1000).^6;
cp_h2o=37.737-41.205*(x(4*N+i-1)/1000)+146.01*(x(4*N+i-1)/1000).^2-217.08*(x(4*N+i-1)/1000).^3+181.54*(x(4*N+i-1)/1000).^4-79.409*(x(4*N+i-1)/1000).^5+14.105*(x(4*N+i-1)/1000).^6;
cp_co=30.429-8.1781*(x(4*N+i-1)/1000)+5.2062*(x(4*N+i-1)/1000).^2+41.974*(x(4*N+i-1)/1000).^3-66.346*(x(4*N+i-1)/1000).^4+37.456*(x(4*N+i-1)/1000).^5-7.6538*(x(4*N+i-1)/1000).^6;
cp_ch4=47.964-178.59*(x(4*N+i-1)/1000)+712.55*(x(4*N+i-1)/1000).^2-1068.7*(x(4*N+i-1)/1000).^3+856.93*(x(4*N+i-1)/1000).^4-358.75*(x(4*N+i-1)/1000).^5+61.321*(x(4*N+i-1)/1000).^6;
cp_h2=21.157+56.036*(x(4*N+i-1)/1000)-150.55*(x(4*N+i-1)/1000).^2+199.29*(x(4*N+i-1)/1000).^3-136.15*(x(4*N+i-1)/1000).^4+46.903*(x(4*N+i-1)/1000).^5-6.4725*(x(4*N+i-1)/1000).^6;
cp_mol=mol_ch4*cp_ch4+mol_h2*cp_h2+mol_co*cp_co+mol_co2*cp_co2+mol_h2o*cp_h2o;
cpf=1000*cp_mol/(mol_ch4*16+ mol_h2*2+mol_co*28+mol_co2*44+mol_h2o*18);
Mf=(16*mol_ch4+2*mol_h2+28*mol_co+44*mol_co2+18*mol_h2o)/1000;%KfPEN=4.25;
C11=KfPEN/(rouf(i)*cpf*hf);
C12=C11;
KaPEN=10.1;
C21=KaPEN/(roua(i)*cpa*ha);
C22=C21;
lmdPEN=22; rouPEN=2000 cpPEN=800; C31=lmdPEN/(rouPEN*cpPEN);
hPEN=0.001 C32=KfPEN/(rouPEN*cpPEN*hPEN);
C33=KaPEN/(rouPEN*cpPEN*hPEN);
dh1= (206*10^3 - 16379.71 + 8.314*(7.951*x(4*N+i-1) - 4.354*10^(-3)*x(4*N+i-1)^2 + 0.7213*10^(-6)*x(4*N+i-1)^3 - 0.097*10^5/x(4*N+i-1)));
dh2= (-41.1*10^3 - 7657.295 + 8.314*(1.86*x(4*N+i-1) - 0.27*10^(-3)*x(4*N+i-1)^2 + 1.164*10^5/x(4*N+i-1)));
dh3= (-241*10^3 - 4098.471 + 8.314*(-1.5985*x(6*N+i-1) + 0.3875*10^(-3)*x(6*N+i-1)^2 - 0.1515*10^5/x(6*N+i-1
Ca=2.04*10^(-3);
C1=3.28*10^(-9);
C2=3.39*10^(-6);
Cir=1.12*10^(-2);
detaHc1=132000;
detaHc2=67100;
detaHa=23700;
detaHir=23000;
vopt=0.7
detaG=-242000+45.8*x(6*N+i-1);
E=-detaG/(2*F);
Ph2_an=x(2*N+i-1)*mol_h2*10^(-5);%
Po2_ca=x(3*N+i-1)*mol_o2*10^(-5);
Pco2_ca=x(3*N+i-1)*mol_co2_ca*10^(-5);
Ph2o_an=x(2*N+i-1)*mol_h2o*10^(-5);
Pco2_an=x(2*N+i-1)*mol_co2*10^(-5);
Nerst=R*x(4*N+i-1)*log((Ph2_an*Pco2_ca*Po2_ca^0.5)/(Ph2o_an*Pco2_an))/(2*F);
(Ph2_an*Pco2_ca*Po2_ca^0.5)/(Ph2o_an*Pco2_an)
Ran=Ca*(Ph2_an)^(-0.5)*exp(detaHa/(R*x(6*N+i-1)));
Rca=C1*(Po2_ca)^(-0.75)*(Pco2_ca)^(0.5)*exp(detaHc1/(R*x(6*N+i-1)))+C2*exp(detaHc2/(R*x(6*N+i-1)))*mol_co2_ca^(-1);
Rir=Cir*exp(detaHir/(R*x(6*N+i-1)));% change the pressure unit from pa to atm
Rtol=Ran+Rca+Rir;
j=10000*(E-vopt-Nerst)/Rtol;%current density unit is A/m2
k0=4274*10^5
Ea=82000 kwgsr=2.35;%
Keq=exp(4276/x(4*N+i-1)-3.961);
P_ch4=x(2*N+i-1)*mol_ch4
P_co=x(2*N+i-1)*mol_co
rate1 = k0*P_ch4*exp(- Ea/(R*x(4*N+i-1))); rate2 = kwgsr*P_co*(1-(mol_co2*mol_h2)/(mol_co*mol_h2o*Keq))/10^5;
rate3 = j/(2*F);
C34=(1/(rouPEN*cpPEN*hPEN))*[-dh3*rate3-j*vopt];
sigma=5.6704*10^(-8); % epsI=0.55;%
epsPEN=0.55;% C35=(1/(rouPEN*cpPEN*hPEN))*(sigma/(1/epsI+1/epsPEN-1));
lmdI=22; rouI=8100 cpI=462 C41=lmdI/(rouI*cpI);
hI=0.0008;
KfI=KfPEN;
C42=KfI/(rouI*cpI*hI);
KaI=KaPEN;
C43=KaI/(rouI*cpI*hI);
C44=(1/(rouI*cpI*hI))*(sigma/(1/epsI+1/epsPEN-1));
sys(i)=-Af*(x(2*N+i)-x(2*N+i-1))/h-2*f_w*nd_fuel*x(i)/(rouf(i)*dhf^2);
sys(N+i)=-Af*(x(3*N+i)-x(3*N+i-1))/h-2*f_w*nd_oxid*x(N+i)/(roua(i)*dhf^2);
sys(2*N+i)=R*x(4*N+i)*(-(x(i)-x(i-1))/(Af*h))/Mf;
sys(3*N+i)=R*x(5*N+i)*(-(x(N+i)-x(N+i-1))/(Af*h))/Ma;
sys(4*N+i)=(-x(i)*x(4*N+i))/(rouf(i)*h*Af)+x(i-1)*x(4*N+i-1)/(rouf(i)*h*Af)+C11*(x(6*N+i)-x(4*N+i))+C12*(x(7*N+i)-x(4*N+i))+C11*(-dh1*rate1-dh2*rate2)/KfPEN;
sys(5*N+i)=(-x(N+i)*x(5*N+i))/(roua(i)*h*Af)+x(N+i-1)*x(5*N+i-1)/(roua(i)*h*Af)+C21*(x(6*N+i)-x(5*N+i))+C22*(x(7*N+i)-x(5*N+i));
sys(6*N+i)=C31*(x(6*N+i)-2*x(6*N+i-1)+x(6*N+i-2))/h/h-C32*(x(6*N+i)-x(4*N+i))-C33*(x(6*N+i)-x(5*N+i))+C35*((x(7*N+i))^4-(x(6*N+i))^4)+C34;
sys(7*N+i)=C41*(x(7*N+i)-2*x(7*N+i-1)+x(7*N+i-2))/h/h-C42*(x(7*N+i)-x(4*N+i))-C43*(x(7*N+i)-x(5*N+i))+C44*((x(7*N+i))^4-(x(6*N+i))^4);
sys(8*N+i)=(-x(i)*x(8*N+i)+x(i-1)*x(8*N+i-1))/(h*Af*rouf(i))-rate1/hf;
sys(9*N+i)=(-x(i)*x(9*N+i)+x(i-1)*x(9*N+i-1))/(h*Af*rouf(i))+(3*rate1+rate2-rate3)/hf;
sys(10*N+i)=(-x(i)*x(10*N+i)+x(i-1)*x(10*N+i-1))/(h*Af*rouf(i))-rate2/hf;
sys(11*N+i)=(-x(i)*x(11*N+i)+x(i-1)*x(11*N+i-1))/(h*Af*rouf(i))+(rate2+rate3)/hf;
sys(12*N+i)=(-x(i)*x(12*N+i)+x(i-1)*x(12*N+i-1))/(h*Af*rouf(i))-(rate1+rate2-rate3)/hf;
sys(13*N+i)=(-x(N+i)*x(13*N+i)+x(N+i-1)*x(13*N+i-1))/(h*Af*roua(i))-rate3/(2*hf);
sys(14*N+i)=(-x(N+i)*x(14*N+i)+x(N+i-1)*x(14*N+i-1))/(h*Af*roua(i));
sys(15*N+i)=(-x(N+i)*x(15*N+i)+x(N+i-1)*x(15*N+i-1))/(h*Af*roua(i))-rate3/(hf);
end
end
end 这个太专业了啊！
页:
[1]