声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2461|回复: 1

[综合讨论] 使用s函数时提示call must be a real vector

[复制链接]
发表于 2008-12-7 18:02 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
提示: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
回复
分享到:

使用道具 举报

发表于 2010-12-7 10:49 | 显示全部楼层
这个太专业了啊!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-28 20:18 , Processed in 0.064320 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表