DAEs (ODE 15s)
1 view (last 30 days)
Show older comments
I have a system of 10 differential equations and 3 algebraic equations.
I used ODE15s to integrate. But, I am not sure why the initial value of the variables defined by the algebraic equations changes in the first calculation and it is not the same as initial condition vector x0.( x(11,x(12) & x(13)) . (mfiles dynamosofcpub2RApril7 and sofcmultipdynpub2Rapril7)
I call this function [t x]=dynamosofcpub2RApril7 in command window. Do you have any suggestions? Even if I use "if" statement in dynamosofcpub2RApril7 to apply the step change in variable u at a specific time, the sysyem does not converge.
Please find the mfiles below.
function [t x]=dynamosofcpub2RApril7
global u
u=10;
M=eye(13);
M(13,13)=0;
M(12,12)=0;
M(11,11)=0;
options=odeset('Mass',M,'MassSingular','yes','RelTol',1e-6,'AbsTol',1e-6);
tspan=[0:0.1:2000];
x0=[951.314293356128,19.9949642189062,12.0282577584507,2.40556847002325,
11427.3505058679,90.2350929738687,10.0833309517251,12.0659722222222,
10.8536790016920,11470.2069812040,0.126623256318938,0.276095749286402,0.139558206138141];
[t x]=ode15s(@sofcmultipdynpub2Rapril7,tspan,x0,options);
end
function f=sofcmultipdynpub2Rapril7(t,x)
% constants
global u
f=zeros(13,1);
% Rload=5;
% cpaircat=34.4;
% cpfuelan=31.4;
% if t<=100
% % Rload=0.014;
% Rload=5;
% else
% Rload=10;
% end;
Rload=u;
dA=10^-4;
dC=0.5*10^-4;
dE=1.8*10^-4;
L=0.04;
B=0.04;
F=96485.3;
naircatin=3.8*10^-2;
nfuelanin=1.39*10^-3;
HA=10^-3;
HC=10^-3;
Taircatin=950;
Tfuelanin=950;
% Taircatin=1100;
% Tfuelanin=1100;
vaircatin=78.98;
vfuelanin=2.88;
rhocp=10^6;
deltaHR=-241830;
gammaA=5.7*10^7;
gammaC=7*10^9;
thetaaC=1.4;
thetacC=0.6;
thetaaA=2;
thetacA=1;
betha1=3.34*10^4;
betha2=1.03*10^4;
EA=140000;
EC=160000;
alpha=25;
Deltaanode=5*10^-5;
Deltacathode=5*10^-5;
cpH2=30;
cpH2o=44;
cpo2=36;
epsiloncat=5*10^-1;
thoucat=3;
epsilonan=5*10^-1;
thouan=3;
deltaan=5*10^-7;
deltacat=5*10^-7;
nuH2=7.07;
nuo2=16.6;
nuN2=17.9;
nuH2o=12.7;
R=8.314;
Mo2=31.994;
MH2=2.02;
MH2o=18.02;
MN2=28.02;
naircatin=3.8*10^-2;
nfuelanin=1.39*10^-3;
yo2in=0.2;
yN2in=0.8;
vaircat=vaircatin;
vfuelan=vfuelanin;
ao2=34.57;
bo2=0.1078e-2;
do2=-784586;
ah2o=34.36;
bh2o=0.0627e-2;
ch2o=0.56012e-5;
ah2=27.67;
bh2=0.3386e-2;
an2=27.17;
bn2=0.4180e-2;
yH2in=0.9;
yH2oin=0.1;
deltaGR=-175933;
deltaSR=-57;
Tref=1300;
%physical parameters variations
Do2N2=1.013*10^-7*(x(1))^1.75*(1/Mo2+1/MN2)^0.5/((x(2)*x(1)*9.86923267*10^-6+x(5)*R*9.86923267*10^-6*(1-x(4)/x(3)))*(nuo2^(1/3)+nuN2^(1/3))^2);
Dko2=97*deltacat*abs(sqrt(x(1)/Mo2));
Do2=1/(1/Do2N2+1/Dko2);
ko2cat=(epsiloncat/(thoucat*Deltacathode))*Do2;
Ho2catin=ao2*(Taircatin-298.15)+bo2*(Taircatin^2/2-298.15^2/2)-do2*(1/Taircatin-1/298.15);
HN2catin=an2*(Taircatin-298.15)+bn2*(Taircatin^2/2-298.15^2/2);
Haircatin=yo2in*Ho2catin+yN2in*HN2catin;
Ho2cat=ao2*(x(5)/x(3)-298.15)+bo2*((x(5)/x(3))^2/2-298.15^2/2)-do2*(1/(x(5)/x(3))-1/298.15);
HN2cat=an2*(x(5)/x(3)-298.15)+bn2*((x(5)/x(3))^2/2-298.15^2/2);
Haircat=x(4)/x(3)*Ho2cat+(1-x(4)/x(3))*HN2cat;
DkH2=97*deltaan*abs(sqrt(x(1)/MH2));
DH2H2o=1.013*10^-7*x(1)^1.75*(1/MH2+1/MH2o)^0.5/((x(6)*x(1)*9.86923267*10^-6+x(7)*x(1)*9.86923267*10^-6)*(nuH2^(1/3)+nuH2o^(1/3))^2);
DH2=1/(1/DH2H2o+1/DkH2);
kH2an=(epsilonan/(thouan*Deltaanode))*DH2;
HH2anin=ah2*(Tfuelanin-298.15)+bh2*(Tfuelanin^2/2-298.15^2/2);
HH2oanin=ah2o*(Tfuelanin-298.15)+bh2o*(Tfuelanin^2/2-298.15^2/2)+ch2o*(Tfuelanin^3/3-298.15^3/3)-241814;
Hfuelanin=yH2in*HH2anin+yH2oin*HH2oanin;
HH2an=ah2*(x(10)/x(8)-298.15)+bh2*((x(10)/x(8))^2/2-298.15^2/2);
HH2oan=ah2o*(x(10)/x(8)-298.15)+bh2o*((x(10)/x(8))^2/2-298.15^2/2)+ch2o*((x(10)/x(8))^3/3-298.15^3/3)-241814;
Hfuelan=x(9)/x(8)*HH2an+(1-x(9)/x(8))*HH2oan;
DkH2o=97*deltaan*abs(sqrt(x(1)/MH2o));
DH2oH2=DH2H2o;
DH2o=1/(1/DH2oH2+1/DkH2o);
kH2oan=(epsilonan/(thouan*Deltaanode))*DH2o;
%Variables definition
No2=ko2cat*(x(4)-x(2)/(R));
NH2=kH2an*(x(9)-x(6)/(R));
NH2o=kH2oan*(x(7)/(R)-x(8)*(1-x(9)/x(8)));
U0=-1/(2*F)*(deltaGR-deltaSR*(x(1)-Tref)+R*x(1)*log((x(7)*x(1))/(x(6)*x(1)*((x(2)*x(1))/(1.013*10^5))^0.5)));
rhoE=1/betha1*exp(betha2/x(1));
% algebraic equations
v=U0-x(12)-x(13)-rhoE*dE*x(11)/(L*B);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%physical properties of the solid
cph2solid=ah2+bh2*x(1);
cpo2solid=ao2+bo2*x(1)+do2*x(1)^(-2);
cpo2cat=ao2+bo2*x(1)+do2*x(1)^(-2);
cpn2cat=an2+bn2*x(5)/x(3);
cpaircat=x(4)/x(3).*cpo2cat+(1-x(4)/x(3)).*cpn2cat;
cph2an=ah2+bh2*x(10)/x(8);
cph2oan=ah2o+bh2o*x(10)/x(8)+ch2o*(x(10)/x(8))^2;
cpfuelan=x(9)/x(8)*cph2an+(1-x(9)/x(8))*cph2oan;
% cpfuelan=31.4;
%Differential equations
f(1)=1/(rhocp*(dA+dC+dE))*((-deltaHR/(2*F)-v)*x(11)/(L*B)+(alpha+cph2solid/(2*F)*x(11)/(L*B))*(x(10)/x(8)-x(1))+(alpha+cpo2solid/(4*F)*x(11)/(L*B))*(x(5)/x(3)-x(1)));
f(2)=R/Deltacathode*(No2-x(11)/(L*B*4*F));
f(3)=1/(L*B*HC)*(naircatin-vaircat*x(3)*B*HC-No2*L*B);
f(4)=1/(L*B*HC)*(naircatin*yo2in-vaircat*x(4)*B*HC-No2*L*B);
f(5)=1/(L*B*HC*cpaircat)*(naircatin*Haircatin-vaircat*x(3)*Haircat*B*HC+L*B*alpha*(x(1)-x(5)/x(3))-No2*Ho2cat*L*B);
f(6)=R/Deltaanode*(NH2-(1/(L*B))*(x(11)/(2*F)));
f(7)=R/Deltaanode*(-NH2o+(1/(L*B))*(x(11)/(2*F)));
f(8)=1/(L*B*HA)*(nfuelanin-vfuelan*x(8)*B*HA+L*B*(NH2o-NH2));
f(9)=1/(L*B*HA)*(nfuelanin*yH2in-vfuelan*x(9)*B*HA-NH2*L*B);
f(10)=1/(L*B*HA*cpfuelan)*(nfuelanin*Hfuelanin-vfuelan*x(8)*Hfuelan*B*HA+alpha*(x(1)-x(10)/x(8))*L*B+L*B*(HH2oan*NH2o-HH2an*NH2));
f(11)=x(11)/(L*B)-gammaA*((x(6)*x(1))/(1.013*10^5))*((x(7)*x(1))/(1.013*10^5))*exp(-EA/(R*x(1)))*(exp(thetaaA*F/(R*x(1))*x(12))-exp(-thetacA*F/(R*x(1))*x(12)));
f(12)=x(11)/(L*B)-gammaC*((x(2)*x(1))/(1.013*10^5)).^0.25*exp(-EC/(R*x(1)))*(exp(thetaaC*F*x(13)/(R*x(1)))-exp(-thetacC*F*x(13)/(R*x(1))));
f(13)=v-Rload*x(11);
end
2 Comments
Rafael Freire
on 11 Apr 2011
At least format the code for readability, looking that is a kind of discouraging.
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!