Temperature profile must increase expoenetial but mine is decreasing,dy(15) is temperature profile diff equation.If ra or rp in code gets changed it changes because rest is constant
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
An adiabtic bathc reactor model is getting developed the function is given here,I have to get an exponential increase when plotting time vs temperature dy(15),but i am getting exponential decrease.Please help me with that
function dy=adiabatic_batch(t,y)
k=deadkineticestimation(y(15));
dy=zeros(17,1);
if(y(6)==0.0 || y(7)==0.0 || y(9)==0.0 || y(10)==0.0 || y(12)==0.0 || y(13)==0.0)
lam3=0.0;
mu3=0.0;
gam3=0.0;
else
lam3=(y(8)*((2*y(8)*y(6))-(y(7)*y(7))))./(y(7)*y(6));
mu3=(y(11)*((2*y(11)*y(9))-(y(10)*y(10))))./(y(10)*y(9));
gam3=(y(14)*((2*y(14)*y(12))-(y(13)*y(13))))./(y(13)*y(12));
end
Hra=50.100;%%kj/mol
Hrp=-23.300;
ra=-((-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(3)*y(1)*y(9))+(k(4)*y(6)*y(4)));
rp=((-1*k(5)*y(5)*y(3))-(k(7)*y(5)*y(6))+(k(8)*y(6)));
mtot=11.734;
cp=2.06;
Mmon=144.14;
Mcat=405.1;
Mcocat=186.34;
Macid=144.22;
Di=315.9; %in mm
%equation 16
density_cat=1200;
density_cocat=820;
density_acid=908.8;
density_pla=k(17);
% Dilution terms included
rcat=-ra;
rcocat=(-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(9)*y(2)*y(6))+(k(10)*y(3)*y(9));
racid=ra;
rmon=-rp;
rini=(-1*k(5)*y(5)*y(3))+(k(1)*y(1)*y(2))-(k(2)*y(3)*y(4))+(k(9)*y(2)*y(6))-(k(10)*y(9)*y(3));
rla0=(k(3)*y(9)*y(1))-(k(4)*y(6)*y(4))+(k(5)*y(5)*y(3))-(k(9)*y(6)*y(2))+(k(10)*y(9)*y(3));
rla1=((k(3)*y(10)*y(1))-(k(4)*y(7)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*y(6))-(k(8)*y(6))-(k(9)*y(7)*y(2))+ ...
(k(10)*y(10)*y(3))-(k(11)*y(7)*y(9))+(k(12)*y(10)*y(6))-(k(13)*y(7)*(y(10)-y(9)))+(0.5*k(14)*y(6)*(y(11)-y(10)))-(k(15)*y(7)*((y(13)-y(12)))+(0.5*k(16)*y(6)*(y(14)-y(13)))-(0.5*k(20)*(y(8)-y(7)))));
rla2=(k(3)*y(11)*y(1))-(k(4)*y(8)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*((2*y(7))+y(6)))+(k(8)*(y(6)-(2*y(7))))- ...
(k(9)*y(8)*y(2))+(k(10)*y(11)*y(3))-(k(11)*y(8)*y(9))+(k(12)*y(11)*y(6))+(0.33333*k(13)*y(6)*(y(7)-lam3))+ ...
(k(14)*y(7)*(y(8)-y(7)))-(k(15)*y(8)*(y(10)-y(9)))+(0.16667*k(16)*y(6)*((2*mu3))-((3*y(11))+y(10)))-((k(23)*y(8)*(y(13)-y(12))))+(0.166*k(24)*y(6)*((2*gam3)-(3*y(14))+y(13)))-(0.1666*k(20)*((4*lam3)-(3*y(8))-y(7)));
rmu0=(-1*k(3)*y(9)*y(1))+(k(4)*y(6)*y(4))+(k(9)*y(6)*y(2))-(k(10)*y(9)*y(3));
rmu1=(-1*k(3)*y(10)*y(1))+(k(4)*y(7)*y(4))+(k(9)*y(7)*y(2))-(k(10)*y(10)*y(3))+(k(11)*y(7)*y(9))-(k(12)*y(10)*y(6))...
+(k(15)*y(7)*(y(10)-y(9)))-(0.5*k(16)*y(6)*(y(11)-y(10)))-(0.5*k(20)*(y(11)-y(10)));
rmu2=(-1*k(3)*y(11)*y(1))+(k(4)*y(8)*y(4))+(k(9)*y(8)*y(2))-(k(10)*y(11)*y(3))+(k(11)*y(8)*y(9))-(k(12)*y(11)*y(6))+ ...
(k(15)*y(8)*(y(10)-y(9)))+(k(16)*y(7)*(y(11)-y(10)))+(0.16667*k(16)*y(6)*((-4*mu3)+(3*y(11))+y(10)))-(0.166*k(20)*((4*mu3)-(3*y(11))-y(10)));
rga0=((k(20)*(y(7)-y(6)))+(k(21)*(y(10)-y(9))));
rga1=(k(13)*y(7)*(y(13)-y(12))-(0.5*k(14)*y(6)*(y(14)-y(13)))-(k(20)*(y(14)-y(13)))...
+(0.5*k(21)*(y(8)-y(7)))+(0.5*k(22)*(y(11)-y(10))));
rga2=((k(13)*y(8)*(y(13)-y(12))+(k(14)*y(7)*(y(14)-y(13)))+(0.16666*k(15)*y(6)*((-4*gam3+(3*y(14)+y(13)))))-(0.333*k(20)*((4*gam3)-(3*y(14))-y(13))))+(0.1666*k(21)*((2*lam3)-(3*y(8))+y(7)))...
+(0.1666*k(22)*((2*mu3)-(3*y(11))+y(10))));
% dilution terms for 1-14
dilution_1=rcat-(y(1)/y(16))*dy(16);
dilution_2=rcocat-(y(2)/y(16))*dy(16);
dilution_3=rini-(y(3)/y(16))*dy(16);
dilution_4=racid-(y(4)/y(16))*dy(16);
dilution_5=rmon-(y(5)/y(16))*dy(16);
dilution_6=rla0-(y(6)/y(16))*dy(16);
dilution_7=rla1-(y(7)/y(16))*dy(16);
dilution_8=rla2-(y(8)/y(16))*dy(16);
dilution_9=rmu0-(y(9)/y(16))*dy(16);
dilution_10=rmu1-(y(10)/y(16))*dy(16);
dilution_11=rmu2-(y(11)/y(16))*dy(16);
dilution_12=rga0-(y(12)/y(16))*dy(16);
dilution_13=rga1-(y(13)/y(16))*dy(16);
dilution_14=rga2-(y(14)/y(16))*dy(16);
%}
a=rmon/k(17);
b=(y(5)*y(16)*Mmon)/((k(17))^2);
diff_mon=((-0.8462)/(1.110865+(7.391*10^-4*y(15)))^2);
diff_pla=((-0.8462)/(1.110865+(7.391*10^-4*y(15)))^2);
c=rcat/density_cat;
d=rcocat/density_cocat;
e=racid/density_acid;
f=(1-(1/density_pla));
g=(mtot-(((y(5)*Mmon)+(y(1)*Mcat)+(y(2)*Mcocat)+(y(4)*Macid))*y(16))/((density_pla)^2));
dy(1)=(-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(3)*y(1)*y(9))+(k(4)*y(6)*y(4))+dilution_1;
dy(2)=(-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(9)*y(2)*y(6))+(k(10)*y(3)*y(9))+dilution_2;
dy(3)=(-1*k(5)*y(5)*y(3))+(k(1)*y(1)*y(2))-(k(2)*y(3)*y(4))+(k(9)*y(2)*y(6))-(k(10)*y(9)*y(3))+dilution_3; % +(k(4)*y(6));
dy(4)=(k(1)*y(1)*y(2))-(k(2)*y(3)*y(4))+(k(3)*y(1)*y(9))-(k(4)*y(6)*y(4))+dilution_4;
dy(5)=(-1*k(5)*y(5)*y(3))-(k(7)*y(5)*y(6))+(k(8)*y(6))+dilution_5;
dy(6)=(k(3)*y(9)*y(1))-(k(4)*y(6)*y(4))+(k(5)*y(5)*y(3))-(k(9)*y(6)*y(2))+(k(10)*y(9)*y(3))+dilution_6;
dy(7)=((k(3)*y(10)*y(1))-(k(4)*y(7)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*y(6))-(k(8)*y(6))-(k(9)*y(7)*y(2))+ ...
(k(10)*y(10)*y(3))-(k(11)*y(7)*y(9))+(k(12)*y(10)*y(6))-(k(13)*y(7)*(y(10)-y(9)))+(0.5*k(14)*y(6)*(y(11)-y(10)))-(k(15)*y(7)*((y(13)-y(12)))+(0.5*k(16)*y(6)*(y(14)-y(13)))-(0.5*k(20)*(y(8)-y(7)))))+dilution_7;
%mu3
dy(8)=(k(3)*y(11)*y(1))-(k(4)*y(8)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*((2*y(7))+y(6)))+(k(8)*(y(6)-(2*y(7))))- ...
(k(9)*y(8)*y(2))+(k(10)*y(11)*y(3))-(k(11)*y(8)*y(9))+(k(12)*y(11)*y(6))+(0.33333*k(13)*y(6)*(y(7)-lam3))+ ...
(k(14)*y(7)*(y(8)-y(7)))-(k(15)*y(8)*(y(10)-y(9)))+(0.16667*k(16)*y(6)*((2*mu3))-((3*y(11))+y(10)))-((k(23)*y(8)*(y(13)-y(12))))+(0.166*k(24)*y(6)*((2*gam3)-(3*y(14))+y(13)))-(0.1666*k(20)*((4*lam3)-(3*y(8))-y(7)))+dilution_8;
dy(9)=(-1*k(3)*y(9)*y(1))+(k(4)*y(6)*y(4))+(k(9)*y(6)*y(2))-(k(10)*y(9)*y(3))+dilution_9;
dy(10)=(-1*k(3)*y(10)*y(1))+(k(4)*y(7)*y(4))+(k(9)*y(7)*y(2))-(k(10)*y(10)*y(3))+(k(11)*y(7)*y(9))-(k(12)*y(10)*y(6))...
+(k(15)*y(7)*(y(10)-y(9)))-(0.5*k(16)*y(6)*(y(11)-y(10)))-(0.5*k(20)*(y(11)-y(10)))+dilution_10;
dy(11)=(-1*k(3)*y(11)*y(1))+(k(4)*y(8)*y(4))+(k(9)*y(8)*y(2))-(k(10)*y(11)*y(3))+(k(11)*y(8)*y(9))-(k(12)*y(11)*y(6))+ ...
(k(15)*y(8)*(y(10)-y(9)))+(k(16)*y(7)*(y(11)-y(10)))+(0.16667*k(16)*y(6)*((-4*mu3)+(3*y(11))+y(10)))-(0.166*k(20)*((4*mu3)-(3*y(11))-y(10)))+dilution_11;
dy(12)=((k(20)*(y(7)-y(6)))+(k(21)*(y(10)-y(9)))+dilution_12);
dy(13)=(k(13)*y(7)*(y(13)-y(12))-(0.5*k(14)*y(6)*(y(14)-y(13)))-(k(20)*(y(14)-y(13)))...
+(0.5*k(21)*(y(8)-y(7)))+(0.5*k(22)*(y(11)-y(10)))+dilution_13);
dy(14)=((k(13)*y(8)*(y(13)-y(12))+(k(14)*y(7)*(y(14)-y(13)))+(0.16666*k(15)*y(6)*((-4*gam3+(3*y(14)+y(13)))))-(0.333*k(20)*((4*gam3)-(3*y(14))-y(13))))+(0.1666*k(21)*((2*lam3)-(3*y(8))+y(7)))...
+(0.1666*k(22)*((2*mu3)-(3*y(11))+y(10)))+dilution_14);
dy(15)=(((-Hra*ra+(-Hrp)*rp)*y(16)))/(mtot*cp);
dy(16)=(((a-(b*diff_mon*dy(15)+c+d+e)*f)-(g*diff_pla*dy(15))));
dy(17)=((4/(3.14*Di^2))*dy(16))*10^6;
end
calling Code:
y0=[1 15 0 0 6000 0 0 0 0 0 0 0 0 0 140 10.121 185 ];
ts=[0 2];
[P,X]=ode15s(@adiabatic_batch,ts,y0);
mw=144*(X(:,8)+X(:,11)+X(:,14))./(X(:,10)+(X(:,7)+X(:,13)));
mn=144*(X(:,7)+X(:,10)+X(:,13))./(X(:,9)+X(:,6)+X(:,12));
conv=1-(X(:,5)/y0(5));
pdi=mw./mn ;
Temp=X(:,15);
plot(P,Temp,'-','LineWidth',2,'DisplayName','3771');
hold on
2 commentaires
Rik
le 20 Déc 2020
This time I edited your question for you. Next time, please use the tools explained on this page to make your question more readable.
Réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!