DAEs (ODE 15s)
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Rafael Freire
am 11 Apr. 2011
At least format the code for readability, looking that is a kind of discouraging.
Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!