DAEs (ODE 15s)
显示 更早的评论
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 个评论
Oleg Komarov
2011-4-11
format the code
Rafael Freire
2011-4-11
At least format the code for readability, looking that is a kind of discouraging.
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
